使用refGenome加上dplyr玩转gtf文件
使用refGenome加上dplyr玩转gtf文件
不是所有人都像我这样喜欢linux的黑白命令行，但是他们仍然是可以处理NGS数据的，比如最常用的gtf格式的基因组注释文件：
library(Rsubread)
# 推荐从ENSEMBL上面下载成套的参考基因组fa及基因注释GTF文件
dir='~/data/project/release1/Genomes/'
## 这个文件有1.4G哦！！！
gtf <- file.path(dir,'Homo_sapiens.GRCh38.82.gtf')
if(!require(refGenome)) install.packages("refGenome")
# create ensemblGenome object for storing Ensembl genomic annotation data
ens <- ensemblGenome()
# read GTF file into ensemblGenome object
read.gtf(ens,useBasedir = F,gtf)
## 耗时2分钟，取决于电脑配置
class(ens)
# counts all annotations on each seqname
tableSeqids(ens)
# returns all annotations on mitochondria
extractSeqids(ens, 'MT')
# summarise features in GTF file
tableFeatures(ens)
# create table of genes
my_gene <- getGenePositions(ens)
dim(my_gene)
# 这样就轻松拿到了所有基因的坐标信息啦
当然，这个gtf是有非常多的值得探索的地方，比如可以完成http://www.biotrainee.com/thread-626-1-1.html 我在生信技能树»生信技能树›互动作业›脚本能力实践›生信人必练的200个数据处理任务›生信编程直播第三题：hg38每条染色体基因,转录本的分布 ！
深度探索gtf来理解基因结构
可以针对上面的my_gene这个简单的数据框进行探索：
> colnames(my_gene)
[1] "id" "seqid"
[3] "source" "start"
[5] "end" "score"
[7] "strand" "frame"
[9] "ccds_id" "exon_id"
[11] "protein_id" "transcript_support_level"
[13] "protein_version" "havana_transcript_version"
[15] "transcript_biotype" "transcript_version"
[17] "transcript_source" "transcript_id"
[19] "havana_gene" "exon_number"
[21] "gene_name" "gene_biotype"
[23] "tag" "havana_gene_version"
[25] "exon_version" "gene_id"
[27] "gene_source" "havana_transcript"
[29] "transcript_name" "gene_version"
下面是一些简单的考题啦，之前是强调学生用shell脚本知识完成， 现在换成R，同样的效果。
- 不同类型的基因的数量，提示：
table(my_gene$gene_biotype) - 不同染色体的基因的数量，并且按照不同基因类型分组后继续统计
- 所有
protein_coding类型的基因的长度分布情况 - 所有的
非protein_coding类型的基因的长度分布情况 - 把全部基因的坐标画在染色体上面
还可以根据gtf提取转录本信息，外显子信息，这样就可以统计更多的基因特征啦。
GTF背景知识
Ensemble（ ensembl.org网站是常用真核生物参考基因组来源之一 ）能够对人类基因自动进行注释，包括人类，小鼠，斑马鱼，猪和大鼠等，也包括来自HAVANA的人工注释信息。
Ensembl是一项生物信息学研究计划，旨在开发一种能够对真核生物基因组进行自动注释(automatic annotation)并加以维护的软件系统。该计划由英国Sanger研究所Wellcome基金会及欧洲分子生物学实验室所属分部欧洲生物信息学研究所共同协作运营。
Ensembl与NCBI的NCBI Map Viewer和UCSC是最为常用基因组检索数据库。
Ensembl 与NCBI Map Viewer和UCSC最大区别表现在以下5点：
- a.Ensembl的基因数据集是依据mRNA和蛋内序列的数据信息白动注释的。数据来源为新的基因组数据，UniProt/SwissProt和UniProt/TrEMBL的蛋白序列，NCBI的RefSeq里的DNA和蛋白序列和EMBL的cDNA序列。
- b.Ensembl是一个开源（Perl API ）的全自动的基因注释软件系统，很多网站都采用Ensembl这套软件系
- c.Ensembl拥存其特有的BioMart功能。BioMart可以依据设定的要求对基 因组进行条件性检索，检索的结果吋以以图表的形式给出。
- d.与其它数据库相整合，比如DAS。
- e.基因组间的比较分析。
基因注释机构
目前从事基因注释的机构组织有很多，这里列出的只是较为常用的几个。
- Ensembl:目的是做出最好的基因注释集。
2.Havana (VEGA):是桑格中心的一个基因注释组织，它的目标和Ensembl—致，因此，结合得也最紧密。 - HGNC -给出人类基因唯一的名字和符号。
- UniProt 主要集中于蛋白质的信息注释。
Ensembl的通用基因注释有两种，一是Ensembl GeneBuild，它是自动化注释，速度快，实时更新，在不同物种上均适用；另一种是Wellcome基金会的 Havana (VEGA)小组的注释，它是手工注释，速度慢，但是准确，它依据的都是已经验证过的mRNA和蛋白序列来注释，比较费时。因此Ensembl基因组数据库 中，会有两种注释。
Havana (VEGA)小组的注释常有以下几种类型：
详细信息：http://vega.sanger.ac.uk/info/about/gene_and_transcript_types.html
- Protein coding： 包括开放阅读框 (ORF).
- Processed transcript：没有开放阅读框（ORF）
- Pseudogene：假基因，是指脱氧核糖核酸（DNA）的碱基序列中，一段与其他生物体内已知的基因序列非常相似的片段。但是这个片段由于移码突变或者无义突变破坏了ORF，无法发挥原有的基因功能，也就是无法制造出蛋白质
- IG gene：免疫球蛋白家族基因
- TR Gene：T细胞受体基因
- TEC (To be Experimentally Confirmed)
人类和小鼠基因组的GTF文件与GENCODE计划发布的gene set文件相同。
The GENCODE project 的目标为对人类和小鼠基因组提供高质量的注释信息和实验确证。
The GENCODE gene sets被其他项目作为参考而广泛使用(如 1000 Genomes).
详细内容：https://www.gencodegenes.org/about.html
还可以结合 Gviz 进行可视化，下回分解
