一、从上游输出文件结果中获取基因有效长度

一般而言，RNA-seq得到原始counts表达矩阵最常用到的上游软件就是featureCounts和Salmon了，在这两类软件的输出结果中，除了基因（或转录本）的counts信息外，也包含了基因有效长度信息，如featureCounts输出文件的Length这一列对应的就是基因有效长度。获取基因有效长度的最简便方法是直接从featureCounts或salmon的输出文件中提取。

featureCounts中基因有效长度Length：即为基因的非冗余外显子长度之和！

featureCounts中基因有效长度Length具体是怎么计算来的？见以下的参考文章：

featureCounts软件中的length是怎样计算的？ //www.greatytc.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种方法：//www.greatytc.com/p/37b5976d39e2
featureCounts软件中的length是怎样计算的？ //www.greatytc.com/p/35b52d309d8e

若有错误，恳请指出，共同进步！