转录组数据分析—Count转换为FPKM或TPM
2023-02-28 本文已影响0人
生信助手
- Count:高通量测序中比对到exon上的reads数。可以使用featureCounts、HTseq-count等软件进行计算
- FPKM:全称为Fragments Per Kilobase Million,Fragments Per Kilobase of exon model per Million mapped fragments(每千个碱基的转录每百万映射读取的fragments)。通俗讲,把比对到的某个基因的Fragment数目,除以基因的长度,其比值再除以所有基因的总长度。
- TPM:TPM的计算方法也同RPKM/FPKM类似,首先计算每个基因的表达值,去除基因长度的影响。随后计算每个基因的表达量的百分比,最后再乘以10^6,TPM可以看作是RPKM/FPKM值的百分比。相当于重新标准化的文库,保证每个样本中所有TPM的总和是相同的。
count值为什么要用FPKM或TPM标准化?
基因长度长或者测序数据量多的基因所得到的count值会越大,但这并不能反映该基因的表达高,而FPKM或TPM法能消除基因长度和测序量差异对计算基因表达的影响,计算得到的基因表达量可直接用于比较不同样品间的基因表达差异和不同基因间表达高低的比较。
获取基因长度
基因长度需要使用R包GenomicFeatures以及基因组注释文件
#######################获取基因长度文件############################
library(GenomicFeatures)
txdb <- makeTxDbFromGFF("F:\\reference\\homo\\Homo_sapiens.GRCh38.104.gtf",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))
n <- data.frame(names = row.names(n),n)
names(n)<-c('genes','length')
write.table(n,'gene_length.txt',col.names=T,row.names=F,quote=F,sep='\t')
Count转换TPM
######################转换为TPM####################################
library(dplyr)
count <- read.table("rawread.txt",header=T,sep='\t',check.names=F)
merge <- left_join(count,length,by="genes")
merge <- na.omit(merge)
rownames(merge)<-merge[,1]
merge<-merge[,-1]
kb <- merge$length / 1000
countdata <- merge[,1:9]
rpk <- countdata / kb
tpm <- t(t(rpk)/colSums(rpk) * 1000000)
tpm <- data.frame(genes = row.names(tpm),tpm)
ensembl = tpm[,1]
geneSymbol <- bitr(ensembl, fromType = "ENSEMBL", #fromType为输入数据ID类型
toType = "SYMBOL", #toType是指目标类型
OrgDb = org.Hs.eg.db)#Orgdb是指对应的注释包,这是人的注释包,相关注释包请更换
tpm <- merge(geneSymbol,tpm,by.x="ENSEMBL",by.y="genes")
write.table(tpm,file="tpm.xls",row.names=F,sep="\t",quote=F)
Count转换FPKM
######################转换为FPKM###################################
fpkm <- t(t(rpk)/colSums(countdata) * 10^6)
fpkm <- data.frame(genes = row.names(fpkm),fpkm)
ensembl = fpkm[,1]
geneSymbol <- bitr(ensembl, fromType = "ENSEMBL", #fromType为输入数据ID类型
toType = "SYMBOL", #toType是指目标类型
OrgDb = org.Hs.eg.db)#Orgdb是指对应的注释包,这是人的注释包,相关注释包请更换
fpkm <- merge(geneSymbol,fpkm,by.x="ENSEMBL",by.y="genes")
write.table(fpkm,file="fpkm.xls",row.names=F,sep="\t",quote=F)
转换的过程中都添加了基因注释,也可以选择不注释基因。
有了标准化的数据,后续有很多分析就可以顺利进行啦,比如时序性分析、GSEA等等