【原】获取基因有效长度的N种方法

在RNAseq的下游分析中，一般都会将上游处理完得到的原始counts数转变为FPKM/RPKM或是TPM来进行后续的展示或分析，其定义和计算公式在 前面的分享是：Counts FPKM RPKM TPM CPM 的转化 提到了。

需要注意一点的是，在计算FPKM/RPKM和TPM时，基因长度一般指的都是基因有效长度effective length，即该基因的外显子总长度或转录本总长度，以此为标准来消除测序造成的基因长度影响才更为准确。参见生信技能树文章：

基因长度之多少 | 生信菜鸟团 (bio-info-trainee.com)

那么问题来了，在计算FPKM/RPKM时，每个基因的基因有效长度数据该如何获取呢？
我总结了几种获取基因有效长度（或非冗余总外显子长度、总转录本长度）的方法，现整理如下：

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

一般而言，RNA-seq得到原始counts表达矩阵最常用到的上游软件就是featureCounts和Salmon了，在这两类软件的输出结果中，除了基因（或转录本）的counts信息外，也包含了基因有效长度信息，如featureCounts输出文件的Length这一列对应的就是基因有效长度。（P.S. 之前一直以为featureCounts的Length只是单纯的基因长度，后来经过多种方法比较后发现其实Length这一列就已经是基因的有效长度了...在文章后面我也会展示这几种方法比较的结果）

因此，最方便的做法就是在下游获取counts矩阵时，将基因有效长度信息也同时提取出来用于后续的基因表达量转化。

1. 针对featureCounts的输出文件

在R中读取featureCounts的输出文件，提取Length和对应的geneid信息，再按照counts中的rowname（geneid）匹配排序，即可进行后续的TPM、FPKM值的计算了。

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，第一列设置为列名 ### counts矩阵的构建 counts <- a1[,7:ncol(a1)] #截取样本基因表达量的counts部分作为counts  rownames(counts) <- a1$Geneid #将基因名作为行名 ### 从featurecounts 原始输出文件counts.txt中提取Geneid、Length(转录本长度)， geneid_efflen <- subset(a1,select = c("Geneid","Length"))        colnames(geneid_efflen) <- c("geneid","efflen")   geneid_efflen_fc <- geneid_efflen #用于之后比较 ### 取出counts中geneid的对应的efflen dim(geneid_efflen) efflen <- geneid_efflen[match(rownames(counts),                                                             geneid_efflen$geneid),"efflen"] ### 计算 TPM #TPM (Transcripts Per Kilobase Million)  每千个碱基的转录每百万映射读取的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)

其中geneid_efflen内容如下

geneid_efflen

2. 针对Salmon的输出文件

Salmon的输出结果以及各内容的解释如下。Salmon Output File Formats - Salmon 1.8.0 documentation
值得注意的是，salmon不仅给出了基因有效长度Length（转录本长度），还给出了EffectiveLength，即经过考虑各种因素矫正后的转录本有效长度。官方更推荐使用EffectiveLength进行后续的分析，它结果中的TPM值也是根据EffectiveLength计算的。

Salmon的输出结果

Salmon的输出结果官方解释

我们一般使用tximport导入salmon的输出文件“quant.sf”（转录本的统计结果）和转录本id与gene symbol对应关系文件，会生成以下几个数据：
"abundance" "counts" "length" "countsFromAbundance"
tximport生成的Length就是EffectiveLength，而"abundance" 就是TPM值，我们提取Length用于后续计算FPKM。注意，由于每个样本中基因的EffectiveLength有差异，我们要提取的实际为EffectiveLength的矩阵(或者可以每行EffectiveLength取均值)。

library(tximport)  #t2s为从gtf文件中提取的transcript_id和symbol的对应关系文件 t2s <- fread("t2s_vm29_gencode.txt", data.table = F, header = F) 
##创建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数取整 rownames(counts) <- rownames(txi_gene$counts)tpm <- txi_gene$abundance  ###abundance为基因的Tpm值 
###提取geneid_efflen_mat geneid_efflen_mat <- txi_gene$length  ###Length为基因的转录本有效长度 
## 计算 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”文件（基因的统计结果，需要在进行salmon时加上参数 -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)
由于处理数据量很大，代码速度运行较慢，因此在以下代码中还调用了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统计长度，最后计算总长度    exons_gene_lens <- parLapply(cl,exons_gene,function(x){sum(width(reduce(x)))})     exons_gene_lens[1:10]        ##转换为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))) #去重复并统计exon长度元素的数量      }) 
    ##转换为dataframe    geneid_efflen <- data.frame(geneid=names(efflen),                                efflen=as.numeric(efflen))    geneid_efflen_gtf2 <- geneid_efflen

所得结果的比较

现在就可以来进行基因有效长度之间的比较啦。
首先看看从gtf文件中获取基因有效长度的两种方法是否有差异。可以发现，仅有极少数efflen有差异，因此这两种方法可以说是几乎没什么差别了：

从gtf文件中获取efflen的比较

再比较一下geneid_efflen_gtf1和geneid_efflen_salmon，发现有一半的efflen是不匹配的？仔细一想这也是可以理解的，因为上游中salmon是对样本中的转录本进行的统计，这说明了样本中有一半的基因并未表达其全部的转录本而已：

salmon和gtf中获取的efflen比较

再将geneid_efflen_gtf1和geneid_efflen_fc进行比较，发现全都能匹配上，这说明featureCounts的Length确实就已经是我们想要的基因有效长度了（即非冗余外显子总长度）：

featureCounts和gtf中获取的efflen比较

总结：

1. 获取基因有效长度的最简便方法是直接从featureCounts或salmon的输出文件中提取。
   但需要注意的是，featureCounts中基因有效长度Length即为基因的非冗余外显子总长度，而salmon中的基因有效长度Length是目标基因的转录本总长度，由于样本中只有部分基因会表达其全部类型的转录本，因此salmon中的转录本总长度会有部分小于非冗余外显子总长度。

2. salmon输出结果中不仅给出了基因有效长度Length（转录本总长度），还给出了经过考虑各种因素矫正后的转录本有效长度EffectiveLength。Salmon官方更推荐使用EffectiveLength进行后续的分析，认为其能更好消除测序时基因长度的影响，它结果中的TPM值也是根据EffectiveLength计算的，后续分析中可以直接采用。

3. 在没有上游原始输出文件的情况下，也可以采取直接从gtf文件中计算的方法，获取每个基因的非冗余外显子总长度得到基因有效长度。还可以保存geneid_efflen便于之后再读取：
   write.csv(geneid_efflen,file = "geneid_efflen_vm25_gencode.csv",row.names = F)

参考资料

基因长度之多少 | 生信菜鸟团 (bio-info-trainee.com)

Counts FPKM RPKM TPM 的转化 - 简书 (jianshu.com)

基因长度并不是end-start - 简书 (jianshu.com)

Htseq Count To Fpkm | KeepNotes blog (bioinfo-scrounger.com)

Salmon Output File Formats - Salmon 1.8.0 documentation

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：

• 数据挖掘（GEO,TCGA,单细胞）2022年6月场，快速了解一些生物信息学应用图表
• 生信入门课-2022年6月场，你的生物信息学第一课