在RNAseq的下游分析中贬丛，一般都會將上游處理完得到的原始counts數(shù)轉(zhuǎn)變?yōu)镕PKM/RPKM或是TPM來進行后續(xù)的展示或分析，其定義和計算公式在我之前的文章中有所總結(jié)Counts FPKM RPKM TPM 的轉(zhuǎn)化 - 簡書 (jianshu.com)。
需要注意一點的是储玫，在計算FPKM/RPKM和TPM時澎羞，基因長度一般指的都是基因有效長度effective length，即該基因的外顯子總長度或轉(zhuǎn)錄本總長度锅锨，以此為標準來消除測序造成的基因長度影響才更為準確螺句。參見生信技能樹文章：基因長度之多少 | 生信菜鳥團 (bio-info-trainee.com)
那么問題來了，在計算FPKM/RPKM時橡类，每個基因的基因有效長度數(shù)據(jù)該如何獲取呢蛇尚？
我總結(jié)了幾種獲取基因有效長度（或非冗余總外顯子長度、總轉(zhuǎn)錄本長度）的方法顾画，現(xiàn)整理如下：
一取劫、從上游輸出文件結(jié)果中獲取基因有效長度
一般而言匆笤，RNA-seq得到原始counts表達矩陣最常用到的上游軟件就是featureCounts和Salmon了，在這兩類軟件的輸出結(jié)果中谱邪，除了基因（或轉(zhuǎn)錄本）的counts信息外炮捧，也包含了基因有效長度信息，如featureCounts輸出文件的Length這一列對應的就是基因有效長度惦银。（P.S. 之前一直以為featureCounts的Length只是單純的基因長度咆课，后來經(jīng)過多種方法比較后發(fā)現(xiàn)其實Length這一列就已經(jīng)是基因的有效長度了...在文章后面我也會展示這幾種方法比較的結(jié)果）
因此，最方便的做法就是在下游獲取counts矩陣時扯俱，將基因有效長度信息也同時提取出來用于后續(xù)的基因表達量轉(zhuǎn)化书蚪。
1. 針對featureCounts的輸出文件
在R中讀取featureCounts的輸出文件，提取Length和對應的geneid信息迅栅，再按照counts中的rowname（geneid）匹配排序殊校，即可進行后續(xù)的TPM、FPKM值的計算了读存。
rm(list=ls())
options(stringsAsFactors = F)
library(tidyverse) # ggplot2 stringer dplyr tidyr readr purrr tibble forcats
library(data.table) #可多核讀取文件
a1 <- fread('counts.txt', header = T, data.table = F)#載入counts为流，第一列設置為列名
### counts矩陣的構(gòu)建
counts <- a1[,7:ncol(a1)] #截取樣本基因表達量的counts部分作為counts
rownames(counts) <- a1$Geneid #將基因名作為行名
### 從featurecounts 原始輸出文件counts.txt中提取Geneid、Length(轉(zhuǎn)錄本長度)让簿，
geneid_efflen <- subset(a1,select = c("Geneid","Length"))
colnames(geneid_efflen) <- c("geneid","efflen")
geneid_efflen_fc <- geneid_efflen #用于之后比較
### 取出counts中g(shù)eneid的對應的efflen
dim(geneid_efflen)
efflen <- geneid_efflen[match(rownames(counts),
geneid_efflen$geneid),
"efflen"]
### 計算 TPM
#TPM (Transcripts Per Kilobase Million) 每千個堿基的轉(zhuǎn)錄每百萬映射讀取的Transcripts
counts2TPM <- function(count=count, efflength=efflen){
RPK <- count/(efflength/1000) #每千堿基reads (“per million” scaling factor) 長度標準化
PMSC_rpk <- sum(RPK)/1e6 #RPK的每百萬縮放因子 (“per million” scaling factor ) 深度標準化
RPK/PMSC_rpk
}
tpm <- as.data.frame(apply(counts,2,counts2TPM))
colSums(tpm)
2. 針對Salmon的輸出文件
Salmon的輸出結(jié)果以及各內(nèi)容的解釋如下敬察。Salmon Output File Formats - Salmon 1.8.0 documentation
值得注意的是，salmon不僅給出了基因有效長度Length（轉(zhuǎn)錄本長度）尔当，還給出了EffectiveLength莲祸，即經(jīng)過考慮各種因素矯正后的轉(zhuǎn)錄本有效長度。官方更推薦使用EffectiveLength進行后續(xù)的分析居凶，它結(jié)果中的TPM值也是根據(jù)EffectiveLength計算的虫给。
我們一般使用tximport導入salmon的輸出文件“quant.sf”（轉(zhuǎn)錄本的統(tǒng)計結(jié)果）和轉(zhuǎn)錄本id與gene symbol對應關(guān)系文件，會生成以下幾個數(shù)據(jù)：
"abundance" "counts" "length" "countsFromAbundance"
tximport生成的Length就是EffectiveLength侠碧，而"abundance" 就是TPM值抹估，我們提取Length用于后續(xù)計算FPKM。注意弄兜，由于每個樣本中基因的EffectiveLength有差異药蜻，我們要提取的實際為EffectiveLength的矩陣(或者可以每行EffectiveLength取均值)。
library(tximport)
#t2s為從gtf文件中提取的transcript_id和symbol的對應關(guān)系文件
t2s <- fread("t2s_vm29_gencode.txt", data.table = F, header = F)
##創(chuàng)建quant.sf所在路徑 導入salmon文件處理匯總
quantdir <- file.path(getwd(),'salmon'); quantdir
files <- list.files(pattern="*quant.sf",quantdir,recursive=T); files #顯示目錄下所有符合要求的文件
files <- file.path(quantdir,files);files
txi_gene <- tximport(files, type = "salmon", tx2gene = t2s)
##提取出counts/tpm表達矩陣
counts <- apply(txi_gene$counts,2,as.integer) #將counts數(shù)取整
rownames(counts) <- rownames(txi_gene$counts)
tpm <- txi_gene$abundance ###abundance為基因的Tpm值
###提取geneid_efflen_mat
geneid_efflen_mat <- txi_gene$length ###Length為基因的轉(zhuǎn)錄本有效長度
## 計算 TPM 替饿、FPKM
if (F) { #可直接從txi的"abundance" 中提取语泽，不用運行
tpm <- data.frame(rownames(counts),row.names = rownames(counts))
for (i in 1:ncol(counts)) {
count <- counts[,i]
efflength <- geneid_efflen_mat[,i]
RPK <- count/(efflength/1000) #每千堿基reads (reads per million) 長度標準化
PMSC_rpk <- sum(RPK)/1e6 #RPK每百萬縮放因子 (“per million” scaling factor ) 深度標準化
tpm00 <- RPK/PMSC_rpk
tpm <- data.frame(tpm,tpm00)
rm(tpm00)
}
tpm <- tpm[,-1]; colnames(tpm) <- colnames(counts); head(tpm)
}
## 計算 fpkm
if(T){
fpkm <- data.frame(rownames(counts),row.names = rownames(counts))
for (i in 1:ncol(counts)) {
count <- counts[,i]
efflength <- geneid_efflen_mat[,i]
PMSC_counts <- sum(count)/1e6 #counts的每百萬縮放因子 (“per million” scaling factor) 深度標準化
FPM <- count/PMSC_counts #每百萬reads/Fragments (Reads/Fragments Per Million) 長度標準化
fpkm00 <- FPM/(efflength/1000)
fpkm <- data.frame(fpkm,fpkm00)
rm(fpkm00)
}
fpkm <- fpkm[,-1]; colnames(fpkm) <- colnames(counts)
}
如果想要提取一般意義上的基因有效長度，需要利用“quant.genes.sf”文件（基因的統(tǒng)計結(jié)果视卢，需要在進行salmon時加上參數(shù) -g 踱卵，后接gtf文件），提取Length這一列的信息。
a2 <- fread("quant.genes.sf",
data.table = F)
geneid_efflen <- subset(a2, select = c("Name","Length"))
colnames(geneid_efflen) <- c("geneid","efflen")
geneid_efflen_salmon <- geneid_efflen #用于之后比較
二惋砂、 從gtf文件中計算獲取基因有效長度
整理了兩種從gtf文件中計算獲取基因有效長度的方法（非冗余外顯子長度之和）妒挎，參考這兩篇文章：
基因長度并不是end-start - 簡書 (jianshu.com)
Htseq Count To Fpkm | KeepNotes blog (bioinfo-scrounger.com)
由于處理數(shù)據(jù)量很大，代碼速度運行較慢西饵，因此在以下代碼中還調(diào)用了parallel包進行多核運算處理
1. 利用GenomicFeatures包導入gtf處理
library(parallel) #并行計算 parApply parLapply parSaplly
cl <- makeCluster(0.75*detectCores()) #設計啟用計算機3/4的核
## 利用GenomicFeatures包導入gtf處理
library(GenomicFeatures)
txdb <- makeTxDbFromGFF("gencode.vM25.chr_patch_hapl_scaff.annotation.gtf.gz",
format="gtf")
exons_gene <- exonsBy(txdb, by = "gene") ###提取基因外顯子
head(exons_gene)
##計算總外顯子長度：用reduce去除掉重疊冗余的部分酝掩，,width統(tǒng)計長度，最后計算總長度
exons_gene_lens <- parLapply(cl,exons_gene,function(x){sum(width(reduce(x)))})
exons_gene_lens[1:10]
##轉(zhuǎn)換為dataframe
geneid_efflen <- data.frame(geneid=names(exons_gene_lens),
efflen=as.numeric(exons_gene_lens))
geneid_efflen_gtf1 <- geneid_efflen
2.利用rtracklayer包導入gtf處理
##利用rtracklayer包import導入處理
gtf <- as.data.frame(rtracklayer::import("gencode.vM25.chr_patch_hapl_scaff.annotation.gtf.gz"))
table(gtf$type)
exon <- gtf[gtf$type=="exon",
c("start","end","gene_id")]
exon_bygeneid <- split(exon,exon$gene_id) #按照每個geneid所含的exon排序成列表
efflen <- parLapply(cl,exon_bygeneid,function(x){
tmp <- apply(x,1,function(y){ y[1]:y[2] }) #輸出exon長度值的所有元素
length(unique(unlist(tmp))) #去重復并統(tǒng)計exon長度元素的數(shù)量
})
##轉(zhuǎn)換為dataframe
geneid_efflen <- data.frame(geneid=names(efflen),
efflen=as.numeric(efflen))
geneid_efflen_gtf2 <- geneid_efflen
所得結(jié)果的比較
現(xiàn)在就可以來進行基因有效長度之間的比較啦眷柔。
首先看看從gtf文件中獲取基因有效長度的兩種方法是否有差異期虾。可以發(fā)現(xiàn)驯嘱，僅有極少數(shù)efflen有差異镶苞，因此這兩種方法可以說是幾乎沒什么差別了：
再比較一下geneid_efflen_gtf1和geneid_efflen_salmon，發(fā)現(xiàn)有一半的efflen是不匹配的宙拉？仔細一想這也是可以理解的宾尚，因為上游中salmon是對樣本中的轉(zhuǎn)錄本進行的統(tǒng)計丙笋，這說明了樣本中有一半的基因并未表達其全部的轉(zhuǎn)錄本而已：
再將geneid_efflen_gtf1和geneid_efflen_fc進行比較谢澈，發(fā)現(xiàn)全都能匹配上，這說明featureCounts的Length確實就已經(jīng)是我們想要的基因有效長度了（即非冗余外顯子總長度）：
總結(jié)：
- 獲取基因有效長度的最簡便方法是直接從featureCounts或salmon的輸出文件中提取御板。
但需要注意的是锥忿，featureCounts中基因有效長度Length即為基因的非冗余外顯子總長度，而salmon中的基因有效長度Length是目標基因的轉(zhuǎn)錄本總長度怠肋，由于樣本中只有部分基因會表達其全部類型的轉(zhuǎn)錄本敬鬓，因此salmon中的轉(zhuǎn)錄本總長度會有部分小于非冗余外顯子總長度。- salmon輸出結(jié)果中不僅給出了基因有效長度Length（轉(zhuǎn)錄本總長度）笙各，還給出了經(jīng)過考慮各種因素矯正后的轉(zhuǎn)錄本有效長度EffectiveLength。Salmon官方更推薦使用EffectiveLength進行后續(xù)的分析，認為其能更好消除測序時基因長度的影響汉额，它結(jié)果中的TPM值也是根據(jù)EffectiveLength計算的庆捺，后續(xù)分析中可以直接采用。
- 在沒有上游原始輸出文件的情況下惶楼，也可以采取直接從gtf文件中計算的方法右蹦，獲取每個基因的非冗余外顯子總長度得到基因有效長度。還可以保存geneid_efflen便于之后再讀燃呔琛：
write.csv(geneid_efflen,file = "geneid_efflen_vm25_gencode.csv",row.names = F)
參考資料
基因長度之多少 | 生信菜鳥團 (bio-info-trainee.com)
Counts FPKM RPKM TPM 的轉(zhuǎn)化 - 簡書 (jianshu.com)
基因長度并不是end-start - 簡書 (jianshu.com)
Htseq Count To Fpkm | KeepNotes blog (bioinfo-scrounger.com)
Salmon Output File Formats - Salmon 1.8.0 documentation
