10.计算FPKM和RPKM
有一些五六年前的学生们都成长为了各个生物信息学相关公司的小领导，而且他们都有了自己的公众号，知乎号，也算是一番人物。最近他们跟我反馈面试找不到或者说很难直接考核筛选到认真干活的生信工程师，挺有意思的。让我想起来了早在生信技能树论坛创立之初我为了引流，而规划的200个生信工程师面试题。值得继续分享：
FPKM和RPKM都是表示基因表达量的两个指标，虽然近年来有很多质疑这两种衡量表达量的声音，但有依然有研究采用。下面我们分别来计算下FPKM和RPKM。
FPKM: Fragments Per Kilobase of exon model per Million mapped fragments
计算方式为：
RPKM: Reads Per Kilobase of exon model per Million mapped reads
计算方式为和FPKM的类类似，区别在于RPKM用于单末端测序，而FPKM用于双末端测序。在进行二代测序时，会将所有的DNA打断称为片段，然后再去测序。单末端测序时，一个片段对应一个Read，但是双末端测序时，一个片段会从两端分别测定一次，因此这两个配对的Reads对应的是同一个片段（偶尔也会有一个片段只对应一个Read的情况，另外一个Read因为某些原因被剔除或者丢失了）。
对于FPKM来说，配对到同一片段上的两个Read只会算作一个Read，所以FPKM以Fragment为准，不以Read数为准，其他计算方式完全一样。
matr<-read.table("/Users/express_matrix/48_cell_matrix/48_cell_matrix.txt",header=T)
## 载入程序包
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb<-TxDb.Hsapiens.UCSC.hg38.knownGene
## 定义基因的长度为非冗余exon长度之和
if (F) {
exon_txdb=exons(txdb) #取出人类基因组中的全部外显子
genes_txdb=genes(txdb)
o=findOverlaps(exon_txdb,genes_txdb) # 找出外显子和基因之间的重叠部分
o
t1=exon_txdb[queryHits(o)]# 将与gene重叠的exon提取出来
t2=genes_txdb[subjectHits(o)]
t1=as.data.frame(t1)
t1$geneid=mcols(t2)[,1]
#lapply:历遍列表向量中的每个元素，并且使用指定的函数来对其元素进行处理。返回向量列表。
# 函数split()可以按照分组因子，把向量，矩阵和数据框进行适当的分组；返回值为列表
g_l=lapply(split(t1,t1$geneid), function(x){
head(x)
tmp=apply(x, 1,function(y){
y[2]:y[3]
})
length(unique(unlist(tmp)))
})
head(g_l)
g_l=data.frame(gene_id=names(g_l), length=as.numeric(g_l))
save(g_l,file ='hg38_g_l.RData')
}
load('hg38_g_l.RData')
## 下面是定义基因长度为最长转录本长度
if(F) {
t_l=transcriptLengths(txdb)
head(t_l)
t_l=na.omit(t_l)
t_l=t_l[order(t_l$gene_id, t_l$tx_len, decreasing = T),]
str(t_l)
t_l=t_l[!duplicated(t_l$gene_id),]
head(t_l)
g_l=t_l[,c(3,5)]
}
head(g_l)
library(org.Hs.eg.db)
s2g=toTable(org.Hs.egSYMBOL)
head(s2g)
gl=merge(g_l,s2g,by='gene_id') # merge函数可以实现对两个数据框根据共同的列名来进行连接。
# 参考counts2rpkm, 定义基因长度为非冗余CDs之和
a<-matr[,7:ncol(matr)]
rownames(a)=matr$Geneid
ng=intersect(rownames(a),gl$symbol)
exprset=a[ng,]
lengths=gl[match(ng,gl$symbol),'length']
exprset[1:4,1:4]
total_count<-colSums(exprset)
head(total_count)
head(lengths)
10^9/(1122*121297)
rpkm<-t(do.call(rbind,
lapply(1:length(total_count),
function(i){
10^9*exprset[,i]/lengths/total_count[i]
})))
rownames(rpkm)=rownames(exprset)
colnames(rpkm)=colnames(exprset)
rpkm[1:4,1:4]
cells_48_rpkm=rpkm
cells_48_matr=exprset
library(stringr)
colnames(cells_48_matr)=str_split(colnames(cells_48_matr),'\\.',simplify=T)[,1]
colnames(cells_48_rpkm)=str_split(colnames(cells_48_rpkm),'\\.',simplify=T)[,1]