clinvar数据库重新解读

在很久以前的我直播基因组活动，我提到过这个数据库： 【直播】我的基因组67：clinvar数据库

不过，那个时候遗传背景知识不够，其实并没有很好的理解它，现在有机会重新学习一下，可以使用以下代码下载并且注释到clinvar数据库

wget https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/clinvar_20180429.vcf.gz
wget https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/clinvar_20180429.vcf.gz.tbi
java -jar  ~/biosoft/SnpEff/snpEff/SnpSift.jar  annotate clinvar_20180429.vcf.gz merge_snpeff.vcf >merge_clinvar.vcf

得到的注释信息的描述是：

##SnpSiftCmd="SnpSift annotate clinvar_20180429.vcf.gz merge_snpeff.vcf"
##INFO=<ID=CLNDISDBINCL,Number=.,Type=String,Description="For included Variant: Tag-value pairs of disease database name and identifier, e.g. OMIM:NNNNNN">
##INFO=<ID=DBVARID,Number=.,Type=String,Description="nsv accessions from dbVar for the variant">
##INFO=<ID=CLNSIGCONF,Number=.,Type=String,Description="Conflicting clinical significance for this single variant">
##INFO=<ID=AF_TGP,Number=1,Type=Float,Description="allele frequencies from TGP">
##INFO=<ID=MC,Number=.,Type=String,Description="comma separated list of molecular consequence in the form of Sequence Ontology ID|molecular_consequence">
##INFO=<ID=AF_EXAC,Number=1,Type=Float,Description="allele frequencies from ExAC">
##INFO=<ID=CLNDN,Number=.,Type=String,Description="ClinVar's preferred disease name for the concept specified by disease identifiers in CLNDISDB">
##INFO=<ID=ORIGIN,Number=.,Type=String,Description="Allele origin. One or more of the following values may be added: 0 - unknown; 1 - germline; 2 - somatic; 4
 - inherited; 8 - paternal; 16 - maternal; 32 - de-novo; 64 - biparental; 128 - uniparental; 256 - not-tested; 512 - tested-inconclusive; 1073741824 - other">
##INFO=<ID=CLNVI,Number=.,Type=String,Description="the variant's clinical sources reported as tag-value pairs of database and variant identifier">
##INFO=<ID=CLNVC,Number=1,Type=String,Description="Variant type">
##INFO=<ID=AF_ESP,Number=1,Type=Float,Description="allele frequencies from GO-ESP">
##INFO=<ID=CLNHGVS,Number=.,Type=String,Description="Top-level (primary assembly, alt, or patch) HGVS expression.">
##INFO=<ID=CLNDNINCL,Number=.,Type=String,Description="For included Variant : ClinVar's preferred disease name for the concept specified by disease identifier
s in CLNDISDB">
##INFO=<ID=CLNVCSO,Number=1,Type=String,Description="Sequence Ontology id for variant type">
##INFO=<ID=CLNSIG,Number=.,Type=String,Description="Clinical significance for this single variant">
##INFO=<ID=CLNDISDB,Number=.,Type=String,Description="Tag-value pairs of disease database name and identifier, e.g. OMIM:NNNNNN">
##INFO=<ID=CLNREVSTAT,Number=.,Type=String,Description="ClinVar review status for the Variation ID">
##INFO=<ID=CLNSIGINCL,Number=.,Type=String,Description="Clinical significance for a haplotype or genotype that includes this variant. Reported as pairs of VariationID:clinical significance.">
##INFO=<ID=ALLELEID,Number=1,Type=Integer,Description="the ClinVar Allele ID">

记录最多的基因是

zcat clinvar_20180429.vcf.gz|perl -alne '{/GENEINFO=(.*?):/;print $1 if $1}'|sort |uniq -c |sort -k 1,1nr >gene.clinvar.freq
head gene.clinvar.freq
   9800 TTN
   9038 BRCA2
   6311 BRCA1
   4987 ATM
   4498 APC
   3312 TSC2
   3114 MSH6
   2888 NF1
   2653 MSH2
   2223 LDLR

可以看到BRCA1/2遥遥领先呀！！！

记录的致病情况最多的基因是

zcat clinvar_20180429.vcf.gz|perl -alne '{/GENEINFO=(.*?):/;$g=$1;/CLNSIG=(.*?);/;print "$g\t$1" if $1}'| sort |uniq -c |sort -k 1,1nr >gene_sign.clinvar.freq
head gene_sign.clinvar.freq
   5218 TTN Uncertain_significance
   3293 BRCA2   Uncertain_significance
   2665 BRCA2   Pathogenic
   2523 ATM Uncertain_significance
   2162 BRCA1   Pathogenic
   1897 APC Uncertain_significance
   1887 TTN Likely_benign
   1817 BRCA1   Uncertain_significance
   1651 MSH6    Uncertain_significance
   1465 BRCA2   Likely_benign
grep Pathogenic gene_sign.clinvar.freq|head
   2665 BRCA2   Pathogenic
   2162 BRCA1   Pathogenic
    691 LDLR    Pathogenic
    594 MSH2    Pathogenic
    571 MLH1    Pathogenic
    544 COL4A5  Pathogenic
    536 NF1 Pathogenic
    469 DMD Pathogenic
    467 MSH6    Pathogenic
    466 APC Pathogenic
    

