Appendix VI. Genome Annotations

在前面的学习中，我们已经看到，几乎所有的生物信息分析都依赖于前人对于基因组的注释。
在本节中，我们提供了一些整理好的人类基因组注释的资源
我们还会系统的对处理基因组注释(如获得特定类型的RNA的gene model和转录本序列)的方法进行介绍。

参考注释资源

下载链接: https://cloud.tsinghua.edu.cn/d/ad22768345664924b202/?p=%2FFiles%2FAppendix%2FVII.%20Genome%20and%20Annotation&mode=list

region

source

# of genes

# of transcripts

annotation file

sequence file

rRNA $^1$

refseq Release 109

refseq.rRNA.gff3.gz

refseq.rRNA.fa.gz

mRNA

gencode v38

19,955

162,022

gencode.v38.mRNA.gff3.gz

gencode.v38.mRNA.fa.gz

tRNA

gencode v38

649

gencode.v38.tRNA.gff3.gz

gencode.v38.tRNA.fa.gz

snRNA

gencode v38

1,901

gencode.v38.snRNA.gff3.gz

gencode.v38.snRNA.fa.gz

snoRNA

gencode v38

943

gencode.v38.snoRNA.gff3.gz

gencode.v38.snoRNA.fa.gz

misc-srpRNA

gencode v38

680

gencode.v38.srpRNA.gff3.gz

gencode.v38.srpRNA.fa.,gz

misc-Y RNA

gencode v38

756

gencode.v38.Y_RNA.gff3.gz

gencode.v38.Y_RNA.fa.gz

misc-other RNAs $^2$

gencode v38

777

gencode.v38.annotation.others.gff3.gz

gencode.v38.others.fa.gz

pseudogene $^3$

gencode v38

14,537

18,016

gencode.v38.pseudogene.gff3.gz

gencode.v38.pseudogene.fa.gz

lncRNA (set 1)

gencode v38

16,888

47,675

gencode.v38.lncRNA.gff3.gz

gencode.v38.lncRNA.fa.gz

lncRNA (set 2)

mitranscriptome

60,763

175,271

mitranscriptome.v2.hg38.lncRNA.gtf.gz

mitranscriptome.v2.hg38.lncRNA.fa.gz

TUCP $^4$

mitranscriptome

3,630

10,978

mitranscriptome.v2.hg38.tucp.gtf.gz

mitranscriptome.v2.hg38.tucp.fa.gz

miRNA $^5$

miRbase 22.1

1,918

2,883

mirbase.22.1.hsa.gff3.gz

mirbase.22.1.pre-miRNA.fa.gz

piRNA

piRbase v3

77,242

hsa.pirbase.gold.v3.0.fa.gz

circRNA set 1

circbase

140,790

cricbase.hg19.fa.gz

circRNA set 2 $^6$

circbase & mioncocirc

34,376

cricbase.hg19.supported.by.oncocirc.fa.gz

promoter $^7$

EncodeBroadHmm

promoter.merged.bed.gz

enhancer $^8$

EncodeBroadHmm

enhancer.merged.bed.gz

intron $^9$

gencode v38 & mitranscriptome

introns.bed.gz

Transposable Element $^9$

TEtranscripts

GRCh38_GENCODE_rmsk_TE.gtf.gz

注:

这里列出的基因组注释均位于hg38坐标上

