05

用GEMINI来探索vcf格式的突变数据

第一次听说这个软件,是一个香港朋友推荐的:http://davetang.org/muse/2016/01/13/getting-started-with-gemini/ 他写的很棒,但是我当初以为是一个类似于SQLite的数据库浏览模式,所以没在意。实际上,我现在仍然觉得这个软件没什么用!

软件官网有详细的介绍:https://gemini.readthedocs.io/en/latest/

而且提供丰富的教程:

We recommend that you follow these tutorials in order, as they introduce concepts that build upon one another.

  • Introduction to GEMINI, basic variant querying and data exploration. html pdf
  • Identifying de novo mutations underlying Mendelian disease html pdf
  • Identifying autosomal recessive variants underlying Mendelian disease html pdf
  • Identifying autosomal dominant variants underlying Mendelian disease html pdf
  • Other GEMINI tools html pdf

软件本身并不提供注释,虽然它的功能的确包括注释,号称可以利用(ENCODE tracks, UCSC tracks, OMIM, dbSNP, KEGG, and HPRD.)对你的突变位点注释,比如你输入1       861389  .       C       T       ,它告诉你这个突变发生在哪个基因,对蛋白改变如何?是否会产生某些疾病?

虽然它本身没有注释功能,但是它会调用snpEFF或者VEP进行注释,你需要自己先学习它们。

1

软件安装:

GEMINI是用python写的,有一个小脚本可以自动完成安装过程:

7.3K May  4 14:44 gemini_install.py

下载这个脚本,然后安装即可

wget https://github.com/arq5x/gemini/raw/master/gemini/scripts/gemini_install.py

python gemini_install.py $tools $data

PATH=$tools/bin:$data/anaconda/bin:$PATH

where $tools and $data are paths writable on your system.

我把$tools用的就是当前文件夹,$data也是当前文件夹下面的gemini文件夹。

这样就会在当前文件夹下面生成两个文件夹,bin是存储程序,gemini是存储数据用的,而且注意要把bin目录的全路径添加到环境变量!

输入数据:

我们可以直接下载软件作者提供的测试数据

首先是22号染色体的所有突变位点经过WEP注释的文件

然后是一个三口直接的突变ped格式数据

数据存放在亚马逊云,所有的教程pdf也在

http://s3.amazonaws.com/gemini-tutorials

如果是你自己的vcf文件,需要自己用VEP注释一下

1

运行命令:

2

结果解读:

产生是chr22.db就是一个数据库格式的文件,但是需要用gemini 来进行查询,个人认为,并没有多大意思!

你只要熟悉mySQL等SQL语言,完全可以自己来!

05

用VEP对vcf格式的突变数据进行注释

VEP是国际三大数据库之一的ENSEMBL提供的,也是非常主流和方便,但它是基于perl语言的,所以在模块方面可能会有点烦人。跟snpEFF一样,也是对遗传变异信息提供更具体的注释,而不仅仅是基于位点区域和基因。如果你熟悉外显子联盟这个数据库EXAC(ExAC.r0.3.sites.vep.vcf.gz),你可以下载它所有的突变记录数据,看看它对每个变异位点到底注释了些什么,它就是典型的用VEP来注释的。

随便一个位点,注释了如此多的信息!~~~

