非R包的ballgown提取FPKM的方法及R中dplyr包和Reduce函数的使用实例

 现在网上查到的分享好多都是用的R中ballgown包进行差异表达分析，有工具就用可是好文明。不过似乎不用这个R包也可以提取FPKM进行分析，虽然麻烦了点，不过也是一个练手的好机会。

 一年后我才发现有个最最最直接简单的方法，stringtie本身提供了-A参数，加上去就可自动生成带有FPKM和TPM的表格。这个故事告诉我们读工具的manual是多么的重要！（2022.3.28）
 统合stringtie生成的gene_abund.tab可以参考这条：

library(dplyr)
path=getwd()
setwd(path)
folders=dir(pattern="gene_abund.tab")
files_path=list()
for(j in 1:length(folders)){
	files_path[j]=paste(path,"/",folders[j],sep="")
}
C=list()
for(i in 1:length(files_path)){
	assign(paste("TPM",i,sep=""),read.table(file = files_path[[i]],header = T,sep = '\t')%>%select(,one_of("Gene.ID","TPM")))
	C[[i]]=get(paste("TPM",i,sep=""))
}
PRJ=Reduce(function(x,y) dplyr::left_join(x,y,by="Gene.ID"),C)
rownames(PRJ)=PRJ[,1]
PRJ=PRJ[,-1]
names(PRJ)=c(folders)
setwd(path)
write.table(PRJ,file="mergeTPM.tsv")

 最好还是按上面的通过stringtie来生成FPKM、TPM，下面这种直接从gtf里面提取的方法并不适用于大多数情况，只有当转录本和基因的名字完全相同时才可以使用，而基本每个基因都对应多个转录本，需要具体看注释情况。

library(dplyr)

path=getwd()
#在/PRJNA*/*/ballgown下执行，即stringtie中-o参数设置的路径下。

setwd(path)
folders=dir(pattern="SRR")
#获取工作目录下名字含有SRR的目录和文件，我这里只有目录。

files_path=list()
for(j in 1:length(folders)){
        files_path[j]=paste(path,"/",folders[j],sep="")
}
#获取各样品目录的绝对路径。

C=list()
for(i in 1:length(files_path)){
        setwd(files_path[[i]])
        assign(paste("t_data",i,sep=""),read.table(file = "t_data.ctab",header = T))
        assign(paste("t_data",i,sep=""),dplyr::select(data.frame(get(paste("t_data",i,sep=""))),9,12))
        C[[i]]=get(paste("t_data",i,sep=""))
}
#循环读取ballgown中各样品目录下的t_data.ctab文件（其中第九列为gene_id,第十二列为FPKM）。
#读取完t_data.ctab文件后利用dplyr包中的select提取第九和十二列并覆盖。
#最后将数据逐个导入列表C中，为下面Reduce循环做准备。

PRJ=Reduce(function(x,y) dplyr::left_join(x,y,by="gene_id"),C)
#最精华的一条代码，通过Reduce实现多个数据框的合并。
#Reduce的格式为：Reduce(f, x, init, right = FALSE, accumulate = FALSE)

rownames(PRJ)=PRJ[,1]
PRJ=PRJ[,-1]
names(PRJ)=c(folders)
#将第一列（序号）作为行名。

setwd(path)
write.table(PRJ,file="mergefpkm.tsv")
#导出回path路径下。

fpkm=read.table(file="mergefpkm.tsv",header=T,row.names=1)
fpkm_log=log2(fpkm+1)
corr=cor(fpkm_log,method="pearson")
write.table(fpkm_log,file="mergefpkm_log.tsv")
write.table(corr,file="cor_fpkm.tsv")
#利用fpkm进行的相关性分析。

 Reduce这个东西感觉还是蛮有用的，可以将一个只能处理两个数据的函数应用于有三个及以上待处理数据的场景，简单的比如三个数累加，复杂一点就是merge合并多个数据框等等，可惜就是有点难懂，下面是官方的help：

 Reduce(f, x, init, right = FALSE, accumulate = FALSE)

 f：a function of the appropriate arity (binary for Reduce, unary for Filter, Find and Position, k-ary for Map if this is called with k arguments). An arbitrary predicate function for Negate.

 x：a vector.

 init：an R object of the same kind as the elements of x.

 right：a logical indicating whether to proceed from left to right (default) or from right to left.

 accumulate：a logical indicating whether the successive reduce combinations should be accumulated. By default, only the final combination is used.

 下面是我自己的理解：
 f指的是要用的函数，可以是简单的“+”法，也可以是复杂的自定义函数，同时也需要指定需要的数据个数；
 x指的是需要导入给上面的f函数进行处理的一系列数据，一般用向量或列表的形式导入；
 init指的是x的初始值？？
 right是否从右到左？？
 accumula是否显示计算过程？？

 似乎这个f指定的格式有好几种，我感觉还是我例子里用的那种好理解一点点。。。

 Reduce(function(x,y) dplyr::left_join(x,y,by=“gene_id”),C)

 function(x,y)和dplyr::left_join(x,y,by="gene_id"整一个相当于参数中的f，C相当于参数中的x。

 本文仅为本人学习记录，可能有错误，代码也可能有许多值得优化的地方，希望各位路过的大佬批评指正。

 做完之后还是想感叹有工具就用是好文明啊！