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

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)

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) 
} 

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 #用于之后比较

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

##利用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