一、从上游输出文件结果中获取基因有效长度
一般而言，RNA-seq得到原始counts表达矩阵最常用到的上游软件就是featureCounts和Salmon了，在这两类软件的输出结果中，除了基因（或转录本）的counts信息外，也包含了基因有效长度信息，如featureCounts输出文件的Length这一列对应的就是基因有效长度。获取基因有效长度的最简便方法是直接从featureCounts或salmon的输出文件中提取。
featureCounts中基因有效长度Length：即为基因的非冗余外显子长度之和！
featureCounts中基因有效长度Length具体是怎么计算来的？见以下的参考文章：
featureCounts软件中的length是怎样计算的？ https://www.jianshu.com/p/35b52d309d8e
featureCounts得到的counts计算 cpm、 tpm、FPKM，代码如下：
# # 读入featureCounts矩阵
expr_df <- read.table( "merged.featureCounts.txt",header=T, row.names=1, check.names=F, sep="\t")
dim(expr_df); names(expr_df)
head(expr_df[,1:7])
#提取基因信息,featureCounts前几列
featureCounts_meta <- expr_df[,1:5] ;head(featureCounts_meta)
#提取counts
expr_df <- expr_df[,6:ncol(expr_df)]
##基因表达量之和>0
expr_df <- expr_df[rowSums(expr_df)>0,]
head(expr_df[,1:6])
## 保存counts矩阵
write.table(expr_df, "merged.Counts.txt",quote=F, sep="\t", row.names=T, col.names=T )
prefix <-"merged_samples" #设置输出文件前缀名
#----- cPM计算 ------
cpm <- t(t(expr_df)/colSums(expr_df) * 1000000) #参考cpm定义
avg_cpm <- data.frame(avg_cpm=rowMeans(cpm))
# 保存
write.table(avg_cpm, paste0(prefix,"_avg_cpm.xls"), quote=F, sep="\t", row.names=T, col.names=T )
write.table(cpm, paste0(prefix,"_cpm.xls"), quote=F, sep="\t", row.names=T, col.names=T )
# ----- TPM计算 ------
# 基因长度，目标基因的外显子长度之和除以1000，单位是Kb，不是bp
kb <- featureCounts_meta$Length / 1000
rpk <- expr_df / kb #每千碱基reads (“per million” scaling factor) 长度标准化
tpm <- t(t(rpk)/colSums(rpk) * 1000000) # 每百万缩放因子 (“per million” scaling factor ) 深度标准化
avg_tpm <- data.frame(avg_tpm=rowMeans(tpm))
# 保存
write.table(avg_tpm, paste0(prefix,"_avg_tpm.xls"),quote=F, sep="\t", row.names=T, col.names=T )
write.table(tpm, paste0(prefix,"_tpm.xls"), quote=F, sep="\t", row.names=T, col.names=T )
# ----- FPKM计算 ------
fpkm <- t(t(rpk)/colSums(expr_df) * 10^6)
head(fpkm[,1:3])
# 保存
write.table(fpkm,file= paste0(prefix, "_fpkm.xls"), quote=F, sep="\t", row.names=T, col.names=T )
# ----- FPKM转化为TPM ------
fpkm_to_tpm = t(t(fpkm)/colSums(fpkm))*10^6
head(fpkm_to_tpm)
参考
获取基因有效长度的N种方法：https://www.jianshu.com/p/37b5976d39e2
featureCounts软件中的length是怎样计算的？ https://www.jianshu.com/p/35b52d309d8e
若有错误，恳请指出，共同进步！