1       861389  .       C       T       5621.53 PASS    AC=4;AC_AFR=0;AC_AMR=0;AC_Adj=4;AC_EAS=0;AC_FIN=0;AC_Het=4;AC_Hom=0;AC_NFE=3;AC_OTH=1;AC_SAS=0;AF=3.300e-05;AN=121216;AN_AFR=10212;AN_AMR=11516;AN_Adj=119730;AN_EAS=8606;AN_FIN=6594;AN_NFE=65414;AN_OTH=890;AN_SAS=16498;BaseQRankSum=2.78;ClippingRankSum=-2.380e-01;DP=1488042;FS=7.913;GQ_MEAN=62.49;GQ_STDDEV=14.73;Het_AFR=0;Het_AMR=0;Het_EAS=0;Het_FIN=0;Het_NFE=3;Het_OTH=1;Het_SAS=0;Hom_AFR=0;Hom_AMR=0;Hom_EAS=0;Hom_FIN=0;Hom_NFE=0;Hom_OTH=0;Hom_SAS=0;InbreedingCoeff=-0.0004;MQ=59.70;MQ0=0;MQRankSum=0.198;NCC=409;QD=15.11;ReadPosRankSum=0.561;VQSLOD=0.392;culprit=FS;DP_HIST=373|361|219|102|34981|16744|5493|1367|498|210|121|54|32|18|13|9|3|3|3|4,0|0|0|0|0|0|0|0|0|0|0|1|0|0|0|0|0|1|1|1;GQ_HIST=26|352|26|24|472|62|71|34|23|29|34|16|44468|8058|2176|2147|1116|370|365|739,0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|4;CSQ=T|ENSG00000187634|ENST00000420190|Transcript|missense_variant|157|68|23|P/L|cCg/cTg||1||1|SAMD11|HGNC|28706|protein_coding|||ENSP00000411579||Q5SV95_HUMAN&I7FV93_HUMAN&A6PWC8_HUMAN|UPI000155D47C|deleterious(0)|probably_damaging(0.999)|2/7|||ENST00000420190.1:c.68C>T|ENSP00000411579.1:p.Pro23Leu||||||||||||||||||,T|ENSG00000268179|ENST00000598827|Transcript|missense_variant|211|211|71|G/R|Ggg/Agg||1||-1|AL645608.1|Clone_based_ensembl_gene||protein_coding|YES||ENSP00000471152||M0R0C9_HUMAN|UPI0000D61E05||probably_damaging(0.98)|6/6|||ENST00000598827.1:c.211G>A|ENSP00000471152.1:p.Gly71Arg||||||||||||||||||,T|ENSG00000187634|ENST00000437963|Transcript|missense_variant|128|68|23|P/L|cCg/cTg||1||1|SAMD11|HGNC|28706|protein_coding|||ENSP00000393181||Q5SV95_HUMAN&I7FV93_HUMAN|UPI000155D47B|deleterious(0)|probably_damaging(0.999)|2/5|||ENST00000437963.1:c.68C>T|ENSP00000393181.1:p.Pro23Leu||||||||||||||||||,T|ENSG00000187634|ENST00000342066|Transcript|missense_variant|151|68|23|P/L|cCg/cTg||1||1|SAMD11|HGNC|28706|protein_coding|YES|CCDS2.2|ENSP00000342313|SAM11_HUMAN|Q5SV95_HUMAN&I7FV93_HUMAN&A6PWC8_HUMAN|UPI0000D61E04|deleterious(0)|probably_damaging(0.999)|2/14|||ENST00000342066.3:c.68C>T|ENSP00000342313.3:p.Pro23Leu||||||||||||||||||,T||ENSR00000528850|RegulatoryFeature|regulatory_region_variant|||||||1||||||regulatory_region|||||||||||||||||||||||||||||||

头文件里面有对每一列的详细介绍,包括突变的标准格式

HGVS.c   --》ENST00000420190.1:c.68C>T

HGVS.p –》ENSP00000411579.1:p.Pro23Leu

还有该突变对蛋白功能的影响,包括sift和polyphen的打分~~~

不多说了,直接介绍该软件如何使用吧!

 

软件安装:

最新版是84:http://useast.ensembl.org/info/docs/tools/vep/script/vep_download.html

然后进入目录用perl的形式来安装这个软件:perl INSTALL.pl 即可

安装时其实有很多参数可以选择的,请仔细阅读介绍;http://useast.ensembl.org/info/docs/tools/vep/script/vep_download.html

 

前提是你已经安装好了两个模块!

perl -e 'use DBD::mysql'

perl -e 'use Archive::Extract'

如果不报错,就证明你已经安装过这些模块,如果报错,去搜索我以前关于perl模块的博客吧,不是很简单的事情。

