Workflow - Structural Variant Calling¶
École de bioinformatique IFB-INSERM 2024
- Gabriele Adam - INRAe
V - Variations structurales (SV) / TP¶
1. detection de SV à partir de données Illumina short-reads¶
En premier, nous allons utiliser Delly (papier ici doi: 10.1093/bioinformatics/bts378, code là: https://github.com/dellytools/delly). Pour commencer nous travailllons sur données short-reads.
module load delly/0.8.3
Comment connaître la commande à utiliser pour lancer delly ?
delly | cat
Et plus spécifiquement, comment lancer delly pour faire du variant calling ?
delly call | cat
a) Indexation¶
Comme vu dans les cours précédents, nous allons commencerpar indexer notre alignement et notre référence.
Load samtools:
module load samtools/1.13
samtools faidx Zymoseptoria_tritici.fa
samtools index mapping_illumina_chr10_500kb.bam
Et maintenant, utilisons delly pour le variant calling:
delly call -g Zymoseptoria_tritici.fa -o SV_calling_illumina.bcf mapping_illumina_chr10_500kb.bam
Si je veux regarder ce résultat:
head SV_calling_illumina.bcf
Ici le format n'est pas VCF mais BCF. Il faut donc le modifier.
SV_calling_illumina.bcf may be a binary file. See it anyway? n
# Conversion en fichier vcf
module load bcftools/1.10.2
bcftools view SV_calling_illumina.bcf > SV_calling_illumina.vcf
Et maintenant, nous pouvons visualiser nos variants.
Etape 1: Regarder seulement le header du fichier:
grep "^#" SV_calling_illumina.vcf
Etape 2: regarder seulement les 3 premiers variants trouvés:
grep -v "^#" SV_calling_illumina.vcf | head -n3
b) comptage du nombre de SVs¶
Combien de variants sont présents dans le fichier?
# Combien de variants ?
grep -v -c "^#" SV_calling_illumina.vcf
Combien de variants sont de bonne qualité?
# Combien de variants de bonne qualité ?
grep -v "^#" SV_calling_illumina.vcf | grep -v -c "LowQual"
Combien de variants de bonne qualité sont de type Deletion ?
grep -v "^#" SV_calling_illumina.vcf | grep -v "LowQual" | grep -c "<DEL>"
Combien de variants de bonne qualité sont de type Duplication ?
grep -v "^#" SV_calling_illumina.vcf | grep -v "LowQual" | grep -c "<DUP>"
Combien de variants de bonne qualité sont de type Inversion?
grep -v "^#" SV_calling_illumina.vcf | grep -v "LowQual" | grep -c "<INV>"
Combien de variants de bonne qualité sont de type inter chromosomique translocation ?
grep -v "^#" SV_calling_illumina.vcf | grep -v "LowQual" | grep -c "<BND>"
Combien de variants de bonne qualité sont de type Insertion ?
grep -v "^#" SV_calling_illumina.vcf | grep -v "LowQual" | grep -c "<INS>"
c) extraction des informations des délétions¶
Extraire les délétions de bonne qualité:
grep -v "^#" SV_calling_illumina.vcf | grep -v "LowQual" | grep "<DEL>"
Récupération du chromosome et de la position
#Récupération du start des variants
grep -v "^#" SV_calling_illumina.vcf | grep -v "LowQual" | grep "<DEL>" | \
cut -f1,2 > delly_del_start.txt
Récupération des autres informations
#Récupération des autres informations
grep -v "^#" SV_calling_illumina.vcf | grep -v "LowQual" | grep "<DEL>" | \
cut -f8 | cut -d ";" -f1,4,5,13 | sed "s/;/\t/g" > delly_del_info.txt
Fusion des 2 fichiers et formattage:
#Fusion des deux fichiers
paste -d '\t' delly_del_start.txt delly_del_info.txt > delly_del.txt
#Formattage et ménage
awk '{OFS="\t" ; print $1,$2,$4,$3,$5,$6}' delly_del.txt | sed "s/END=//g" > delly_del.tsv
Nettoyage:
rm delly_del_info.txt delly_del_start.txt delly_del.txt
2. detection de SV à partir de données MinION long-reads¶
Ici pour les données long reads, nous allons utiliser un autre caller: sniffles (papier ici https://www.nature.com/articles/s41592-018-0001-7, code là:https://github.com/fritzsedlazeck/Sniffles)
module load sniffles/1.0.11
sniffles --help
Pas besoin de génome de référence!
a) indexation du bam¶
samtools index mapping_minion_chr10_500kb.bam
b) Variant calling¶
sniffles -l 100 -m mapping_minion_chr10_500kb.bam -v SV_calling_minion.vcf
grep '^#' SV_calling_minion.vcf
c) comptage du nombre de SVs¶
Combien de SV sont positionées sur le chr10?
cat SV_calling_minion.vcf | grep ^chr_10 | wc -l
Combien sont des délétions?
cat SV_calling_minion.vcf | grep ^chr_10 | grep "DEL" | wc -l
Combien sont des duplications?
cat SV_calling_minion.vcf | grep ^chr_10 | grep "DUP" | wc -l
Combien d'entre elles sont des inversions?
cat SV_calling_minion.vcf | grep ^chr_10 | grep "INV" | wc -l
Combien d'entre elles sont des translocations?
cat SV_calling_minion.vcf | grep ^chr_10 | grep -w "TRA" | wc -l
Combien d'entre elles sont des insertions?
cat SV_calling_minion.vcf | grep ^chr_10 | grep "INS" | wc -l
d) extraction des positions des délétions¶
cat SV_calling_minion.vcf | grep ^chr_10 | grep "DEL" | head
Extraire le chromosome et la postion de début de la délétion
cat SV_calling_minion.vcf | grep ^chr_10 | grep "DEL" | cut -f -1,2 > sniffles_del_start.txt
Extraire la position de fin de la délétion
cat SV_calling_minion.vcf | grep ^chr_10 | grep "DEL" | cut -d ";" -f 4 | cut -d "=" -f 2 > sniffles_del_stop.txt
Extraire les informations complémentaires des délétions
cat SV_calling_minion.vcf | grep ^chr_10 | grep "DEL" | cut -f 8 | cut -d ";" -f 1 > sniffles_del_infos.txt
Coller tout ensemble
paste sniffles_del_start.txt sniffles_del_stop.txt sniffles_del_infos.txt > sniffles_del.tsv
rm sniffles_del_start.txt sniffles_del_stop.txt sniffles_del_infos.txt
3. Comparaison des résultats de Delly et Sniffles¶
Ici nous allons observer les résultats de Delly et de Sniffles
cat delly_del.tsv
cat sniffles_del.tsv
Les délétions sont elles aux mêmes positions?