在RNAseq的下游分析中，一般都會將上游處理完得到的原始counts數(shù)轉(zhuǎn)變?yōu)镕PKM/RPKM或是TPM來進(jìn)行后續(xù)的展示或分析，其定義和計(jì)算公式在我之前的文章中有所總結(jié)Counts FPKM RPKM TPM 的轉(zhuǎn)化 - 簡書 (jianshu.com)。
需要注意一點(diǎn)的是，在計(jì)算FPKM/RPKM和TPM時(shí)，基因長度一般指的都是基因有效長度effective length，即該基因的外顯子總長度或轉(zhuǎn)錄本總長度，以此為標(biāo)準(zhǔn)來消除測序造成的基因長度影響才更為準(zhǔn)確。參見生信技能樹文章：基因長度之多少 | 生信菜鳥團(tuán) (bio-info-trainee.com)

那么問題來了，在計(jì)算FPKM/RPKM時(shí)，每個(gè)基因的基因有效長度數(shù)據(jù)該如何獲取呢？
我總結(jié)了幾種獲取基因有效長度（或非冗余總外顯子長度、總轉(zhuǎn)錄本長度）的方法，現(xiàn)整理如下：

一、從上游輸出文件結(jié)果中獲取基因有效長度

一般而言，RNA-seq得到原始counts表達(dá)矩陣最常用到的上游軟件就是featureCounts和Salmon了，在這兩類軟件的輸出結(jié)果中，除了基因（或轉(zhuǎn)錄本）的counts信息外，也包含了基因有效長度信息，如featureCounts輸出文件的Length這一列對應(yīng)的就是基因有效長度。（P.S. 之前一直以為featureCounts的Length只是單純的基因長度，后來經(jīng)過多種方法比較后發(fā)現(xiàn)其實(shí)Length這一列就已經(jīng)是基因的有效長度了...在文章后面我也會展示這幾種方法比較的結(jié)果）

因此，最方便的做法就是在下游獲取counts矩陣時(shí)，將基因有效長度信息也同時(shí)提取出來用于后續(xù)的基因表達(dá)量轉(zhuǎn)化。

1. 針對featureCounts的輸出文件

在R中讀取featureCounts的輸出文件，提取Length和對應(yīng)的geneid信息，再按照counts中的rowname（geneid）匹配排序，即可進(jìn)行后續(xù)的TPM、FPKM值的計(jì)算了。

featurecounts輸出文件格式

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，第一列設(shè)置為列名

### counts矩陣的構(gòu)建
counts <- a1[,7:ncol(a1)] #截取樣本基因表達(dá)量的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的對應(yīng)的efflen
dim(geneid_efflen)
efflen <- geneid_efflen[match(rownames(counts),
                              geneid_efflen$geneid),
                        "efflen"]

### 計(jì)算 TPM
#TPM (Transcripts Per Kilobase Million)  每千個(gè)堿基的轉(zhuǎn)錄每百萬映射讀取的Transcripts
counts2TPM <- function(count=count, efflength=efflen){
  RPK <- count/(efflength/1000)       #每千堿基reads (“per million” scaling factor) 長度標(biāo)準(zhǔn)化
  PMSC_rpk <- sum(RPK)/1e6        #RPK的每百萬縮放因子 (“per million” scaling factor ) 深度標(biāo)準(zhǔn)化
  RPK/PMSC_rpk                    
  }  
tpm <- as.data.frame(apply(counts,2,counts2TPM))
colSums(tpm)

其中g(shù)eneid_efflen內(nèi)容如下

geneid_efflen

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進(jìn)行后續(xù)的分析，它結(jié)果中的TPM值也是根據(jù)EffectiveLength計(jì)算的。

Salmon的輸出結(jié)果

Salmon的輸出結(jié)果官方解釋

我們一般使用tximport導(dǎo)入salmon的輸出文件“quant.sf”（轉(zhuǎn)錄本的統(tǒng)計(jì)結(jié)果）和轉(zhuǎn)錄本id與gene symbol對應(yīng)關(guān)系文件，會生成以下幾個(gè)數(shù)據(jù)：
"abundance" "counts" "length" "countsFromAbundance"
tximport生成的Length就是EffectiveLength，而"abundance" 就是TPM值，我們提取Length用于后續(xù)計(jì)算FPKM。注意，由于每個(gè)樣本中基因的EffectiveLength有差異，我們要提取的實(shí)際為EffectiveLength的矩陣(或者可以每行EffectiveLength取均值)。

library(tximport) 

#t2s為從gtf文件中提取的transcript_id和symbol的對應(yīng)關(guān)系文件
t2s <- fread("t2s_vm29_gencode.txt", data.table = F, header = F) 

##創(chuàng)建quant.sf所在路徑  導(dǎo)入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表達(dá)矩陣
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)錄本有效長度

## 計(jì)算 TPM 、FPKM
  if (F) { #可直接從txi的"abundance"  中提取，不用運(yùn)行
    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) 長度標(biāo)準(zhǔn)化
      PMSC_rpk <- sum(RPK)/1e6        #RPK每百萬縮放因子 (“per million” scaling factor ) 深度標(biāo)準(zhǔn)化
      tpm00 <- RPK/PMSC_rpk  
      tpm <- data.frame(tpm,tpm00)
      rm(tpm00)
    }
    tpm <- tpm[,-1];  colnames(tpm) <- colnames(counts);  head(tpm)
   
  }
  
  ## 計(jì)算 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) 深度標(biāo)準(zhǔn)化
      FPM <- count/PMSC_counts        #每百萬reads/Fragments (Reads/Fragments Per Million) 長度標(biāo)準(zhǔn)化
      fpkm00 <- FPM/(efflength/1000)
      fpkm <- data.frame(fpkm,fpkm00)
      rm(fpkm00)
    }
    fpkm <- fpkm[,-1];  colnames(fpkm) <- colnames(counts)
  }