By default the script will install the cache files in the ".vep" subdirectory of the user's home area. Using this option users can configure where cache files are installed.

我不想把cache文件放在默认的$HOME/.vep/下面,所以我安装的时候稍微做了更改

 

下载完了软件,接下来就要下载注释用的数据库啦!

它支持非常多的物种的注释,我这里拿人类做例子咯:ftp://ftp.ensembl.org/pub/release-82/variation/VEP/

我下载的是ftp里面的82 版本: wget ftp://ftp.ensembl.org/pub/release-82/variation/VEP/homo_sapiens_refseq_vep_82_GRCh37.tar.gz

有6.1G,所以会有点耗时~

下载完毕后直接用tar –zxvf解压即可使用啦!

我安装软件的时候指定了cache目录,而不是默认的$HOME/.vep/

 

Download the archive file for your species 

Extract the archive in your cache directory. By default the VEP uses $HOME/.vep/ as the cache directory, where $HOME is your UNIX home directory. 

mv homo_sapiens_vep_84.tar.gz ~/.vep/ cd ~/.vep/tar xfz homo_sapiens_vep_84.tar.gz

Run the VEP with the--cache option

所以要把下载的6.1G数据库放在我自己的cashe目录

如果你安装VEP的时候用的默认安装参数,就需要把自己下载的6.1G文件放在  ~/.vep/ 目录下面

参考:http://davetang.org/wiki2/index.php?title=VEP

输入数据:

它支持好几种输入格式数据:

  • BED: a simple tab-delimited format containing 3-12 columns of data. The first 3 columns contain the coordinates of the feature. If available, the VEP will use the 4th column of the file as the identifier of the feature.
  • GFF: a format for describing genes and other features. If available, the VEP will use the "ID" field as the identifier of this feature.
  • GTF: treated in an identical manner to GFF.
  • VCF: a format used to describe genomic variants. The VEP will use the 3rd column of the file as the identifier.
  • bigWig: a format for storage of dense continuous data. The VEP uses the value for the given position as the "identifier". Note that bigWig files contain their own indices, and do not need to be indexed by tabix.

Any other files can be easily converted to be compatible with the VEP; the easiest format to produce is a BED-like file containing coordinates and an (optional) identifier:

其实重点就是给出你的突变的坐标即可,在哪条染色体,什么位置!

我们可以拿snpEFF里面的example文件夹里面的数据来做测试。

 

 

运行命令:

可以直接进入安装目录(VEP_ensembl/ensembl-tools-release-84/scripts/variant_effect_predictor)运行那个主程序

variant_effect_predictor.pl -i example_GRCh37.vcf \

--cache --assembly GRCh37 \

--offline --force_overwrite

或者用全路径的形式去调用这个程序

参数非常复杂,详细介绍见:http://useast.ensembl.org/info/docs/tools/vep/script/vep_options.html

一般用标准参数就好啦,而且还有一些插件,其中我比较喜欢dbNSFP and LOFTEE plugins,这也是EXAC里面用过的。

3

 

结果解读:

这个非常复杂,对结果理解了多少,就是我们对软件理解了多少。

具体大家看readme吧,注释信息太多了,按需索取:

直接看EXAC(ExAC.r0.3.sites.vep.vcf.gz)文件里面近一亿条突变记录也能慢慢理解!

 

参考:http://gemini.readthedocs.io/en/latest/content/functional_annotation.html

05

用snpEFF对vcf格式的突变数据进行注释

这个软件比较重要,尤其是对做遗传变异相关研究的,很多人做完了snp-calling后喜欢用ANNOVAR来进行注释,但是那个注释还是相对比较简单,只能得到该突变位点在基因的哪个区域,那个基因这样的信息,如果想了解更具体一点,就需要更加功能化的软件了,snpEFF就是其中的佼佼者,而且是java平台软件,非常容易使用!而且它的手册写的非常详细:http://snpeff.sourceforge.net/SnpEff_manual.html

官网是:http://snpeff.sourceforge.net/

