如果你下載了TCGA數(shù)據(jù)庫里的count數(shù)據(jù)馆类，想把它轉(zhuǎn)換成FPKM應該怎么做呢？
參考文章：
１．Htseq Count To Fpkm
２．【原創(chuàng)】R語言實戰(zhàn)：read counts如何轉(zhuǎn)化為TPM和FPKM, TPM和FPKM相互轉(zhuǎn)化

我主要是參考上面兩篇文章的代碼，因為這兩篇文章的代碼都不完整，但是合起來正好是一個完整的轉(zhuǎn)換過程。

首先看一下什么是FPKM吧：
FPKM (Fragments Per Kilobase Million)：Fragment Per Kilobase of transcript per Million mapped reads（每1百萬個map上的reads中map到外顯子的每1K個堿基上的Fragments個數(shù)）梯浪。

公式：

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官網(wǎng)里的下載網(wǎng)址眠砾，進去以后虏劲，選擇這個：

下載后解壓，這個過程就不贅述了。

計算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)))})

得到這個長度后，打開看一下：

看一下數(shù)據(jù)類型：

> class(exons_gene_lens)
[1] "list"
> length(exons_gene_lens)
[1] 60483　　#有60483個基因

是個列表堡掏，我們需要把它轉(zhuǎn)成dataframe格式：

> exons_gene_lens1 <- as.data.frame(exons_gene_lens)
> dim(exons_gene_lens1)
[1]     1 60483

結(jié)果轉(zhuǎn)完后一看应结，長這樣：

再給它轉(zhuǎn)置一下：

> exons_gene_lens1 <- t(exons_gene_lens1)
> dim(exons_gene_lens1)
[1] 60483     1

看一下轉(zhuǎn)置之后的，順眼多了：

把count矩陣和exon長度表合并

#load my counts
> counts = read.csv("TCGA_HNSCC_original_counts.csv",header = TRUE, sep = ",")
head(counts)

看看這個count矩陣和上面我們得到的exon長度的dataframe的基因名有什么區(qū)別泉唁？這個矩陣的基因名我們是處理過的摊趾，小數(shù)點后面的都刪掉了，所以為了合并游两，exon的基因名也處理一下，另外exon長度表里的第一個基因漩绵，在我的count表里沒有贱案，總基因數(shù)差一個，所以還要把我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)
#把基因名的小數(shù)點后面的都去掉
> 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列的數(shù)值都除以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)