如果你下载了TCGA数据库里的count数据，想把它转换成FPKM应该怎么做呢？
参考文章：
１．Htseq Count To Fpkm
２．【原创】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)
网友评论