如果你下载了TCGA数据库里的count数据,想把它转换成FPKM应该怎么做呢?
参考文章:
1.Htseq Count To Fpkm
2.【原创】R语言实战:read counts如何转化为TPM和FPKM, TPM和FPKM相互转化
我主要是参考上面两篇文章的代码,因为这两篇文章的代码都不完整,但是合起来正好是一个完整的转换过程。
首先看一下什么是FPKM吧:
FPKM (Fragments Per Kilobase Million):Fragment Per Kilobase of transcript per Million mapped reads(每1百万个map上的reads中map到外显子的每1K个碱基上的Fragments个数)。
公式:
FPKM= read counts / (mapped reads (Millions) * exon length) #这里的exon length的单位是kb
这里要先得到的是每个exon的长度,为了求长度,我们先要有基因组注释文件,你可以在这里下载:https://gdc.cancer.gov/about-data/data-harmonization-and-generation/gdc-reference-files。这是GDC官网里的下载网址,进去以后,选择这个:
下载后解压,这个过程就不赘述了。
计算exon长度,并且保存成dataframe
> library(GenomicFeatures)
> txdb <- makeTxDbFromGFF("h38_GENCODE_v22_GTF.gtf",format="gtf")
#通过exonsBy获取每个gene上的所有外显子的起始位点和终止位点,然后用reduce去除掉重叠冗余的部分
#最后计算长度
> exons_gene <- exonsBy(txdb, by = "gene")
> exons_gene_lens <- lapply(exons_gene,function(x){sum(width(reduce(x)))})
得到这个长度后,打开看一下:
看一下数据类型:
> class(exons_gene_lens)
[1] "list"
> length(exons_gene_lens)
[1] 60483 #有60483个基因
是个列表,我们需要把它转成dataframe格式:
> exons_gene_lens1 <- as.data.frame(exons_gene_lens)
> dim(exons_gene_lens1)
[1] 1 60483
结果转完后一看,长这样:
再给它转置一下:
> exons_gene_lens1 <- t(exons_gene_lens1)
> dim(exons_gene_lens1)
[1] 60483 1
看一下转置之后的,顺眼多了:
把count矩阵和exon长度表合并
#load my counts
> counts = read.csv("TCGA_HNSCC_original_counts.csv",header = TRUE, sep = ",")
head(counts)
看看这个count矩阵和上面我们得到的exon长度的dataframe的基因名有什么区别?这个矩阵的基因名我们是处理过的,小数点后面的都删掉了,所以为了合并,exon的基因名也处理一下,另外exon长度表里的第一个基因,在我的count表里没有,总基因数差一个,所以还要把我exon长度表里的第一行删掉:
#删除第一行
> exons_gene_lens2 <- exons_gene_lens1[-1,]
> View(exons_gene_lens2)
> exons_gene_lens2 <- as.data.frame(exons_gene_lens2)
> rownames(exons_gene_lens2)
#把基因名的小数点后面的都去掉
> xc <- gsub("\\.(\\.?\\d*)","",rownames(exons_gene_lens2))
> xc
> rownames(exons_gene_lens2) = xc
#把exon长度表的列名设置为“Length”
> colnames(exons_gene_lens2) = "Length"
> View(exons_gene_lens2)
#合并count矩阵和基因长度的dataframe
> count_with_length <- cbind(counts,exons_gene_lens2)
> View(count_with_length)
> count_with_length <- count_with_length[,-1]
到这里,我们就把count矩阵和exon表合并了,exon length是新count表里最后一列。
计算FPKM
#先把length除以1000,就是上面公式里说的单位要kb
> kb <- count_with_length$Length/1000
> kb
#把新count矩阵里的前546列的数值都除以kb
> countdata <- count_with_length[,1:546]
> rpk <- countdata/kb
> rpk
#FPKM计算
> fpkm <- t(t(rpk)/colSums(countdata) * 10^6)
> head(fpkm)
#保存FPKM矩阵
> write.csv(fpkm,file="TCGA_HNSCC_count2FPKM.csv",sep="\t",quote=F)