1       889455  .       G       A       .       .        ## 假设我们的vcf文件里面记录的突变是这个,那么我们可以用snpEFF进行注释,注释得到的信息非常完全!

信息用|符号分割,所有很容易用脚本提取需要的信息

ANN=A|stop_gained|HIGH|NOC2L|ENSG00000188976|transcript|ENST00000327044|protein_coding|7/19|c.706C>T|p.Gln236*|756/2790|706/2250|236/749||,A|downstream_gene_variant|MODIFIER|NOC2L|ENSG00000188976|transcript|ENST00000487214|processed_transcript||n.*865C>T|||||351|,A|downstream_gene_variant|MODIFIER|NOC2L|ENSG00000188976|transcript|ENST00000469563|retained_intron||n.*878C>T|||||4171|,A|non_coding_exon_variant|MODIFIER|NOC2L|ENSG00000188976|transcript|ENST00000477976|retained_intron|5/17|n.2153C>T||||||;LOF=(NOC2L|ENSG00000188976|6|0.17);NMD=(NOC2L|ENSG00000188976|6|0.17)

包括突变类型是:non_coding_exon_variant

突变在各种转录本上面,在每个转录本的第几个碱基呀,哪个氨基酸的改变呀,氨基酸第几位呀!

标准突变表示形式是:

突变发生在NOC2L这个基因上面,它的ensembl 数据库ID是ENSG00000188976

 

其余的看头文件自己慢慢理解:

"Functional annotations: 'Allele | Annotation | Annotation_Impact | Gene_Name | Gene_ID | Feature_Type | Feature_ID | Transcript_BioType | Rank | HGVS.c | HGVS.p | cDNA.pos / cDNA.length | CDS.pos / CDS.length | AA.pos / AA.length | Distance | ERRORS / WARNINGS / INFO'

 

软件安装:

选择最新版软件下载:https://sourceforge.net/projects/snpeff/files/

wget https://sourceforge.net/projects/snpeff/files/snpEff_latest_core.zip

因为是java软件,unzip 解压之后就可以直接使用,当然前提是你有java平台。

1

输入数据:

首先下载用来做注释的数据库:java -jar snpEff.jar download GRCh37.75(自己选择需要的版本)

1

软件下载很快,但是数据库下载就需要一定时间啦,去喝杯咖啡吧。

然后软件本身会提供example文件,里面就是一堆各种各样的vcf数据,而且还提供了运行命令,非常简单(examples.sh) ,这些就是我们的输入数据啦!

运行命令:

运行也很简单:java -Xmx4G -jar snpEff.jar -i vcf -o vcf GRCh37.75 example.vcf > example_snpeff.vcf

指定输入输出格式都是vcf,然后指定刚才下载的必备数据库,然后输入输出文件即可!

也可以调用全路径,如果你写在脚本里面的话!

java -Xmx4G -jar path/to/snpEff/snpEff.jar \

-c path/to/snpEff/snpEff.config \

GRCh37.69 \

path/to/example.vcf > example_snpeff.vcf

 

结果解读:

这个非常复杂,对结果理解了多少,就是我们对软件理解了多少。

具体大家看readme吧,注释信息太多了,按需索取:

  1. chromosome_number_variation
  2. exon_loss_variant
  3. frameshift_variant
  4. stop_gained
  5. stop_lost
  6. start_lost
  7. splice_acceptor_variant
  8. splice_donor_variant
  9. rare_amino_acid_variant
  10. missense_variant
  11. inframe_insertion
  12. disruptive_inframe_insertion
  13. inframe_deletion
  14. disruptive_inframe_deletion
  15. 5_prime_UTR_truncation+exon_loss_variant
  16. 3_prime_UTR_truncation+exon_loss
  17. splice_branch_variant
  18. splice_region_variant
  19. splice_branch_variant
  20. stop_retained_variant
  21. initiator_codon_variant
  22. synonymous_variant
  23. initiator_codon_variant+non_canonical_start_codon
  24. stop_retained_variant
  25. coding_sequence_variant
  26. 5_prime_UTR_variant
  27. 3_prime_UTR_variant
  28. 5_prime_UTR_premature_start_codon_gain_variant
  29. upstream_gene_variant
  30. downstream_gene_variant
  31. TF_binding_site_variant
  32. regulatory_region_variant
  33. miRNA
  34. custom
  35. sequence_feature
  36. conserved_intron_variant
  37. intron_variant
  38. intragenic_variant
  39. conserved_intergenic_variant
  40. intergenic_region
  41. coding_sequence_variant
  42. non_coding_exon_variant
  43. nc_transcript_variant
  44. gene_variant
  45. chromosome

