Workflow - Small Variant calling and Variant Annotation¶

École de bioinformatique IFB-INSERM 2024

  • Vivien DESHAIES - AP-HP

III - Variant Calling / TP¶

GATK HaplotypeCaller :

  • GATK (Genome Analysis ToolKit) est une suite d’outils développée par le Broad Institute
  • Bonne documentation (Best Practices)
  • Permet la gestion d’analyse de plusieurs échantillons (format gVCF)
  • Comporte une étape de réassemblage et réalignement local des indel.
  • Algorithme bayésien (modèles statistiques pour estimer la probabilité de chaque génotype possible, en prenant en compte les différents biais pouvant introduire du bruit dans les données)
In [ ]:
module load gatk4/4.2.3.0

GATK documentations### 1. Single-sample variant calling GATK HaplotypeCaller avec sortie VCF

1. Indexation du génome¶

Créaction du dictionnaire de sequence (GATK)

In [ ]:
gatk CreateSequenceDictionary \
    --REFERENCE genome/Bos_taurus.UMD3.1.dna.toplevel.6.fa \
    --OUTPUT    genome/Bos_taurus.UMD3.1.dna.toplevel.6.dict

Créaction de l'index fai

In [ ]:
module load samtools/1.13
In [ ]:
samtools faidx genome/Bos_taurus.UMD3.1.dna.toplevel.6.fa

1. Single-sample variant calling¶

GATK HaplotypeCaller avec sortie VCF

In [ ]:
#gatk HaplotypeCaller --help
In [ ]:
mkdir -p GATK/vcf

add image path !!!!

In [ ]:
gatk HaplotypeCaller --java-options '-Xmx8G' \
    --input SRR1262731_extract.sort.md.filt.onTarget.bam \
    --reference genome/Bos_taurus.UMD3.1.dna.toplevel.6.fa \
    --min-base-quality-score 18  \
    --minimum-mapping-quality 30 \
    --emit-ref-confidence "NONE" \
    --output GATK/vcf/SRR1262731_extract_GATK.vcf \
    --intervals additionnal_data/QTL_BT6.bed
In [ ]:
ls -ltrh GATK/vcf/
In [ ]:
head GATK/vcf/SRR1262731_extract_GATK.vcf

2. Multi-sample variant calling¶

GATK HaplotypeCaller en mode GVCF

En 3 étapes (=> 3 outils):

  • HaplotypeCaller en mode gvcf
  • CombineGVCFs
  • GenotypeGVCFs

add image!!!

a) HaplotypeCaller en mode gvcf¶

Détection de variants GATK avec sortie gVCF

In [ ]:
# Création d’un répertoire pour l’appel des variants
mkdir -p GATK/gvcf
In [ ]:
gatk HaplotypeCaller --java-options '-Xmx8G' \
    --input SRR1262731_extract.sort.md.filt.onTarget.bam \
    --reference genome/Bos_taurus.UMD3.1.dna.toplevel.6.fa \
    --min-base-quality-score 18  \
    --minimum-mapping-quality 30 \
    --emit-ref-confidence "GVCF" \
    --output GATK/gvcf/SRR1262731_extract_GATK.g.vcf \
    --intervals additionnal_data/QTL_BT6.bed
In [ ]:
ls -ltrh GATK/gvcf/
In [ ]:
grep '^##GVCF' GATK/gvcf/SRR1262731_extract_GATK.g.vcf
In [ ]:
grep -v '^#' GATK/gvcf/SRR1262731_extract_GATK.g.vcf | head

Documentation gvcf : https://gatk.broadinstitute.org/hc/en-us/articles/360035531812-GVCF-Genomic-Variant-Call-Format

b) CombineGVCFs¶

Fusion des fichiers gVCFs en un seul gVCF

In [ ]:
gatk CombineGVCFs --java-options '-Xmx8G' \
    --variant   GATK/gvcf/SRR1262731_extract_GATK.g.vcf \
    --variant   additionnal_data/SRR1205992_extract_GATK.g.vcf \
    --variant   additionnal_data/SRR1205973_extract_GATK.g.vcf \
    --reference genome/Bos_taurus.UMD3.1.dna.toplevel.6.fa \
    --intervals additionnal_data/QTL_BT6.bed \
    --output    GATK/gvcf/pool_GATK.g.vcf

b) GenotypeGVCFs¶

Détection de variants simultanée sur les 3 échantillons du gVCF

In [ ]:
gatk GenotypeGVCFs --java-options '-Xmx8G' \
    --variant   GATK/gvcf/pool_GATK.g.vcf \
    --reference genome/Bos_taurus.UMD3.1.dna.toplevel.6.fa \
    --output    GATK/vcf/pool_GATK.vcf
In [ ]:
grep -v '^#' GATK/vcf/pool_GATK.vcf | head -n1

IV - Variant Annotation / TP (SnpEff)¶

In [ ]:
module load snpeff/4.3.1t
In [ ]:
snpEff databases | grep Bos_taurus
In [ ]:
#snpEff build | head

http://pcingola.github.io/SnpEff/snpeff/inputoutput/#eff-field-vcf-output-files

1. Création de la base de données SnpEff¶

In [ ]:
echo "BosTaurus.genome" > genome/snpeff.config # <genome_name>.genome
In [ ]:
mkdir -p genome/BosTaurus
In [ ]:
cp genome/Bos_taurus.UMD3.1.dna.toplevel.6.fa genome/BosTaurus/sequences.fa
In [ ]:
cp genome/Bos_taurus.UMD3.1.93.chromosome.6.gff3 genome/BosTaurus/genes.gff
In [ ]:
echo -e "BosTaurus\nSnpEff4.3t" > genome/BosTaurus.db
In [ ]:
snpEff build -c genome/snpeff.config -gff3 -v BosTaurus -dataDir . 2>log/snpEff_build.out

2. Annotation avec notre base de données¶

In [ ]:
# usage : snpEff eff  -c <snpeff_config>  <genome_version> <input_vcf>
snpEff eff -c genome/snpeff.config \
           -dataDir . \
           BosTaurus \
           -s snpeff_res.html \
           GATK/vcf/pool_GATK.vcf > GATK.annot.vcf
In [ ]:
grep -v '^#' GATK.annot.vcf | head -n1

3. Filtre (SnpSift)¶

In [ ]:
module load snpsift/4.3.1t
SnpSift filter -h | cat # affiche l’aide

Documentation SnpSift filter : https://pcingola.github.io/SnpEff/snpsift/filter/

Garder les variants codant qui ne sont pas des synonymes :

In [ ]:
# usage : SnpSift filter [options] 'expression' <input.vcf>
SnpSift filter -Xmx8G \
    "(ANN[*].EFFECT != 'synonymous_variant') && (ANN[*].BIOTYPE = 'protein_coding')" \
    GATK.annot.vcf > GATK.annot.coding.nosyn.vcf
In [ ]:
grep -v '^#' GATK.annot.coding.nosyn.vcf | grep 'synonymous_variant'

Sélectionner notre variant d’intérêt parmi les variants hétérozygotes ayant un impact (missense)

In [ ]:
SnpSift filter -Xmx8G \
 "ANN[*].EFFECT = 'missense_variant' & isHet( GEN[2] ) & isVariant( GEN[2] ) \
 & isRef( GEN[0] ) & isRef( GEN[1] ) " \
 GATK.annot.coding.nosyn.vcf > GATK.interest_variants.vcf
In [ ]:
grep -v '^#' GATK.interest_variants.vcf