如果想要提取一般意義上的基因有效長度，需要利用“quant.genes.sf”文件（基因的統(tǒng)計(jì)結(jié)果，需要在進(jìn)行salmon時(shí)加上參數(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文件中計(jì)算獲取基因有效長度

整理了兩種從gtf文件中計(jì)算獲取基因有效長度的方法（非冗余外顯子長度之和），參考這兩篇文章：
基因長度并不是end-start - 簡書 (jianshu.com)
Htseq Count To Fpkm | KeepNotes blog (bioinfo-scrounger.com)
由于處理數(shù)據(jù)量很大，代碼速度運(yùn)行較慢，因此在以下代碼中還調(diào)用了parallel包進(jìn)行多核運(yùn)算處理

1. 利用GenomicFeatures包導(dǎo)入gtf處理

library(parallel) #并行計(jì)算  parApply parLapply parSaplly 
cl <- makeCluster(0.75*detectCores())  #設(shè)計(jì)啟用計(jì)算機(jī)3/4的核

## 利用GenomicFeatures包導(dǎo)入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)

    ##計(jì)算總外顯子長度：用reduce去除掉重疊冗余的部分，,width統(tǒng)計(jì)長度，最后計(jì)算總長度
    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包導(dǎo)入gtf處理

##利用rtracklayer包import導(dǎo)入處理 
    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)   #按照每個(gè)geneid所含的exon排序成列表
    
    efflen <- parLapply(cl,exon_bygeneid,function(x){
      tmp <- apply(x,1,function(y){  y[1]:y[2]  }) #輸出exon長度值的所有元素            
      length(unique(unlist(tmp))) #去重復(fù)并統(tǒng)計(jì)exon長度元素的數(shù)量
      }) 

    ##轉(zhuǎn)換為dataframe
    geneid_efflen <- data.frame(geneid=names(efflen),
                                efflen=as.numeric(efflen))
    geneid_efflen_gtf2 <- geneid_efflen

所得結(jié)果的比較

現(xiàn)在就可以來進(jìn)行基因有效長度之間的比較啦。
首先看看從gtf文件中獲取基因有效長度的兩種方法是否有差異。可以發(fā)現(xiàn)，僅有極少數(shù)efflen有差異，因此這兩種方法可以說是幾乎沒什么差別了：

從gtf文件中獲取efflen的比較

再比較一下geneid_efflen_gtf1和geneid_efflen_salmon，發(fā)現(xiàn)有一半的efflen是不匹配的？仔細(xì)一想這也是可以理解的，因?yàn)樯嫌沃衧almon是對樣本中的轉(zhuǎn)錄本進(jìn)行的統(tǒng)計(jì)，這說明了樣本中有一半的基因并未表達(dá)其全部的轉(zhuǎn)錄本而已：

salmon和gtf中獲取的efflen比較

再將geneid_efflen_gtf1和geneid_efflen_fc進(jìn)行比較，發(fā)現(xiàn)全都能匹配上，這說明featureCounts的Length確實(shí)就已經(jīng)是我們想要的基因有效長度了（即非冗余外顯子總長度）：

featureCounts和gtf中獲取的efflen比較

總結(jié)：

獲取基因有效長度的最簡便方法是直接從featureCounts或salmon的輸出文件中提取。
但需要注意的是，featureCounts中基因有效長度Length即為基因的非冗余外顯子總長度，而salmon中的基因有效長度Length是目標(biāo)基因的轉(zhuǎn)錄本總長度，由于樣本中只有部分基因會表達(dá)其全部類型的轉(zhuǎn)錄本，因此salmon中的轉(zhuǎn)錄本總長度會有部分小于非冗余外顯子總長度。

salmon輸出結(jié)果中不僅給出了基因有效長度Length（轉(zhuǎn)錄本總長度），還給出了經(jīng)過考慮各種因素矯正后的轉(zhuǎn)錄本有效長度EffectiveLength。Salmon官方更推薦使用EffectiveLength進(jìn)行后續(xù)的分析，認(rèn)為其能更好消除測序時(shí)基因長度的影響，它結(jié)果中的TPM值也是根據(jù)EffectiveLength計(jì)算的，后續(xù)分析中可以直接采用。

在沒有上游原始輸出文件的情況下，也可以采取直接從gtf文件中計(jì)算的方法，獲取每個(gè)基因的非冗余外顯子總長度得到基因有效長度。還可以保存geneid_efflen便于之后再讀取：
write.csv(geneid_efflen,file = "geneid_efflen_vm25_gencode.csv",row.names = F)

參考資料
基因長度之多少 | 生信菜鳥團(tuán) (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

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av

獲取基因有效長度的N種方法

獲取基因有效長度的N種方法

一、從上游輸出文件結(jié)果中獲取基因有效長度

1. 針對featureCounts的輸出文件

2. 針對Salmon的輸出文件

二、從gtf文件中計(jì)算獲取基因有效長度

1. 利用GenomicFeatures包導(dǎo)入gtf處理

2.利用rtracklayer包導(dǎo)入gtf處理

所得結(jié)果的比較

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九 欧美,1769亚洲,黄色成人av

獲取基因有效長度的N種方法

一、從上游輸出文件結(jié)果中獲取基因有效長度

1. 針對featureCounts的輸出文件

2. 針對Salmon的輸出文件

二、 從gtf文件中計(jì)算獲取基因有效長度

1. 利用GenomicFeatures包導(dǎo)入gtf處理

2.利用rtracklayer包導(dǎo)入gtf處理

所得結(jié)果的比較

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av

二、從gtf文件中計(jì)算獲取基因有效長度