http://snpeff.sourceforge.net/VCFannotationformat_v1.0.pdf

 

14

居然可以下载千人基因组计划的所有数据bam,vcf数据

它有两个ftp站点存储所有的数据!
直接看最新版的数据,共有NA编号开头的1182个人,HG开头的1768个人!
每个人的目录下面都有 四个数据文件夹
Oct 01 2014 00:00    Directory alignment
Oct 01 2014 00:00    Directory exome_alignment
Oct 01 2014 00:00    Directory high_coverage_alignment
Oct 01 2014 00:00    Directory sequence_read
这些数据实在是太丰富了!
也可以直接看最新版的vcf文件,记录了这两千多人的所有变异位点信息!
可以直接看到所有的位点,具体到每个人在该位点是否变异!
不过它的基因型信息是通过MVNcall+SHAPEIT这个程序call出来的,具体原理见:http://www.ncbi.nlm.nih.gov/pubmed/23093610
它有两个ftp站点存储所有的数据!
直接看最新版的数据,共有NA编号开头的1182个人,HG开头的1768个人!
每个人的目录下面都有 四个数据文件夹
Oct 01 2014 00:00    Directory alignment
Oct 01 2014 00:00    Directory exome_alignment
Oct 01 2014 00:00    Directory high_coverage_alignment
Oct 01 2014 00:00    Directory sequence_read
这些数据实在是太丰富了!
也可以直接看最新版的vcf文件,记录了这两千多人的所有变异位点信息!
可以直接看到所有的位点,具体到每个人在该位点是否变异!
不过它的基因型信息是通过MVNcall+SHAPEIT这个程序call出来的,具体原理见:http://www.ncbi.nlm.nih.gov/pubmed/23093610
clipboard

 

05

使用oncotator做突变注释

功能:vcf格式突变数据进一步注释成maf格式

做过癌症数据分析的童鞋都知道,TCGA里面用maf格式来记录突变!那么maf格式的数据是如何得来的呢,我们都知道,做完snp-calling一般是得到vcf格式的突变记录数据文件,然后再用annovar或者其它蛋白结构功能影响预测软件注释一下,还远达不到maf的近100条记录。

而大名鼎鼎的broad institute就规定了maf格式的突变注释文件,他就是利用了十几个常见的已知数据库来注释我们得到的vcf突变记录,通常是对somatic的突变才注释成maf格式的数据!
大名鼎鼎的broadinstitute出品的突变注释工具:http://www.ncbi.nlm.nih.gov/pubmed/25703262
本身也是一个在线工具:
集成了下面所有的分析资源
而且还提供了API

Genomic Annotations

  • Gene, transcript, and functional consequence annotations using GENCODE for hg19.
  • Reference sequence around a variant.
  • GC content around a variant.
  • Human DNA Repair Gene annotations from Wood et al.

Protein Annotations

  • Site-specific protein annotations from UniProt.
  • Functional impact predictions from dbNSFP.

Cancer Variant Annotations

Non-Cancer Variant Annotations

因为要下载的数据有点多,我这里就不用自己的电脑测试了,安装过程也很简单的!

 

17

转载-VCF格式详解

CHROM(chromosome):染色体

POS - position:参考基因组variant碱基位置,如果是INDEL(插入缺失),位置是INDEL的第一个碱基位置