仍然是BRCA1/2遥遥领，有趣的是TTN基因的致病位点压根就排不上号，其位点虽然在clinVar数据库收录多，但意义不明确，或者并不致病。

看看是否有被专家审核

zcat clinvar_20180429.vcf.gz|perl -alne '{/GENEINFO=(.*?):/;$g=$1;/CLNSIG=(.*?);/;$t=$1;/CLNREVSTAT=(.*?);/;print "$g\t$t\t$1" if $1}'| sort |uniq -c |sort -k 1,1nr >gene_sign_review.clinvar.freq
[jianmingzeng@jade anno]$ head gene_sign_review.clinvar.freq
   4200 TTN Uncertain_significance  criteria_provided,_single_submitter
   2065 BRCA2   Pathogenic  reviewed_by_expert_panel
   1915 BRCA2   Uncertain_significance  criteria_provided,_single_submitter
   1666 BRCA1   Pathogenic  reviewed_by_expert_panel
   1639 TTN Likely_benign   criteria_provided,_single_submitter
   1597 ATM Uncertain_significance  criteria_provided,_single_submitter
   1251 APC Uncertain_significance  criteria_provided,_single_submitter
   1200 TTN Conflicting_interpretations_of_pathogenicity    criteria_provided,_conflicting_interpretations
   1174 TSC2    not_provided    no_assertion_provided

可以看到之前的Pathogenic的BRCA2突变记录是2665个，但是现在增加了 reviewed_by_expert_panel 条件后下降到了 2065个，还有600个是因为：

 grep BRCA2 gene_sign_review.clinvar.freq|grep Pathogenic
   2065 BRCA2   Pathogenic  reviewed_by_expert_panel
    436 BRCA2   Pathogenic  criteria_provided,_single_submitter
     94 BRCA2   Pathogenic  criteria_provided,_multiple_submitters,_no_conflicts
     70 BRCA2   Pathogenic  no_assertion_criteria_provided 

距离我写完这个教程，clinVar又更新了！！！！

使用annovar注释自己的vcf文件到clinVar数据库

annovar的安装及使用教程自行查看我在生信菜鸟团的博客或者生信技能树论坛分享的。

dir=/home/jianmingzeng/biosoft/ANNOVAR/annovar
db=$dir/humandb/ 
ls $db

wget https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/archive_2.0/2018/clinvar_20180603.vcf.gz 

perl $dir/annotate_variation.pl  --downdb --webfrom annovar --buildver hg38 
clinvar_20180603 $db
# perl $bin  --downdb --webfrom annovar --buildver hg38 gnomad_genome $db
mkdir annovar_results 
 $dir/convert2annovar.pl -format vcf4old highQ.vcf  1> highQ.avinput  2>/dev/null 
perl $dir/annotate_variation.pl  -buildver hg38 -filter -dbtype clinvar_20180603  --outfile annovar_results/highQ_clinvar  highQ.avinput  $db

perl $dir/table_annovar.pl   \
-buildver hg38 \
highQ.avinput  $db \
-out test \
-remove -protocol \
refGene,clinvar_20170905 \
-operation g,r \
-nastring NA 

注释结果如下：

981K Aug 13 16:52 highQ_clinvar.hg38_clinvar_20180603_dropped
11M Aug 13 16:52 highQ_clinvar.hg38_clinvar_20180603_filtered
105K Aug 13 16:52 highQ_clinvar.invalid_input
1.1K Aug 13 16:52 highQ_clinvar.log

然后我发现 clinvar_20180603 这个版本的数据库其实是有问题的。

有趣的是使用annovar下载的clinVar数据库与其NCBI官网的文件对snp的书写方式其实是不一样的.

比如官网下载文件如下：

13  32340301    266906  TG  T   .   .   ALLELEID=261321;CLNDISDB=MedGen:C2675520,OMIM:612555;CLNDN=Breast-ovarian_cancer,_familial_2;CLNHGVS=NC_000013.11:g.32340303delG;CLNREVSTAT=reviewed_by_expert_panel;CLNSIG=Pathogenic;CLNVC=Deletion;CLNVCSO=SO:0000159;GENEINFO=BRCA2:675;MC=SO:0001589|frameshift_variant;ORIGIN=1;RS=886040621

但是annovar下载的如下:

13  32340301    32340301    T   -   24364   Neoplasm_of_the_breast|Hereditary_cancer-predisposing_syndrome|Familial_cancer_of_breast|Hereditary_breast_and_ovarian_cancer_syndrome|Fanconi_anemia\x2c_complementation_group_D1|Breast-ovarian_cancer\x2c_familial_2|Pancreatic_cancer_2|not_specified|BRCA2-Related_Disorders|not_provided  Human_Phenotype_Ontology:HP:0100013\x2cMeSH:D001943\x2cMedGen:C1458155\x2cOrphanet:ORPHA180250\x2cSNOMED_CT:126926005|MedGen:C0027672\x2cSNOMED_CT:699346009|MedGen:C0346153\x2cOMIM:114480\x2cOrphanet:ORPHA227535\x2cSNOMED_CT:254843006|MedGen:C0677776\x2cOrphanet:ORPHA145|MedGen:C1838457\x2cOMIM:605724\x2cOrphanet:ORPHA319462|MedGen:C2675520\x2cOMIM:612555|MedGen:C3150546\x2cOMIM:613347|MedGen:CN169374|MedGen:CN239275|MedGen:CN517202    reviewed_by_expert_panel    Pathogenic