从果子老师的从TCGA数据中提取lncRNA并进行下游分析一文中，我学到了如何从TCGA表达谱数据中提取RNA。老师的教程已经写的很详细了，我再补充一点：如何下载GTF注视文件呢？打开终端
wegt -m ftp://ftp.ensembl.org/pub/release-90/gtf/homo_sapiens/
然而，一个bug卡了我一天，这都怪我没学好分子生物学。本科是学林学的，整天在林子里拈花惹草，硕士莫名进了一个癌生物学湿实验室，又莫名的做了生信。研一把课全翘了留在实验室写代码，偶尔去上课就捧着Mac air躲在班级最后一排鬼鬼祟祟，我室友后来说发白光的小苹果像是在嘲笑老师。
这个bug是，在结束果子老师的教程全文后，我将lncRNA和mRNA的表达谱提出来，把他们的gene symbol 做交集，竟发现有一万个重复的基因。
于是，我回去查bug。
GTF注释文件如下：
提取mRNA:
mRNA_exprSet<- gtf_df %>% dplyr::filter(type=="gene",gene_biotype=="protein_coding")
type是基因，gene_biotype是蛋白质编码基因，这没毛病
提取lncRNA:
ncRNA <-c("processed_transcript","sense_overlapping","lincRNA","3prime_overlapping_ncRNA","sense_intronic","bidirectional_promoter_lncRNA","non_coding")
LncRNA_exprSet<- gtf_df %>% dplyr::filter(type=="transcript",transcript_biotype %in% ncRNA)
这时将type设为transcript，果子老师在原文中解释说，有一部分编码基因，转录成非编码RNA了，“一个编码基因也有可能转录出非编码基因的，对么?”
这句话我拿去问师兄们，他们一致说，怎么可能。
编码基因(protein_coding gene)在转录时，外显子重新编码到一起，就算有可变剪切体，也是编码出多种蛋白质，怎么可能转录出非编码RNA呢？
为了验证蛋白质编码基因能否转录出非编码RNA，我将代码改为
ncRNA <-c("processed_transcript","sense_overlapping","lincRNA","3prime_overlapping_ncRNA","sense_intronic","bidirectional_promoter_lncRNA","non_coding")
LncRNA_exprSet<- gtf_df %>% dplyr::filter(type=="transcript", gene_biotype="protein_coding",transcript_biotype %in% ncRNA)
type=="transcript" 说明我要的是转录本
gene_biotype="protein_coding" 说明基因类型要的是蛋白质编码基因
transcript_biotype %in% ncRNA 是要指定的几种ncRNA类型
这逻辑没毛病
跑了之后，
标红处可以看出，蛋白质编码基因，对应着的转录本类型有且只有processed_transcript
那processed_transcript又是什么呢？
百度肯定不会告诉我们答案，呵，百度，你对力量一无所知。
正确上网后
原网址：http://vega.archive.ensembl.org/info/about/gene_and_transcript_types.html
可见，processed_transcript不含ORF，包含LncRNA, ncRNA和unclassified processed transcript。由于transcript_biotype中已经含有"sense_overlapping", "lincRNA", "3prime_overlapping_ncRNA"等很多种，所以我推测， 此处的processed_transcript就是未分类的转录本。
因此果子老师的那句话是对的，至少GTF告诉我们，一个编码基因也有可能转录出非编码基因的。
由于我的分子生物学底子非常弱，还请各位前辈多多指教。
那么，如果想从TCGA里提lncRNA，不要蛋白质编码基因转录出来的该怎么写
LncRNA_exprSet<- gtf_df %>%
dplyr::filter(type=="transcript",gene_biotype!="protein_coding",
transcript_biotype %in% ncRNA) %>%
dplyr::select(c(gene_name,gene_id,transcript_biotype)) %>%
dplyr::distinct()
mRNA_exprSet<- gtf_df %>% dplyr::filter(type=="gene",gene_biotype=="protein_coding") %>%
dplyr::select(c(gene_name,gene_id,gene_biotype))
网友评论