RPKM概念和三种计算方法
WGCNA的输入到底是什么格式

基因长度

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

最后，感谢生信技能树的分享！