变异信息那些事（下）

2019-01-03 本文已影响44人刘小泽

刘小泽写于19.1.3 又安心学了一天🤠学习时光过得真快
之前写到了VCF介绍（https://www.jianshu.com/p/3045aadbd723）和基本操作（https://www.jianshu.com/p/9aa8ca6e565c）
假设找到变异，然后呢？

变异位点的注释

我们得到变异位点，但仅仅是知道了它们在基因组上的位置信息和相关的碱基信息。那么还存在许许多多的疑问没有解决：

这个位点是在基因上吗？是内含子还是外显子区域？这个突变对基因功能产生了什么影响？对于转录翻译有没有影响？除了研究的样本，还有没有其他样本也出现了这个变异？有的话是什么人种，又是什么病例？

这些问题都要靠变异注释来解决

一般来说，变异注释分为：突变频率注释、变异的蛋白功能危害注释、剪切位点突变危害注释、突变相关的疾病注释

突变频率注释

做这个内容的数据库有许多，其中比较重要的有dbSNP、1000人基因组项目（1000 Genome）、ExAC、gnomeAD

dbSNP（The single-nucleotide polymorphism database）：http://www.ncbi.nlm.nih.gov/SNP/ NCBI与人类基因组研究所合作建立，包含了SNP、短重复序列、微卫星标记等来源、检测方法、基因型信息、上下游序列、人群分布频率等
1000G （千人基因组项目） 研究时限：2008-2015年；汇集30个人种、3904个样本WGS和WES测序结果。目前已被ANNOVAR收纳为变异位点在正常人群中进行突变频率注释的数据库，实际分析中也应该将1000G的不同人群作为control组进行疾病关联分析

它的构建总共分为4个阶段：
一、Pilot phase 【A map of human genome variation from population-scale sequencing Nature 467, 1061–1073 (28 October 2010)】
image.png
二、Phase one 【An integrated map of genetic variation from 1,092 human genomes Nature 491, 56–65 (01 November 2012)】

三、Phase two

四、Phase three 【A global reference for human genetic variation Nature 526, 68–74 (01 October 2015)； An integrated map of structural variation in 2,504 human genomes Nature 526, 75–81 (01 October 2015)】

参考：http://www.internationalgenome.org/about

含有两个下载镜像：ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/

ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/

其中Phase3所有的样本信息下载：ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/working/20130606_sample_info/
ExAC（Exome Aggregation Consortium）：整合了60706个人的WES测序数据及相关遗传信息，包含超过1000万种基因变异信息 http://exac.broadinstitute.org/ 。包括了AFR（African）、AMR（Admixed American）、EAS（East Asian）、FIN（Finnish）、NFE（NON-finnish European）、SAS（South Asia）等种群的突变频率（AF）信息http://exac.broadinstitute.org/faq
image.png
gnomeAD：（Genome Aggregation Database）博得研究所支持建立，包含了千人基因组、ESP数据库以及绝大部分的ExAC数据库。目前有125,748个外显子数据和15,708个基因组数据 http://gnomad.broadinstitute.org/，这些数据来自大型人群测序和疾病研究项目

变异的蛋白功能危害注释

PROVEAN：（Protein Variation Effect Analyzer）http://provean.jcvi.org/index.php 用来预测SNP或者InDel是否影响蛋白质的生物功能，不仅可以对CDS区域的非同义突变进行预测，还可以对CDS区域的非移码InDel对蛋白功能的影响进行预测，并将结果大致分为：危害、可以容忍、无害
SIFT：（Sorting Intolerant From Tolerant）https://sift.bii.a-star.edu.sg/ 根据氨基酸在蛋白序列中的保守程度来预测氨基酸的变化对蛋白功能造成的影响。其中保守程度是比对进化关系较近的蛋白序列得到，分值（SIFT-score）表示突变对蛋白序列的影响，分值越小越严重 ，一般认为：SIFT值小于0.05为有害（D：Deleterious），大于0.05表示容忍（T：Tolerance）
Polyphen2_HAVR: (Polymorphism Phenotyping v2) http://genetics.bwh.harvard.edu/pph2/dokuwiki/downloads 根据HumanVar数据库预测突变对蛋白的影响，来诊断孟德尔遗传病。分值表示SNP导致蛋白结构或功能改变的可能性，越大越严重
Polyphen2_HDIV：根据HumanDiv数据库预测，分值越大越严重
LRT：也是基于序列保守性进行预测（像SIFT和Polyphen）http://www.genetics.wustl.edu/jflab/lrt_query.html 。对每一个测试的密码子，LRT将来自31个物种的氨基酸进行比对来预测突变的危害。结果的有害突变（D：Deleterious）表示：突变来自高度保守的密码子；突变氨基酸在其他比对的真核哺乳动物中不存在。中性突变（N: Neutral）表示：突变发生在非高度保守的密码子；突变的氨基酸至少在一个进行比对的真核哺乳动物中发现

剪切位点突变危害注释

