转录组定量后计算FPKM

1、安装GenomicFeatures

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("GenomicFeatures")

但是
报错了。

ERROR: failed to lock directory C:\Users\xiejinghe\Documents\R\win-library\4.0 
removeing 00LOCK

解决办法
进入C:\Users\xiejinghe\Documents\R\win-library\4.0
删掉00LOCK文件夹
重新安装。

2、获取每个基因的长度

library(GenomicFeatures)
setwd("E:/4/")
#注释文件获取
txdb <- makeTxDbFromGFF("gencode.v22.annotation.gtf", format="gtf")
exons_gene <- exonsBy(txdb,by="gene")
exons_gene_lens <- lapply(exons_gene, function(x){sum(width(reduce(x)))})
#查看下
class(exons_gene_lens)
length(exons_gene_lens)
#转换成data frame
exons_gene_lens1 <- as.data.frame(exons_gene_lens)
class(exons_gene_lens1)
dim(exons_gene_lens1)
#转置
exons_gene_lens1 <- t(exons_gene_lens1)
dim(exons_gene_lens1)

##此时基因名字后面还存在小数点，所以需要转换一下。
#然后替换名字
#也就可以不替换，直接跳过之一步骤。根据你自己的需求。
xc <- gsub("\\.\\d*","",rownames(fpkm))
rownames(exons_gene_lens1) <- xc

##将文件输出
write.csv(exons_gene_lens1, file="gene_Length.csv")

3、基因名字的转换（根据需求选择）

#基因名字转换，将gene_Length第一列ensemble转换成symble
#根据自己需求
#perl 03.GeneName2Symble.pl

4、计算FPKM

##经过名字转化之后，重新载入基因长度的数据。
exons_gene_lens2 <- read.csv("gene_Length.csv", header=T)
#exons_gene_lens2 <- read.table("gene_Length.txt", header=T)
colnames(exons_gene_lens2) <- c("id","Length")
View(exons_gene_lens2)

##载入reads counts的原始数据。
readscounts <- read.table("mRNAmatrix.txt", header=T)

##如果你自己的矩阵名字更刚才的gene_Length.csv一样，就不需要。
##如果不一样，要么转化你自己的，要么转换gene_Length.csv的。
xc <- gsub("\\.(\\.?\\d)","",readscounts$id)
readscounts$id <- xc
readscounts[1:10,1:2]

###将基因长度和原始矩阵合并
##取出两个样本之间都有的，也就是交集。
mycounts <- merge(exons_gene_lens2, readscounts, by="id", all=F)
dim(mycounts)
mycounts[1:3,1:3]

##将第一列变成列名
##然后删除第一列。
rownames(mycounts)<-mycounts[,1]
mycounts<-mycounts[,-1]
mycounts[1:3,1:3]

##TPM计算
kb <- mycounts$Length / 1000
write.csv(kb, file="kb.csv")
##筛选出样本的表达数据
# 这一步骤是筛选出所有样本的表达量数据，第一列是length，所以不要。
countdata <- mycounts[,2:1000]
##计算每个值的rpk
rpk <- countdata / kb

##计算TPM
tpm <- t(t(rpk)/colSums(rpk) * 1000000)
head(tpm)

##计算FPKM
fpkm <- t(t(rpk)/colSums(countdata) * 10^6) 
head(fpkm)
write.csv(fpkm, file="mRNA_fpkm.csv")

##将FPKM转化成TPM
fpkm_to_tpm = t(t(fpkm)/colSums(fpkm))*10^6
head(fpkm_to_tpm)