标签归档:VCF

将VCF文件中的突变拆分成SNP和INDEL

VCFTOOLS

得到SNP

vcftools --vcf X.vcf --remove-indels --out X.snps --recode --recode-INFO-all

得到INDEL

vcftools --vcf X.vcf --keep-only-indels --out X.indel --recode --recode-INFO-all

GATK

得到SNP

java -jar GenomeAnalysisTK.jar \
    -T SelectVariants \
    -R reference.fasta \
    -V X.vcf   \
    -selectType SNP \
    -o X.snps.vcf

得到INDEL

java -jar GenomeAnalysisTK.jar \
    -T SelectVariants \
    -R reference.fasta \
    -V X.vcf  \
    -selectType INDEL \
    -o X.indel.vcf

本着有轮子不造轮子的原则,可以用VCFTOOLS和GATK来实现,当然如果想自己拆分的话,可以根据VCF中是否有SNP和INDEL的tag标签,或者根据ALT和REF中的碱基长度是否一致来实现拆分。

参考:
https://www.biostars.org/p/48204/
https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_variantutils_SelectVariants.php

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################

BCFTOOLS支持的表达式

翻译https://samtools.github.io/bcftools/bcftools-man.html#expressions

有效的表达式可能包含

  • 数字常量,字符常量,文件名
  • 1, 1.0, 1e-4
    "String"
    @file_name
  • 算数运算符:
  • +,*,-,/
  • 比较运算符:
  • == (同=), >, >=, < =, <, !=
  • 正则运算符"~"及其非值"!~"。大小写敏感,加"/i"这大小写不敏感:
  • INFO/HAYSTACK ~ "needle"
    INFO/HAYSTACK ~ "NEEDless/i"
  • 括号:
  • (, )
  • 逻辑运算符:
  • && (同&), ||,  |
  • VCF/BCF中的INFO tags, FORMAT tags, column names:
  • INFO tags: INFO/DP or DP
    FORMAT tags: FORMAT/DV, FMT/DV, or DV
    column names:  FILTER, QUAL, ID, POS, REF, ALT[0]
  • 1 (或0) 检测flag的存在 (或不存在) :
  • FlagA=1 && FlagB=0

    继续阅读

    VCF文件中的原始突变过滤–filter raw variants in vcf

    Hard filter突变的传统过滤方式

    此时VCF文件中的突变,与刚开始下机得到的FASTQ文件类似,称为raw data。此时的突变集合中,有很多假阳性突变,这些突变需要在突变分析之前过滤掉。
    传统的过滤方式,直接根据每个突变的注释信息,进行过滤。最直接和最常见的是根据DP标签过滤,即根据该突变位点的测序深度进行过滤。通常,深度越低,支持该突变的reads数目越少,该突变越不可信。还可以根据前面提到的QUAL质量分值进行过滤,分值越低越不可信。Forward reads和Reverse reads的比例。通过,设定一定的阈值,看这些注释信息是高于还是低于该阈值。

    这种直接根据突变信息进行过滤的方式,GATK成为hard filter。GATK常用的hard filter参数还有DP < 10,QD < 2.0,FS > 60.0,MQ < 10.0,MQRankSum < -12.5,ReadPosRankSum < -8.0等方式。这些阈值通过GATK的VariantFiltration工具进行过滤,突变满足其中一条,就会被过滤。 这种过滤方式直接根据特定阈值就将突变过滤掉,考虑的维度较少,真实突变可能因为某一注释没有到达阈值而被错误的过滤掉。如果为了保留这些真实的突变,而放松阈值,又可能同样将假阳性突变保留。 Hard filter常用于panel测序,panel测序得到的突变位点较少,不足以通过机器学习的方式进行过滤。

    VQSR突变质量分值再校准

    虽然这步叫做Varaint Quality Score Recalibrate,但该步并没有再校准突变质量分值,而是重新生成了一个叫做VQSLOD (for variant quality score log-odds)的分值,且认为该分值是被充分校准的。

    VQSR不同于传统的hard filter。VQSR整合突变的多种维度的信息,如下图,通过机器学习的方式,得到好的突变good variant和坏的突变bad variant的特征范式,对突变进行过滤,从而挑选出聚簇在一起的好的突变。需要注意的是,VQSR只利用INFO的注释信息,而没有利用样本特异性的SAMPLE信息。

    通常用于训练的数据集,是通过分析得到的raw variant data数据集与hapmap、omni、1000Genome等(hapma和omni中的突变通常被认为是高质量的、验证过的,1000 Genome中的突变可信度次之)进行比较,得到raw variant data数据集中高可信的突变集合,集合中的突变注释信息将被用于机器学习训练模型。同理,低质量假阳性突变的数据也同样。通过机器学习,会对训练数据集高可信的突变的分布进行建模和聚类。同时,也对所有突变进行建模,计算VQSLOD值。通常离聚类形成的簇核心近的突变,VQSLOD值越高,越可信。

    在最终过滤突变的时候,GATK VQSR并没有一个直接的阈值进行过滤,而是通过引入sensitivity灵敏度的概念。比如过滤的阈值是99%,则意味着训练数据集中有99%的突变也位于hapmap中,此时的VQSLOD值就是用于突变过滤的阈值。VQSLOD值高于该阈值的突变会被标为PASS,低于的则被过滤,并在FILTER列标记。
    因为SNP和INDEL的特征不同,VQSR通常将SNP和INDEL分开进行过滤。

    #####################################################################
    #版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
    #Author: Jason
    #####################################################################

    oncotator对VCF进行注释,并生成MAF格式文件

    MAF格式Mutation Annotation Format (MAF) ,是TCGA组织对突变进行注释的格式。一些和癌症分析相关的软件,经常要求MAF格式文件作为输入。而现在经过GATK或samtools检测出突变的格式一般为VCF格式,的注释软件,即使经过SNPEff和annovar注释(当然还有VEP),结果依然为VCF格式或者tab分割的文件等。

    MAF中每一列是一种注释信息,由于包含的注释信息太多(详见格式),单纯的通过写脚本转换SNPEff或者annovar的注释文件,会变得非常麻烦而且考虑的问题可能不完全(有人实现过,通过Ensembl的VEP对VCF注释,然后转换,可以在github上搜索到)。

    这里介绍注释软件oncotator,可以注释VCF文件,并直接生成MAF格式,相当于将VCF格式转换成MAF格式。Broad institute开发的,用起来放心哈。

    首先安装python virtualenv

    virtualenv可以创建一个虚拟的python环境,在虚拟环境中安装的包,不影响机器已经安装的python包。

    curl -O https://pypi.python.org/packages/source/v/virtualenv/virtualenv-14.0.6.tar.gz
    tar -xvzf virtualenv-14.0.6.tar.gz 
    cd virtualenv-14.0.6/
    python setup.py install --prefix=~/pypackages/
    cd ..
    rm -rf virtualenv-14.0.6*

    安装oncotator依赖

    git clone  https://github.com/broadinstitute/oncotator.git
    cd oncotator
    mkdir env
    bash scripts/create_oncotator_venv.sh -e ./env

    安装oncotator

    source ./env/bin/activate
    python setup.py install
    deactivate

    下载注释数据文件

    wget -c http://www.broadinstitute.org/~lichtens/oncobeta/oncotator_v1_ds_Jan262015.tar.gz
    tar -xvzf  oncotator_v1_ds_Jan262015.tar.gz

    注释内容如下,包括基因组注释,蛋白注释,癌症突变注释等,用的基因组是GENCODE的hg19版本。
    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

    Observed cancer mutation frequency annotations from COSMIC.
    Cancer gene and mutation annotations from the Cancer GenCensus.
    Overlapping mutations from the Cancer Cell Line Encyclopedia.
    Cancer gene annotations from the Familial Cancer Database.
    Cancer variant annotations from ClinVar.

    Non-Cancer Variant Annotations

    Common SNP annotations from dbSNP.
    Variant annotations from 1000 Genomes.
    Variant annotations from NHLBI GO Exome Sequencing Project (ESP).

    运行oncotator

    #进入虚拟环境env
    source ./env/bin/activate
    # WGS = whole genome sequencing, WXS = whole exome sequencing
    # 通过annotate-manual设置MAF注释,比如样本barcode,样本来源等,至于氨基酸突变等,会自动进行注释
    oncotator --input_format=VCF --db-dir=./oncotator_v1_ds_Jan262014 --output_format=TCGAMAF --tx-mode=EFFECT --collapse-number-annotations --annotate-default='Sequencer:Illumina HiSeq' --annotate-manual='Sequence_Source:WGS' --annotate-manual='Tumor_Sample_Barcode:Tumor'  --annotate-manual='Matched_Norm_Sample_Barcode:normal'  input.vcf output.maf hg19
    deactivate

    结果

    继续阅读

    Sort VCF by Chr and Pos根据染色体位置对VCF进行排序

    对VCF文件中的突变按照染色体和位置进行排序,下面是本人的总结,其中利用bash命令的方法不依赖其他的工具或包。htslib前文中也提到过。

    1, Use bash

    bash raw.vcf

    chr_order="chrMnchr1nchr2nchr3nchr4nchr5nchr6nchr7nchr8nchr9nchr10nchr11nchr12nchr13nchr14nchr15nchr16nchr17nchr18nchr19nchr20nchr21nchr22nchrXnchrY"
    
    cat "$1" | grep "^#" > .header.vcf
    cat "$1" | grep -v "^#" | sort -k1,1 -k2,2n > .pre.sorted.vcf
    echo -e $chr_order | while read line
    do
        cat .pre.sorted.vcf | grep "^$line"$'t' >> .header.vcf
    done
    mv .header.vcf  sorted.vcf && rm .header.vcf .pre.sorted.vcf

    2,Use awk and sed

    (awk ‘/^#/{print}!/^#/{exit}’ raw.vcf;sed ‘/^#/d’ raw.vcf|awk -F”\t” ‘($1~/^[0-9]+$/){sub(“^chr”,””,$0);print $0}’|sort -k1,1n -k2,2n|awk ‘{print “chr”$0}’ ;sed ‘/^#/d’ raw.vcf| awk -F”\t” ‘($1!~/^[0-9]+$/){sub(“^chr”,””,$0);print $0}’|sort -k1,1d -k2,2n|awk ‘{print “chr”$0}’) > sort.vcf

    3, Use Picard

    Sorts one or more VCF files. This tool sorts the records in VCF files according to the order of the contigs in the header/sequence dictionary and then by coordinate. It can accept an external sequence dictionary.

    java -jar picard.jar SortVcf I=unsort.vcf O=sorted.vcf

    4,Use vcf-sort (in vcftools)

    cat file.vcf  |  vcf-sort > sorted.vcf

    继续阅读