rRNA数量这里指的是refseq中注释的转录本数量，对应着5S，5.8S，18S，28S 4种rRNA和45S rRNA前体。refseq对rRNA在基因组上不同位置的多个拷贝都给出了单独的注释。
"other RNAs"这里指的是gencode中除了Y RNA和srpRNA之外的，被注释为misc_RNA的RNA，如RN7SK，Vault RNA等。
pseudogene指的是protein coding gene的假基因(rRNA等非编码RNA也会有假基因)。
mitranscriptome基于大量 TCGA 中癌症病人的 RNA-seq 数据对肿瘤转录组进行了组装和注释。TUCP (Transcript with Uncertain Coding Potential) 是mitranscriptome自己定义的一些不确定是否能编码蛋白的RNA。
miRNA基因数量指的是有注释的miRNA前体的数量，转录本数量对应成熟miRNA的数量。
circRNA set 2指的是在circbase和mioncocirc两个数据库中均被注释的环状RNA。
不同于转录本，promoter的边界通常很难严格的定义。ENCODE项目从ChIP数据出发获得了不同细胞系的组蛋白修饰等表观遗传信号(Jason Ernst, et al.,2011 )，据此用利用隐马尔可夫模型对基因组进行了segmentation，将每一个hidden state定义为一种染色质的状态，据此对每一个segment给出相应的注释。这里的"promoter"指的是在ENCODE HMM的至少一个细胞系对应的track中被注释成"Active_Promoter","Weak_Promoter"或"Poised_Promoter"的区域。在前面的章节(iii.3.chip-seq,iii.4.sequence-motif)中，我们曾经把TSS上下游的的一段区域定义为promoter，这也是一种常见的做法。
这里的"enhancer"指的是在ENCODE HMM的至少一个细胞系对应的track中被注释成"Strong_Enhancer"或"Weak_Enhancer"的区域。
这里的"intron"指的是gencode v38中的mRNA,pseudogene,lncRNA以及mitranscriptome中的lncRNA和TUCP去掉外显子之后剩下的部分。
我们这里列出的是用于分析转座元件表达量的工具TEtranscripts提供的注释。

获取人类参考基因组

目前来说，在生物信息学分析中最常用的人类基因组版本是GRCh38/hg38，一些相对早一些的研究用的是GRCh37/hg19。对于常规的生物信息分析，建议大家优先使用GRCh38/hg38。大家可以从gencode的ftp站点下载hg38的primary assembly:

http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_38/GRCh38.primary_assembly.genome.fa.gz

从其他很多地方也可以下载到hg38的序列。以下罗列了另外几个常用的链接:

ensembl: http://ftp.ensembl.org/pub/release-104/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
UCSC: https://hgdownload.cse.ucsc.edu/goldenpath/hg38/bigZips/hg38.fa.gz
NCBI - refseq: https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.26_GRCh38/GCF_000001405.26_GRCh38_genomic.fna.gz

hg38参考基因组中包含的不只有25条序列(22条常染色体，XY两条性染色体以及线粒体基因组),还包括了unplaced contigs等其他的序列。从这些来源下载到的参考基因组在主要的染色体上序列都是相同的，但是有时序列的名称不同。有人整理了这些序列id之间的对应关系，请大家参考https://github.com/dpryan79/ChromosomeMappings。

细心的同学可能会注意到ensembl的ftp站点除了"primary assembly"，还会提供一种"top level assembly"。primary assembly是top level assembly的一个子集，top level assembly除了主要的染色体和unplaced contigs等之外，还会包括一些haplotype和patch的序列。相应的名词解释请参考Assembly Terminology。对于我们前面介绍的常规数据分析，这些haplotype和patch序列基本是用不上的，所以通常会下载primary assembly。

获取常见基因/RNA类型的注释信息

同一个基因在不同基因组版本的上的坐标通常是不同的，所以大家一定要注意基因组注释的坐标相对于的是哪个版本的参考基因组。至于如何将别人提供的基于不同基因组版本的注释统一到我们自己用的版本上，我们将在后文进行介绍。

对于人类基因组,在常规的分析中大家常常用到的是gencode提供的注释。目前(2022.04.02)它已经更新到了第39版。我们同样从它的ftp站点下载gff3格式的注释文件:

https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_39/gencode.v39.annotation.gff3.gz

我们用zcat gencode.v39.annotation.gff3.gz | head看一看这个文件的开头几行:

##gff-version 3
#description: evidence-based annotation of the human genome (GRCh38), version 39 (Ensembl 105)
#provider: GENCODE
#contact: gencode-help@ebi.ac.uk
#format: gff3
#date: 2021-09-02
##sequence-region chr1 1 248956422
chr1	HAVANA	gene	11869	14409	.	+	.	ID=ENSG00000223972.5;gene_id=ENSG00000223972.5;gene_type=transcribed_unprocessed_pseudogene;gene_name=DDX11L1;level=2;hgnc_id=HGNC:37102;havana_gene=OTTHUMG00000000961.2
...

