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
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
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
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