如果突变发生在剪切位点附近，我们可以判断它对剪切的危害。可以用的软件有：DbscSNV、Spidex、MaxEntScan

DbscSNV：属于VEP（Variant Effect Predictor）插件，http://asia.ensembl.org/info/docs/tools/vep/script/vep_plugins.html#plugins_existing

由AdaBoost与Random Forest开发，它根据突变前后分值的变化来预测剪切位点的突变危害性
Spidex：http://www.openbioinformatics.org/annovar/spidex_download_form.php 基于深度学习，因此预测的剪切变异可能距离常规的剪切位点比较远（这一点和DbscSNV不同）
MaxEntScan：对5'剪切位点附近的6bp内含子与3bp的编码区(http://genes.mit.edu/burgelab/maxent/Xmaxentscan_scoreseq.html)以及3‘ 剪切位点附近20bp的内含子与3bp的编码区内突变进行预测(http://genes.mit.edu/burgelab/maxent/Xmaxentscan_scoreseq_acc.html)，按照突变前后分值变化来得到结论（认为有危害：突变后比突变前分值降低15%以上）

突变相关的疾病注释

OMIM：（Online Mendelian Inheritance in Man）https://www.omim.org/在线人类孟德尔遗传信息数据库，包含了遗传性的基因疾病信息与表型信息，目前收录了16000多个基因词条和5400多表型词条
HGMD：（The Human Gene Mutation Database）1996年创立的人类基因突变数据库，目前包括240,269个变异，覆盖9976个基因。收集的突变包含了SNP、InDel、CNV、SV、基因重组等，可以说是遗传病变异检测金标准数据库。有两个版本，一个是免费的学术public版，但更新慢（http://www.hgmd.cf.ac.uk/ac/index.php）；另一个是收费可试用的Professional版（https://www.qiagenbioinformatics.com/products/human-gene-mutation-database/），包含的变异数量也更多
ClinVar：2013年创立，是一个已报道突变与疾病表型关联数据库，https://www.ncbi.nlm.nih.gov/clinvar/。数据主要来源是OMIM、dbSNP、locus specific database等开源数据库，对变异位点的审核比较缺乏，因此会包含报道中冲突的致病位点
HGMD and ClinVar: Avoiding the Knowledge Blind Spot

练一个工具--snpEff

conda安装

$ conda install -y snpeff

看看snpeff目前有什么数据库

# 目前有42791个数据库
$ snpEff databases > listing.txt

找到Ebola相关数据库

$ cat listing.txt | grep Homo_sapiens
#GRCh37.75                                                   Homo_sapiens                                                                                               http://downloads.sourceforge.net/project/snpeff/databases/v4_3/snpEff_v4_3_GRCh37.75.zip

下载数据库

$ snpEff download GRCh37.75  
# 或者
$ wget -c http://downloads.sourceforge.net/project/snpeff/databases/v4_3/snpEff_v4_3_GRCh37.75.zip

进行注释

$ snpEff GRCh37.75 subset_hg19.vcf > subset_hg19.anno.vcf

结果

主要还是看官方manual，得到的新注释的vcf中最明显的变化就是INFO列增加了一个字段ANN，默认ANN中又会给出几种信息

Allele：列出突变碱基
Annotation：列出Sequence Ontology中的条目，表示变异的后果或者影响（effect or consequence），例如intron_variant；如果有多个，用&连接intron_variant&nc_transcript_variant
Putative impact：变异位点的危害程度大小，四个取值：HIGH、MODERATE、LOW、MODIFIER
image.png
Gene name：HGNC官方基因名
Gene ID
Feature type：feature信息（如transcript, motif, miRNA等等），如果是组织特异性信息，可以添加细胞类型或者组织信息等（用冒号隔开），如H3K4me3:HeLa-S3
Feature ID ：根据type来决定ID，比如type是transcript，那么就是Transcript ID，另外还有Motif ID、miRNA、ChipSeq peak、Histone mark等
Transcript biotype：Ensembl数据库的转录本类型（Coding / Noncoding）
Rank：变异位点出现在基因区域时，会给出位点在exon/intron的第几位。例如，变异位点出现在某基因的第2个exon上，而这个基因共有10个exon，因此就显示2/10
HGVS.c：根据HGVS（http://www.hgvs.org/）标准命名的基因水平变异
HGVS.p：根据HGVS标准命名的蛋白水平变异（前提是变异位点在编码区）。如果Transcript ID在feature ID中表示出来了，这里就可以省略
cDNA_position（可选cDNA_len）：变异位点在cDNA的位置（或cDNA的总长度）
CDS_position （可选CDS_len）：变异位点在CDS的位置（或CDS长度）
Protein_positoin (可选Protein_len)：位点在AA的位置（或AA总长度）
Distance：变异位点与最接近的feature的距离，例如位点在exon，会给出与最近的内含子的距离；位于基因间区会给出与最近基因的距离

image.png
Errors，Warnings：注释可靠性评估【见官网】