第二行明确写着这个注释是相对于human genome (GRCh38)坐标的。

在常规的二代测序分析中，用gencode提供的这一个注释文件已经足够了。更多的信息请参考https://www.gencodegenes.org/human/。

在RNA的生物信息分析中，我们除了protein coding的转录本，即mRNA，还会考虑一些non coding transcript，常见的有rRNA， tRNA，lncRNA，srpRNA，snRNA， snoRNA，Y RNA，miRNA和piRNA等。根据具体目的的不同，有时我们需要的是这些RNA在参考基因组上对应的genomic intervals，有时我们需要的是这些转录本自身的序列。例如，例如在我们的第一个大作业中，我们直接把小RNA测序的结果mapping到转录本序列上，所以需要的是序列信息。由于我们关注的绝大部分RNA都已经被组装到了hg38上，所以获得了前者的信息，就等价于获得了后者的信息。在下一节中，我们会介绍如何利用已有的工具从hg38中提取注释文件中给出的转录本。

gencode不止提供了mRNA的注释，还提供了lncRNA，多种非编码RNA，以及各种类型的pesudogene(假基因)的注释。对于上述的tRNA，lncRNA，srpRNA， snRNA， snoRNA和Y RNA我们仍然使用gencode的注释。如果我们只关注其中的一类RNA，可以根据gff文件第9列中的"gene_type"键值对进行筛选。

我们知道，脊椎动物有5S，5.8S，18S，28S四种rRNA。45S的前体先被RNA polymerase I转录出来，再加工成28S，5.8S和18S rRNA。5S rRNA由RNA polymerase III独立的转录产生。编码45S rRNA和5S rRNA的DNA序列在人类基因组中都有多个拷贝，所以这些rDNA有时会被注释成基因组中的一类重复序列(repeats)。对于rRNA，我们会发现，gencode.v39.annotation.gff3.gz提供了5Sr RNA和5.8S rRNA和两种线粒体rRNA(16S, 12S)的注释，并没有提供16S和28S rRNA的注释。refseq对 rRNA的注释更完整一些，所以这里我们下载refseq提供的注释，再筛选出rRNA对应的注释:

wget https://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/annotation/GRCh38_latest/refseq_identifiers/GRCh38_latest_genomic.gff.gz
gunzip GRCh38_latest_genomic.gff.gz
cat GRCh38_latest_genomic.gff | awk 'BEGIN{FS="\t";OFS="\t";}$0~"^#"{print;next;}$3=="rRNA"{print}' > refseq.rRNA.gff

在gencode.v39.annotation.gff3.gz中提供了线粒体tRNA的注释，其他染色体编码的tRNA在另外一个注释文件gencode.v39.tRNAs.gtf.gz中，需要我们单独下载:

wget http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_39/gencode.v39.tRNAs.gtf.gz

还有一些基因组注释的工作专门关注于特定类型的RNA，所以在专门研究这些类型的RNA时我们会优先参考这些来源的注释信息。

下载miRbase对人miRNA的注释:

# hsa means homo sapiens
wget -O hsa.mirbase.gff3 https://www.mirbase.org/ftp/CURRENT/genomes/hsa.gff3

下载piRbase提供的piRNA序列

wget -O hsa.pirbase.gold.v3.0.fa.gz http://bigdata.ibp.ac.cn/piRBase/download/v3.0/fasta/hsa.gold.fa.gz

有时我们会有一些更个性化的需求，如分析肿瘤中的long noncoding RNA，分析已知的环状RNA等等，这时我们也可以利用前人的注释工作。

mitranscriptome (2015, Matthew K Iyer et al.)从TCGA的肿瘤RNA-seq数据出发，注释了很多的转录本，并对转录本分成了protein coding，TUCP(transcript of unknown coding potential，这一类别是这项研究自己定义的，主要指一些编码蛋白的能力不确定的转录本)，lncRNA，read trough和pseudogene几类。这些转录本中有一部分是和gencode重合的，也有一部分在gencode中没有注释。这项研究主要关注于lncRNA，所以如果希望分析一些在肿瘤中表达，但在gencode中有没有注释的lncRNA，可以参考mitranscriptome的lncRNA注释:

