如何把从TCGA下载的HTseq count转换成FPKM

如果你下载了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)
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
禁止转载,如需转载请通过简信或评论联系作者。
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 203,271评论 5 476
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,275评论 2 380
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 150,151评论 0 336
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,550评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,553评论 5 365
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,559评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,924评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,580评论 0 257
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,826评论 1 297
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,578评论 2 320
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,661评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,363评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,940评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,926评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,156评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,872评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,391评论 2 342

推荐阅读更多精彩内容