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

 

十二 29

用Mutation-Assessor软件来看突变位点对基因或者蛋白功能的影响

这是一个在线工具,非常好用,对snp位点进行注释,看看该突变是否影响蛋白功能,一定要收藏!!!

官网:http://mutationassessor.org/

也应该是有standalone版本,我没有去找,不过网页版就很好用了,只需要复制粘贴进去自己想分析的数据,按照一定的格式即可,比如:
clipboard
该软件就能分析出该突变位点发生在BRCA2这个基因上面,对氨基酸的改变也能写出来,对蛋白功能的改变等选项都是可以自由定制化的。
输入数据非常广泛:
The server accepts list of variants, one variant per line, plus optional text describing your variants,
in genomic coordinates, "+" strand assumed :
<genome build>,<chromosome>,<position>,<reference allele>,<substituted allele>
Genome build is optional (build 18 assumed), accepted values: 'hg18' and 'hg19'
Examples:

hg19,13,32912555,G,T   BRCA2
hg18,7,55178574,G,A   GBM
7,55178574,G,A   GBM

or in protein space: <protein ID> <variant> <text>, where <protein ID> can be :

1. Uniprot protein accession (i.e. EGFR_HUMAN)
2. NCBI Refseq protein ID (i.e. NP_005219)

examples:

EGFR_HUMAN R521K
EGFR_HUMAN R98Q Polymorphism
EGFR_HUMAN G719D disease
NP_000537 G356A
NP_000537 G360A dbSNP:rs35993958
NP_000537 S46A Abolishes phosphorylation

ID types can be mixed in one list in any way.

 

16

3000多份水稻全基因组测序数据共享-主要是突变数据

感觉最近接触的生物信息学知识越多,越对大数据时代的到来更有同感了。现在的研究者,其实很多都可以自己在家里做了,大量的数据基本都是公开的, 但是一个人闭门造车成就真的有限,与他人交流的思想碰撞还是蛮重要的。

这里面列出了3000多份水稻全基因组测序数据,都共享在亚马逊云上面,是全基因组的双端测序数据,共3,024个水稻数据,比对到了五种不同的水稻参考基因组上面,而且主要是用GATK来找差异基因的。
而且,数据收集者还给出了一个snp calling的标准流程
我以前也是用这样的流程
SNP Pipeline Commands

1. Index the reference genome using bwa index

   /software/bwa-0.7.10/bwa index /reference/japonica/reference.fa

2. Align the paired reads to reference genome using bwa mem. 
   Note: Specify the number of threads or processes to use using the -t parameter. The possible number of threads depends on the machine where the command will run.

   /software/bwa-0.7.10/bwa mem -M -t 8 /reference/japonica/reference.fa /reads/filename_1.fq.gz /reads/filename_2.fq.gz > /output/filename.sam

3. Sort SAM file and output as BAM file

   java -Xmx8g -jar /software/picard-tools-1.119/SortSam.jar INPUT=/output/filename.sam OUTPUT=/output/filename.sorted.bam VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=TRUE

4. Fix mate information

   java -Xmx8g -jar /software/picard-tools-1.119/FixMateInformation.jar INPUT=/output/filename.sorted.bam OUTPUT=/output/filename.fxmt.bam SO=coordinate VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=TRUE

5. Mark duplicate reads

   java -Xmx8g -jar /software/picard-tools-1.119/MarkDuplicates.jar INPUT=/output/filename.fxmt.bam OUTPUT=/output/filename.mkdup.bam METRICS_FILE=/output/filename.metrics VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=TRUE MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000

6. Add or replace read groups

   java -Xmx8g -jar /software/picard-tools-1.119/AddOrReplaceReadGroups.jar INPUT=/output/filename.mkdup.bam OUTPUT=/output/filename.addrep.bam RGID=readname PL=Illumina SM=readname CN=BGI VALIDATION_STRINGENCY=LENIENT SO=coordinate CREATE_INDEX=TRUE