wget https://mitranscriptome.org/download/mitranscriptome.gtf.tar.gz
tar xvzf mitranscriptome.gtf.tar.gz
cd mitranscriptome.gtf
zcat mitranscriptome.v2.gtf.gz > mitranscriptome.v2.hg19.gtf

虽然mitranscriptome文档里没有明确写明，但是这个gtf文件是在hg19坐标上的(!!)。如果我们想要对hg19上的基因注释和hg38上的基因注释进行比较，或者根据已经mapping到hg38上的bam文件count　mitranscriptome中有注释的基因表达，首先需要统一二者的坐标。我们这里将mitranscriptome.v2.gtf的坐标转换到hg38上，具体做法如下:

首先我们需要下载一个chain file，这个文件提供了不同基因组版本之间坐标的对应关系:

wget  http://hgdownload.soe.ucsc.edu/goldenPath/hg19/liftOver/hg19ToHg38.over.chain.gz

然后我们再用CrossMap这个工具进行坐标转换:

CrossMap.py gff hg19ToHg38.over.chain.gz mitranscriptome.v2.hg19.gtf mitranscriptome.v2.hg38.gtf

这时mitranscriptome.v2.hg38.gtf就已经在hg38的坐标上了，我们可以用cufflinks和stringtie自带的gffcompare(也可以独立的安装，请参考https://github.com/gpertea/gffcompare)来比较它和gencode v39的重合情况，或者用它和map到hg38上的bam文件count基因表达。

gffcompare -r gencode.v39.annotation.gff3 -o mitranscriptome.v2.hg38.query.gencode.v39 mitranscriptome.v2.hg38.gtf
# files with prefix mitranscriptome.v2.hg38.query.gencode.v39 contains information for overlaps of transcripts in these two annotations

circBase中注释了多个物种的一些已知的环状RNA，我们可以直接下载它提供的Putative spliced circRNA sequences:

wget http://www.circbase.org/download/human_hg19_circRNAs_putative_spliced_sequence.fa.gz

mioncocirc注释了人类肿瘤转录组中的环状RNA，也可供大家参考。

利用基因组注释获取转录本序列

假设我们有一组genomic intervals，每一个genomic interval都对应着一个转录本，我们可以很容易的自己实现一个python/perl/R脚本读取fasta文件，再提取指定的genomic interval。我们还可以把genomic interval存储为bed格式，再使用bedtools getfasta -s，更简洁的实现这个目的。

如果存在splicing的情况，我们同样可以实现一个脚本，parse gtf文件，把每一个exon对应的genomic interval按parent transcripts分组(每一个对应exon的record的第九列中会有一个键值对指明parent transcript id)，再把同一个转录本的exon按顺序拼接到一起。但是自己实现一个parse gtf文件的脚本还是相对比较繁琐的，也容易出错。一个更方便的做法是使用gffread这个工具。gffread和gffcompare一样是cufflinks和stringtie自带的工具，它同样也可以独立的安装，请参考https://github.com/gpertea/gffread。

samtools faidx GRCh38.primary_assembly.genome.fa #indexing fasta file to make the next step more efficient
gffread -W -M -F -G -A -O -E -g GRCh38.primary_assembly.genome.fa -w gencode.v39.tx.fa -d collapsed.info gencode.v39.annotation.gff3

以上命令的输入是参考基因组GRCh38.primary_assembly.genome.fa，基因组注释文件gencode.v39.annotation.gff3；输出文件是转录本序列gencode.v39.tx.fa，以及对完全重叠的转录本进行合并的信息collapsed.info。 gffread有大量的可选参数。除了从基因组中提取转录本序列，gffread另外一个很实用的功能是把gtf/gff文件转换成bed12格式(每行对应一个可以包含多个exon的转录本)。有兴趣的同学可以通过gffread -h详细了解。

PreviousAppendix V. Software and Tools

Last updated 9 months ago