2021.04.22【RNA-seq流程】丨count值转换为FPKM值优化2.0

• 解决每次转换需要设置样本数和基因数目
• 实现基因count值与length精准匹配

• 大概半年前，我写过一篇将HTseq生成的基因COUNT值转换为FPKM值文章，用于对count的入门级均一化处理。随着项目越做越多，逐渐发现了之前写的脚本的局限性。比如，每次换算都需要设置包括样品数，基因数目等参数。另外，以前的转换脚本哪怕使用同一个gff文件，定量得到的基因数目和makeTxDbFromGFF中exonsby、transcriptsby定义的基因数目可能是不一样的，这就导致count列和length列不能一一对应进行计算，而产生差错。count如何精准换算成FPKM就成了非常重要的一个问题。因此，我用了2天时间（R确实不太熟悉）研究了多个解决方法， 最后通过匹配方式将gene_count与gene_length通过gene名称匹配起来，完成转换。

• R version 3.6.1
• 依赖包GenomicFeatures

• 最开始的方向是尝试统一基因数目，比如将htseq和makeTxDbFromGFF读取相同的gtf或者gff文件。在发现有不少其他类型的注释后（比如lnc_RNA、cDNA_match、region），又试图抓取exon、mRNA、transcript、CDS这些主要类型，可是基因数目仍然不同（有些exon已经注释到其他类型中，定量的时候无法识别）。于是，在确定之前那个思路已经行不通的时候，就开始想着另外一种方式了，这就是精确匹配。
• 匹配思路其实最初就想到了，但是匹配计算量较大，并且当时测试的注释文件是得到与定量相同的基因数目，可以直接进行计算。随着使用的不同物种注释文件多了，慢慢才发现这个问题。所以现在才进行匹配处理。需要注意的是，匹配前要对基因名称中的特殊字符进行替换，这是R对通配符的识别及替换导致的，因此我们也需要进行替换去匹配，在这里我使用的是merge()函数，关于该函数的使用方法可以参考merge函数匹配数据

#!/usr/bin/R

library(GenomicFeatures)
txdb <- makeTxDbFromGFF("../ref/GCF_000002035.6_GRCz11_genomic.gff",format="auto")
exons_gene <- exonsBy(txdb, by = "gene")
exons_gene_lens <- lapply(exons_gene,function(x){sum(width(reduce(x)))})
n=t(as.data.frame(exons_gene_lens))
write.table(n,'gene_length.txt',col.names=F,row.names=T,quote=F,sep='\t')
length=read.table("gene_length.txt",header=F,sep='\t',check.names=F)
names(length)<-c('gene','length')

#导入所有reads的count值,对特殊符号进行处理
count <- read.table(file="all_count.txt",header = T,sep = "\t")

count[,1] = gsub(':','.',count[,1])
count[,1] = gsub('-','.',count[,1])
count[,1] = gsub(' ','.',count[,1])
count[,1] = gsub('/','.',count[,1])

merge<-merge(count,length,by = 'gene') #匹配gene_count与gene_length
dim(merge)
count <- merge[,1:(dim(merge)[2]-1)]
gene_num <- dim(merge)[1]
sample_num <- dim(merge)[2]-2 #减去gene_name和gene_length列
#从第二列开始，将每个样本的count循环转化成FPKM
i <- 2

repeat{
        mapped_reads <- sum(merge[1:gene_num,i])#计算每个样本的mapped reads数
        FPKM <- merge[1:gene_num,i]/(10^-9*mapped_reads*merge[1:gene_num,dim(merge)[2]])#计算FPKM值
        FPKM <- pmax(FPKM,0)#去掉矫正带来的负值
        count = data.frame(count[1:gene_num,],FPKM)#添加到count表格后面i
        i <- i + 1

        if(i > sample_num+1){
        break
        }
}

#生成表格列名称
head(count)
count_colname <- read.table("all_count.txt",header = F,nrow = 1,as.is=TRUE)
FPKM_colname <- paste(count_colname[1,],"_FPKM",sep="")
FPKM_colname
colname <- c(count_colname,FPKM_colname)
col_name <- colname[-which(colname=="gene_FPKM")]#删除gene_FPKM
col_name
names(count) <- col_name
head(count)

#生成表格
write.csv(count,"ALL_FPKM.csv")
write.table(count[,c(1,(sample_num+2):(sample_num*2+1))],"FPKM_table.txt",row.names = FALSE, quote = FALSE, sep = "\t")
#/usr/bin/bash sed 's/,/\t/g' FPKM_table_tmp.txt > FPKM_table.txt
head(read.csv("ALL_FPKM.csv"))

结果展示

• 首先可以看一下merge（）的结果

• 

• 首先，通过匹配基因名称将length添加到最后一列，前面每个样品的count值只需要和最后一列的值进行计算。

• 其次，通过dim()识别行列，确定了样品数量和基因数目，避免了每次调试。

• 最后，生成出来的gene数目是定量和makeTxDbFromGFF定义的交集。尽管进行了一些特殊符号的换算，但仍然有gene没有识别出来。

总结

• 此次的优化版本解决了匹配计算的问题，当然，在匹配过程中带来了一些新的问题。之后仍然需要在项目中不断测试。

• 欢迎加群交流，遇见二维码过期可添加VX：bbplayer2021 ，备注 申请加入生信交流群。

•