ID - identifier: variant的ID。比如在dbSNP中有该SNP的id,则会在此行给出;若没有,则用’.'表示其为一个novel variant。

REF - reference base(s):参考碱基,染色体上面的碱基,必须是ATCGN中的一个,N表示不确定碱基

ALT - alternate base(s):与参考序列比较发生突变的碱基

QUAL - quality: Phred格式(Phred_scaled)的质量值,表 示在该位点存在variant的可能性;该值越高,则variant的可能性越大;计算方法:Phred值 = -10 * log (1-p) p为variant存在的概率; 通过计算公式可以看出值为10的表示错误概率为0.1,该位点为variant的概率为90%。

FILTER - _filter status: 使用上一个QUAL值来进行过滤的话,是不够的。GATK能使用其它的方法来进行过滤,过滤结果中通过则该值为”PASS”;若variant不可靠,则该项不为”PASS”或”.”。

INFO - additional information:  这一行是variant的详细信息,具体如下:

DP-read depth:样本在这个位置的reads覆盖度。是一些reads被过滤掉后的覆盖度。          DP4:高质量测序碱基,位于REF或者ALT前后

MQ:表示覆盖序列质量的均方值RMS Mapping Quality

FQphred值关于所有样本相似的可能性

AF1 AF(Allele Frequency) 表示Allele的频率,AF1为第一个ALT allele 发生频率的可能性评估

AC1AC表示Allele(等位基因)的数目,AC1为对第一个ALT allele count的最大可能性评估

AN:AN(Allele Number) 表示Allele的总数目

IS插入缺失或部分插入缺失的reads允许的最大数量

ACAC(Allele Count) 表示该Allele的数目

G3ML 评估基因型出现的频率

HWE:chi^2基于HWE的测试p值和G3

CLR在受到或者不受限制的情况下基因型出现可能性log值

UGT:最可能不受限制的三种基因型结构

CGT:最可能受限制三种基因型的结构

PV4四种P值得误差,分别是(strand、baseQ、mapQ、tail distance bias)

INDEL:表示该位置的变异是插入缺失

PC2非参考等位基因的phred(变异的可能性)值在两个分组中大小不同

PCHI2后加权chi^2,根据p值来测试两组样本之间的联系

QCHI2:Phred scaled PCHI2.

PR置换产生的一个较小的PCHI2

QBD:Quality by Depth,测序深度对质量的影响

RPB序列的误差位置(Read Position Bias)

MDV:样本中高质量非参考序列的最大数目

VDB:Variant Distance Bias,RNA序列中过滤人工拼接序列的变异误差范围

GT样品的基因型(genotype)。两个数字中间用’/'分 开,这两个数字表示双倍体的sample的基因型。0 表示样品中有ref的allele; 1 表示样品中variant的allele; 2表示有第二个variant的allele。因此: 0/0 表示sample中该位点为纯合的,和ref一致; 0/1 表示sample中该位点为杂合的,有ref和variant两个基因型; 1/1 表示sample中该位点为纯合的,和variant一致。

GQ基因型的质量值(Genotype Quality)。Phred格式(Phred_scaled)的质量值,表示在该位点该基因型存在的可能性;该值越高,则Genotype的可能性越 大;计算方法:Phred值 = -10 * log (1-p) p为基因型存在的概率。

GL三种基因型(RR RA AA)出现的可能性,R表示参考碱基,A表示变异碱基

DV高质量的非参考碱基

SPphred的p值误差线

PL:指定的三种基因型的质量值(provieds the likelihoods of the given genotypes)。这三种指定的基因型为(0/0,0/1,1/1),这三种基因型的概率总和为1。和之前不一致,该值越大,表明为该种基因型的可能 性越小。 Phred值 = -10 * log (p) p为基因型存在的概率。

FORMAT BC1-1-base.sorted.bam这两行合起来提供了’ BC1-1-base′这个sample的基因型的信息。’ BC1-1-base′代表这该名称的样品,是由BAM文件中的@RG下的 SM 标签决定的。