转录组定量后计算FPKM
2021-02-08 本文已影响0人
谢京合
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)