一、從上游輸出文件結(jié)果中獲取基因有效長度
一般而言，RNA-seq得到原始counts表達(dá)矩陣最常用到的上游軟件就是featureCounts和Salmon了，在這兩類軟件的輸出結(jié)果中，除了基因（或轉(zhuǎn)錄本）的counts信息外，也包含了基因有效長度信息，如featureCounts輸出文件的Length這一列對應(yīng)的就是基因有效長度。獲取基因有效長度的最簡便方法是直接從featureCounts或salmon的輸出文件中提取。
featureCounts中基因有效長度Length：即為基因的非冗余外顯子長度之和！
featureCounts中基因有效長度Length具體是怎么計算來的？見以下的參考文章：
featureCounts軟件中的length是怎樣計算的？ http://www.itdecent.cn/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)]
##基因表達(dá)量之和>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" #設(shè)置輸出文件前綴名
#----- 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計算 ------
# 基因長度，目標(biāo)基因的外顯子長度之和除以1000，單位是Kb，不是bp
kb <- featureCounts_meta$Length / 1000
rpk <- expr_df / kb #每千堿基reads (“per million” scaling factor) 長度標(biāo)準(zhǔn)化
tpm <- t(t(rpk)/colSums(rpk) * 1000000) # 每百萬縮放因子 (“per million” scaling factor ) 深度標(biāo)準(zhǔn)化
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轉(zhuǎn)化為TPM ------
fpkm_to_tpm = t(t(fpkm)/colSums(fpkm))*10^6
head(fpkm_to_tpm)
參考
獲取基因有效長度的N種方法：http://www.itdecent.cn/p/37b5976d39e2
featureCounts軟件中的length是怎樣計算的？ http://www.itdecent.cn/p/35b52d309d8e
若有錯誤，懇請指出，共同進(jìn)步！