7. Create index and dictionary for reference genome

   /software/samtools-1.0/samtools faidx /reference/japonica/reference.fa
   
   java -Xmx8g -jar /software/picard-tools-1.119/CreateSequenceDictionary.jar REFERENCE=/reference/japonica/reference.fa OUTPUT=/reference/reference.dict

8. Realign Target 

   java -Xmx8g -jar /software/GenomeAnalysisTK-3.2-2/GenomeAnalysisTK.jar -T RealignerTargetCreator -I /output/filename.addrep.bam -R /reference/japonica/reference.fa -o /output/filename.intervals -fixMisencodedQuals -nt 8

9. Indel Realigner

   java -Xmx8g -jar /software/GenomeAnalysisTK-3.2-2/GenomeAnalysisTK.jar -T IndelRealigner -fixMisencodedQuals -I /output/filename.addrep.bam -R /reference/japonica/reference.fa -targetIntervals /output/filename.intervals -o /output/filename.realn.bam 

10. Merge individual BAM files if there are multiple read pairs per sample

   /software/samtools-1.0/samtools merge /output/filename.merged.bam /output/*.realn.bam

11. Call SNPs using Unified Genotyper

   java -Xmx8g -jar /software/GenomeAnalysisTK-3.2-2/GenomeAnalysisTK.jar -T UnifiedGenotyper -R /reference/japonica/reference.fa -I /output/filename.merged.bam -o filename.merged.vcf -glm BOTH -mbq 20 --genotyping_mode DISCOVERY -out_mode EMIT_ALL_SITES
09

对vcf突变数据与公开发表的进行比对

当我们对NGS数据call了snp之后一般会输出成vcf格式的数据,一行代表一个突变,例如
20      2451451 .       G       T       1939.77 .
AC=1;AF=0.500;AN=2;BaseQRankSum=-10.134;DP=239;Dels=0.00;FS=2.276;HaplotypeScore=0.0000;MLEAC=1;MLEAF=0.500;MQ=60.00;MQ0=0;MQRankSum=-0.258;QD=8.12;ReadPosRankSum=0.823;SOR=0.870
GT:AD:DP:GQ:PL  0/1:150,89:239:99:1968,0,3874
#前几列记录着该突变发生在第几号染色体以及该染色体的哪个坐标,我们的参考基因组在该位点是什么碱基,我们测到的突变成了什么碱基。
最后两列是测序深度以及正负测序深度,或者ref和allele的测序深度。
只有第8列是最复杂的,可以有高达几百个数据信息,取决于我们用什么样的软件来call的snp,以及call了snp之后用什么样的软件做的注释。
接下来我们还需要探究我们找到的突变是否在其它以及公开发表的数据库里面被找到过,所以可以下载非常多的公共数据库进行比对,我所见过的有一下一些,估计完全下载有0.5T
dbsnp144 (这个是ncbi提供的最权威的啦)
cgi69
ExAC.vcf.gz(这个是broadinstitute提供的外显子联盟)
Cosmic_v73.ann.vcf.gz (这个是癌症突变信息集)
finalTCGA.vcf.gz (TCGA计划也是癌症相关的)
icgc.vcf.gz
dbNSFP2.6vcf
SCLP.ann.vcf.gz
CCLE.ann.vcf.gz
ESP6500-SIv2.vcf.gz (Variants from the Exome Sequencing Project (ESP))
adni-sum
safs-sum.indel.vcf.gz
gonl.vcf.gz
ssm.vcf.gz
ssi.vcf.gz
uk10k.vcf.gz
1000g-ph3v5.gff.gz  (千人基因组计划)
gwasCatalog.gff.gz  \
phewascatalog.gff.gz  \
dbgap-gwas.gff.gz  \
interproDomain.gff.gz \
clinvar.gff.gz \
RegulomeDB.gff.gz \
CancerGAMAdb.gff.gz \