不是所有人都像我这样喜欢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每条染色体基因,转录本的分布 !
可以针对上面的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提取转录本信息,外显子信息,这样就可以统计更多的基因特征啦。
Ensemble( ensembl.org网站是常用真核生物参考基因组来源之一 )能够对人类基因自动进行注释,包括人类,小鼠,斑马鱼,猪和大鼠等,也包括来自HAVANA的人工注释信息。 Ensembl是一项生物信息学研究计划,旨在开发种能够对真核生物基因组进行自动注释(automatic annotation)并加以维护的软件系统。该计划由英国Sanger研究所Wellcome基金会及欧洲分子生物学实验室所属分部欧洲生物信息学研究所共同协作运营。
Ensembl与NCBI的NCBI Map Viewer和UCSC是最为常用基因组检索数据库。
目前从事基因注释的机构组织有很多,这里列出的只是较为常用的几个。
Ensembl的通用基因注释有两种,一是Ensembl GeneBuild,它是自动化注释,速度快,实时更新,在不同物种上均适用;另一种是Wellcome基金会的 Havana (VEGA)小组的注释,它是手工注释,速度慢,但是准确,它依据的都是已经验证过的mRNA和蛋白序列来注释,比较费时。因此Ensembl基因组数据库 中,会有两种注释。
详细信息:http://vega.sanger.ac.uk/info/about/gene_and_transcript_types.html
人类和小鼠基因组的GTF文件与GENCODE计划发布的gene set文件相同。 The GENCODE project 的目标为对人类和小鼠基因组提供高质量的注释信息和实验确证。 The GENCODE gene sets被其他项目作为参考而广泛使用(如 1000 Genomes). 详细内容:https://www.gencodegenes.org/about.html
还可以结合 Gviz 进行可视化,下回分解