基因长度
image.png
RPKM = counts数 /（基因长度×文库大小）
（所以测到的counts数和文库大小对于RPKM的影响较大）
1 counts矩阵
library(limma)
rt = read.table("gene matrix.txt",header=T,sep="\t",check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
a=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
a=avereps(a) #对于重复行取平均值
a=a[rowSums(a)>1,] #保证基因至少在2个样本中表达
a<-as.data.frame(a)
a[1:4,1:4]
a
2 基因长度
2.1 定义基因长度为非冗余exon长度之和
TxDb.Hsapiens.UCSC.hg19.knownGene would be a TxDb object, of Homo sapiens data from UCSC build hg19 based on the knownGene Track.
#小鼠
library("TxDb.Mmusculus.UCSC.mm10.knownGene") #该函储存转录本信息，包括基因的外显子信息
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
#人
library("TxDb.Hsapiens.UCSC.hg19.knownGene")
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
if(F){
exon_txdb=exons(txdb) #包含了所有外显子的坐标和外显子编码id
exon_txdb
genes_txdb=genes(txdb) #包含了所有基因的坐标和基因ID
genes_txdb
o = findOverlaps(exon_txdb,genes_txdb) #找到二者之间的交叉
o #找到exon的第18个和gene的第6909个对应
# queryHits subjectHits
# <integer> <integer>
# [1] 18 6909
# [2] 19 6909
# [3] 20 6909
# [4] 21 6909
t1=exon_txdb[queryHits(o)] #将所有对应上的外显子信息提取出来
t2=genes_txdb[subjectHits(o)] #将所有基因信息提取出来
t1=as.data.frame(t1)
t1$geneid=mcols(t2)[,1] #mcols就是取t2的gene_id,给ti加上一列gene_id
#lapply : 遍历列表向量内的每个元素，并且使用指定函数来对其元素进行处理。返回列表向量。
#函数split()可以按照分组因子，把向量，矩阵和数据框进行适当的分组;
#它的返回值是一个列表，代表分组变量每个水平的观测。
g_l = lapply(split(t1,t1$geneid),function(x){
# x=split(t1,t1$geneid)[[1]]
head(x)
tmp=apply(x,1,function(y){
y[2]:y[3]
}) #按行将一个基因中所有外显子的坐标点全部展示出来
length(unique(unlist(tmp))) #对所有坐标点进行去重复的操作，去除外显子坐标存在重合的问题
# sum(x[,4])
})
head(g_l)
g_l=data.frame(gene_id=names(g_l),length=as.numeric(g_l))
save(g_l,file = 'step7-g_l.Rdata')
}
load(file = 'step7-g_l.Rdata') #基因编号和对应外显子的长度
head(g_l)
library(org.Mm.eg.db) #包含gene_id的symbol 小鼠
s2g=toTable(org.Mm.egSYMBOL)
head(s2g)
g_l=merge(g_l,s2g,by='gene_id') #把g_l,s2g两个数据框以'gene_id'为连接进行拼接
gene_id对应的基因名的外显子长度
2.2 定义基因长度为最长转录本长度
#小鼠
library("TxDb.Mmusculus.UCSC.mm10.knownGene") #该函储存转录本信息，包括基因的外显子信息
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
#人
library("TxDb.Hsapiens.UCSC.hg19.knownGene")
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
if(F){
t_l=transcriptLengths(txdb) #transcriptLengths提取基因对应转录本长度
head(t_l)
t_l=na.omit(t_l)
head(t_l)
t_l=t_l[order(t_l$gene_id,t_l$tx_len,decreasing = T),] #按照降序，先给gene_id排序，再给转录本长度排序
head(t_l)
t_l=t_l[!duplicated(t_l$gene_id),]#去重，只保存第一次出现的gene_id，刚好就保留了最长转录本
head(t_l)
g_l=t_l[,c(3,5)] #得到的就是gene_id和对应的最长转录本
}
head(g_l)
#小鼠
library(org.Mm.eg.db) #包含gene_id的symbol
s2g=toTable(org.Mm.egSYMBOL)
#人
#library(org.Hs.eg.db) #包含gene_id的symbol
#s2g=toTable(org.Hs.egSYMBOL)
head(s2g)
g_l=merge(g_l,s2g,by='gene_id') #把g_l,s2g两个数据框以'gene_id'为连接进行拼接
最长转录本长度
3 计算RPKM
#基因长度用上述的最长转录本长度
a[1:4,1:4]
#取a数据框的行名与g_l数据框的symbol列的交集，因为有基因长度的，我们才能计算RPKM
ng=intersect(rownames(a),g_l$symbol)
# 有了counts矩阵和对应的基因长度信息，就很容易进行各种计算了：
exprSet=a[ng,] #保留有长度信息的基因
lengths=g_l[match(ng,g_l$symbol),2] #按照ng的顺序，提取了对应转录本长度
head(lengths)
head(rownames(exprSet))
exprSet[1:4,1:4]
total_count<- colSums(exprSet) #每个样本中所有基因counts数求和是文库大小
head(total_count)
head(lengths)
#计算第四个样本中第一个基因的RPKM值,因为K M，所有最后*10^9。
#total_count[4]
#lengths[1]
#1*10^9/(2434*123319) #RPKM就是counts数/（基因长度×文库大小）
#批量计算RPKM值
rpkm <- t(do.call( rbind,
lapply(1:length(total_count),
function(i){
10^9*exprSet[,i]/lengths/total_count[i]
}) ))
#do.call用法
#dat <- list(matrix(1:25, ncol = 5), matrix(4:28, ncol = 5), matrix(21:45, ncol=5))
#dat_bind <- do.call(cbind,dat)
#dat_bind <- do.call(rbind,dat)
rpkm[1:4,1:4]
rownames(rpkm)<-rownames(exprSet)
colnames(rpkm)<-colnames(exprSet)
rpkm<-as.data.frame(rpkm)
write.table(rpkm,file="RPKM.txt",quote=F,sep = "\t")
rpkm
最后，感谢生信技能树的分享！
网友评论