Annexe 5 : Impacts de certaines actions - Dioufamad/SNPs_Calling GitHub Wiki

Prologue :

1 - Vision :

*L'idée de cette partie est de montrer comment certaines précautions prises dans le pipeline aurait eu un impact si elles n'étaient pas prise.

  • le point de comparaison est fait avec les sorties normalisées du caller et de faire vt peek.
  • on normalise pour enlever tous les appels identiques mais resprésentés juste différemment dans le vcf de sortie.
  • vt peek permet d'avoir un résumé assez fourni des différents types de variants (on se focalisera surtout sur observer le nombre de SNP, MNP, CLUMPED,et les taux de ts/tv)
  • une possibilité sera de lancer bk et pindel sur ces sorties ou les bam ayant donnés ces sorties pour aller pécher indels et sv
  • on pourrait étendre en soumettant ces sorties/bam le ayant donné aux sv callers de ma sélection
  • L'objectif est de poser ce que peut se passer pour chaque action et en même temps voir ce que nos softs puisent selon le traitement des données qu'on leur fourni

2 - Méthode de comparaison

Préparation des données

NB : pour manipuler des deux vcf souvent il est demandé la version compressée gz du vcf. pour le créer, faire : Download & install with make the latest tabix jar at http://sourceforge.net/projects/samtools/files/tabix/. see here for the installation help (http://iobio.io/2015/09/03/install-run-tabix/). Once you have successfully compiled Tabix (and bgzip), you can compress your vcf files with the following command: $ /path/to/tabix/bgzip myvcf.vcf $ /path/to/tabix/tabix -p vcf myvcf.vcf.gz These two commands will create two new files: the compressed vcf (.vcf.gz) and the index (.vcf.gz.tbi). Now you are ready to run vcf.iobio or gene.iobio, selecting these two files when prompted. The manual page for Tabix and BGZip can be found at http://samtools.sourceforge.net/tabix.shtml.

Comparaison

en le faisant avec la suite vt, pas besoin de .gz et les vcf sont obtenus avec leur index. on peut y aller directement.

  • creer le fichier de comparison :
 touch trial_compare.txt
  • normalisation du vcf (a faire pour chacun des vcf à comparer)
vt normalize raw.vcf -r ../../test8_0_REFERENCE_with_faidx/MSU7.fa -o trial1.normalized.vcf &>> trial_compare.txt
  • eliminer les duplicats : (a faire pour chacun des vcf à comparer)
vt uniq trial1.normalized.vcf -o trial1.normalized.uniq.vcf &>> trial_compare.txt
  • summary du contenu du vc : (à faire pour chacun des vcf à comparer)
vt peek trial1.normalized.uniq.vcf &>> trial_compare.txt

**Démarche : ** we have these comparisons to make

nom_comparaison = Compare_fixMBQS_vs_allowMBQSsinceIR_vs_allowMBQSsinceRTC_onr900q.txt
vcf1 = ../test8_3_1_HaplotypeCaller/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realignedbam_raw_snp_indels.vcf
vcf2 = ../test8_3_1_HaplotypeCaller/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realignedbam_misencBQSallowed_raw_snp_indels.vcf
vcf3 = ../test8_3_1_HaplotypeCaller/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realignedbam_aMBQSsinceRTC.bam_raw_snp_indels.vcf

nom_comparaison = Compare_ppnodup_vs_aranodup_onr900q.txt
vcf1 = ../test8_3_1_HaplotypeCaller/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realignedbam_aMBQSsinceRTC.bam_raw_snp_indels.vcf
vcf2 = ../test8_3_1_HaplotypeCaller/AX227_vsMSU7_bwamem_ara_nodup_sorted_idxed_realignedbam_aMBQSsinceRTC.bam_raw_snp_indels.vcf

nom_comparaison = Compare_HCwithIR_vs_HCwoIR_onr900q.txt
vcf1 = ../test8_3_1_HaplotypeCaller/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realignedbam_aMBQSsinceRTC.bam_raw_snp_indels.vcf
vcf2 = ../test8_3_1_HaplotypeCaller/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed.bam_noindelrealign_raw_snp_indels.vcf

nom_comparaison = Compare_cleanedfq_vs_rawfq_onr900q.txt
vcf1 = ../test8_3_1_HaplotypeCaller/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realignedbam_aMBQSsinceRTC.bam_raw_snp_indels.vcf
vcf2 = ../test8_3_1_HaplotypeCaller/r900_AX227_raw_vsMSU7_bwamem_pp_nodup_sorted_idxed_realigned_aMBQSsinceRTC_raw_snp_indels.vcf

nom_comparaison = Compare_HConr900q_vs_HConbionfoq.txt
vcf1 = ../test8_3_1_HaplotypeCaller/r900_AX227_raw_vsMSU7_bwamem_pp_nodup_sorted_idxed_realigned_aMBQSsinceRTC_raw_snp_indels.vcf
vcf2 = ../test8_3_1_HaplotypeCaller/bioinfq_AX227_raw_vsMSU7_bwamem_pp_nodup_sorted_idxed_realigned_aMBQSsinceRTC_raw_snp_indels.vcf

nom_comparaison = Compare_HC_vs_UG_onbioinfoq.txt
vcf1 = ../test8_3_1_HaplotypeCaller/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realignedbam_aMBQSsinceRTC.bam_raw_snp_indels.vcf
vcf2 = ../test8_3_1_UnifiedGenotyper/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realigned_aMBQSsinceRTC_byUG_raw_snp_indels.vcf

the unique set of vcf files used are these :

vcf1 = ../test8_3_1_HaplotypeCaller/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realignedbam_raw_snp_indels.vcf
vcf2 = ../test8_3_1_HaplotypeCaller/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realignedbam_misencBQSallowed_raw_snp_indels.vcf
vcf3 = ../test8_3_1_HaplotypeCaller/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realignedbam_aMBQSsinceRTC.bam_raw_snp_indels.vcf

nom_comparaison = Compare_ppnodup_vs_aranodup_onr900q.txt
vcf2 = ../test8_3_1_HaplotypeCaller/AX227_vsMSU7_bwamem_ara_nodup_sorted_idxed_realignedbam_aMBQSsinceRTC.bam_raw_snp_indels.vcf

nom_comparaison = Compare_HCwithIR_vs_HCwoIR_onr900q.txt
vcf2 = ../test8_3_1_HaplotypeCaller/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed.bam_noindelrealign_raw_snp_indels.vcf

nom_comparaison = Compare_cleanedfq_vs_rawfq_onr900q.txt
vcf2 = ../test8_3_1_HaplotypeCaller/r900_AX227_raw_vsMSU7_bwamem_pp_nodup_sorted_idxed_realigned_aMBQSsinceRTC_raw_snp_indels.vcf

nom_comparaison = Compare_HConr900q_vs_HConbionfoq.txt
vcf2 = ../test8_3_1_HaplotypeCaller/bioinfq_AX227_raw_vsMSU7_bwamem_pp_nodup_sorted_idxed_realigned_aMBQSsinceRTC_raw_snp_indels.vcf
nom_comparaison = Compare_HC_vs_UG_onbioinfoq.txt
vcf2 = ../test8_3_1_UnifiedGenotyper/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realigned_aMBQSsinceRTC_byUG_raw_snp_indels.vcf

we just have to treat them and the rest is reporting comparisons to a table of our making :

to prepare the comparison files only these three commandes are to be done for each vcf to compare (each time replace with the root name of the vcf to be analysed):

vt normalize AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realigned_aMBQSsinceRTC_byUG_raw_snp_indels.vcf -r ../test8_0_REFERENCE_with_faidx/MSU7.fa -o AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realigned_aMBQSsinceRTC_byUG_raw_snp_indels.normalized.vcf &> AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realigned_aMBQSsinceRTC_byUG_raw_snp_indels.normalize_stdout.txt

vt uniq AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realigned_aMBQSsinceRTC_byUG_raw_snp_indels.normalized.vcf -o AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realigned_aMBQSsinceRTC_byUG_raw_snp_indels.normalized.uniq.vcf &> AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realigned_aMBQSsinceRTC_byUG_raw_snp_indels.uniq_stdout.txt

vt peek AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realigned_aMBQSsinceRTC_byUG_raw_snp_indels.normalized.uniq.vcf &> AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realigned_aMBQSsinceRTC_byUG_raw_snp_indels.peek_stdout.txt

I-Problème 1 : -fix_misencoded_quality_scores ou -allowPotentiallyMisencodedQuals ou --useOriginalQualities :

  • l'option fix_misencoded_quality_scores permet de changer l'encodage en soustraiyant 31 et trvaillant avec l'ancien encodage Sanger. Si ceci n'est pas fait et que HC travaille avec l'ancien encodage et que c'est la nouvelle qui est présente, il y'aura bcp de callings ratés. voyons si ceci est le cas : comparer 2 calling non filtré (celui obtenu sans le fix et celui avec le fix). -Résultat : en voulant lancer la meme cmd de l'HC mais avec le fix des quality scores des bases, on a : ##### ERROR MESSAGE: Bad input: while fixing mis-encoded base qualities we encountered a read that was correctly encoded; we cannot handle such a mixture of reads so unfortunately the BAM must be fixed with some other tool
  • Ceci veut dire que les quality scores des bases sont deja reparées. ceci a déja été fait en effet avec l'indelrealigner en commençant à le faire depuis la construction du target intervals avec targetcreator. ce que l'on peut faire c'est faire un un targetcreator et un indelrealigner sans rectifier les quality scores des bases et faire un calling HC dessus et comparer. Nous avons les 3 options suivantes pour bien gérer les quality scores et fair travailler le caller sur le bon encodage : --fix_misencoded_quality_scores / -fixMisencodedQuals Fix mis-encoded base quality scores

--allow_potentially_misencoded_quality_scores / -allowPotentiallyMisencodedQuals Ignore warnings about base quality score encoding . --useOriginalQualities / -OQ Use the base quality scores from the OQ tag

  • Choix des options à comparer : Nous n'avons pas encore fait de recalibrage alors on ne testera pas l'impact de -OQ car les qualités qu'il utilisera seront juste celles qui sont là. Nous testerons donc les 2 options (-allowPotentiallyMisencodedQuals et -fixMisencodedQuals) pour génerer un bam non fixed, faire un calling sur celui-ci, obtenir un vcf, le normaliser avec vt et enfin comparer les vcf normalisés

I-1-allowing the misencodedbase quality scores since the realigner instead of since the realignertargetcreator

  • comparaison des deux options :
java -Xmx12g -jar /usr/local/gatk-3.3/GenomeAnalysisTK.jar -T IndelRealigner -R ../test8_0_REFERENCE_with_faidx/MSU7.fa -I ../test6_sortsam/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed.bam -targetIntervals ../test8_1_RealignerTargetCreator/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed.intervals -o ./AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realigned_misencBQSallowed.bam -allowPotentiallyMisencodedQuals
(1 bam + 1 bai)
puis faire le calling avec HC : 
java -Xmx12g -jar /usr/local/gatk-3.3/GenomeAnalysisTK.jar -T HaplotypeCaller -R ../test8_0_REFERENCE_with_faidx/MSU7.fa -I ../test8_2_IndelRealigner/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realigned_misencBQSallowed.bam -o ./AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realignedbam_misencBQSallowed_raw_snp_indels.vcf -rf BadCigar -allowPotentiallyMisencodedQuals

More SNPs detected (+10 000). less indels, complex variants unchanged.
Reason why : long ago, the gatk tools werent really able to use the phred64 quality format instead allowing them to substract 31 out of every quality value made the scale comes back to the one of phred33. but by doing so, lower values are considered by the tools and they off certain detections while they are still good just because of that lower value. a new option in gatk helps working with the actual phred64 and to "save" those detections.

Lesson : we will keep working with the allow option and the phred64 then. Steps to update : adopter la nouvelle commande pour l'indelrealigner faites pr generer le allowedBQS et faire un HC.

I-2-allowing the misencoded base quality score since realignertargetcreator instead of since the realigner

creation du target intervals avec un allowedMBQS :
java -Xmx12g -jar /usr/local/gatk-3.3/GenomeAnalysisTK.jar -T RealignerTargetCreator -R ../test8_0_REFERENCE_with_faidx/MSU7.fa -I ../test6_sortsam/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed.bam -o ./AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_aMBQSsinceRTC.intervals -allowPotentiallyMisencodedQuals
NB : aMBQSsinceRTC : missencodedbasequality scores are allowed since the realignertargetcreator (the first thing done was allowing them just since the realigner, we will compare afterwards with that one)
puis faire un indelrealigner avec allowedMBQS : 
java -Xmx12g -jar /usr/local/gatk-3.3/GenomeAnalysisTK.jar -T IndelRealigner -R ../test8_0_REFERENCE_with_faidx/MSU7.fa -I ../test6_sortsam/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed.bam -targetIntervals ../test8_1_RealignerTargetCreator/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_aMBQSsinceRTC.intervals -o ./AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realigned_aMBQSsinceRTC.bam -allowPotentiallyMisencodedQuals
(1 bam + 1 bai)
puis faire le calling avec HC : 
java -Xmx12g -jar /usr/local/gatk-3.3/GenomeAnalysisTK.jar -T HaplotypeCaller -R ../test8_0_REFERENCE_with_faidx/MSU7.fa -I ../test8_2_IndelRealigner/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realigned_aMBQSsinceRTC.bam -o ./AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realignedbam_aMBQSsinceRTC.bam_raw_snp_indels.vcf -rf BadCigar -allowPotentiallyMisencodedQuals

Your job 1373730 ("amad_hc_ara_aMBQSsinceRTC") has been submitted. output in test8_3_1_HaplotypeCaller_x_amad_hc_ara_aMBQSsinceRTC. move it to the normal folder

  • sdtout de HC avec allowedMBQS since RTC (sortie 11.1) (get it from the old_jobs folder)
  • pour bcp de softs utilisés ds la suite, il est nécessaire et meme preferable en termes de taille d'avoir un .gz et aussi indexer pour accélérer les analyses. on a le choix entre un bzip puis un tabix ou vt view -o .gz puis vt index .vcf. on choisis de faire un bzip puis un tabix pour rester comparable à ce que l'on a fait plus haut. on utilisera vt pour les autres analyses. $ /path/to/tabix/bgzip myvcf.vcf $ /path/to/tabix/tabix -p vcf myvcf.vcf.gz
  • on fait un vt peek pour avoir un aperçu de notre calling

II - Problème 2 : selection via les sam-flags

II-1- impact de ne pas restreindre les alignements à seulement les reads proprement alignés en pair

  • ne pas restrindre les alignements aux reads parfaiment alignés en pair (pas d'option -f 0x2)
  • et ne pas garder les duplicats toujours (option -F 0x400 incluse)
[adiouf@node21 test5_1_sview_ara_nodup]$ samtools view -b -S -h -F 0x400 -o AX227_vsMSU7_bwamem_ara_nodup.bam ../test4_2_bwamem/AX227_vsMSU7_bwamem.sam
[adiouf@node21 test5_1_sview_ara_nodup]$ ls
AX227_vsMSU7_bwamem_ara_nodup.bam

NB : ara : all reads aligned

  • sort + index en une fois
java -Xmx12g -jar /usr/local/picard-tools-2.5.0/picard.jar SortSam VALIDATION_STRINGENCY=SILENT CREATE_INDEX=TRUE SORT_ORDER=coordinate INPUT=../test5_1_sview_ara_nodup/AX227_vsMSU7_bwamem_ara_nodup.bam OUTPUT=./AX227_vsMSU7_bwamem_ara_nodup_sorted_idxed.bam
  • stats sur lalignement
samtools flagstat ../test6_sortsam/AX227_vsMSU7_bwamem_ara_nodup_sorted_idxed.bam > AX227_vsMSU7_bwamem_ara_nodup_sorted_idxed.samtoolsFlagstat

rapport généré :

6823393 + 0 in total (QC-passed reads + QC-failed reads)
12507 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
6657909 + 0 mapped (97.57% : N/A)
6810886 + 0 paired in sequencing
3405443 + 0 read1
3405443 + 0 read2
5547448 + 0 properly paired (81.45% : N/A)
6592010 + 0 with itself and mate mapped
53392 + 0 singletons (0.78% : N/A)
444824 + 0 with mate mapped to a different chr
142771 + 0 with mate mapped to a different chr (mapQ>=5)
  • prétraitement au calling : realignertargetcreator NB : c'est à partir d'ici que nous devons appliquer la lesson retenue d'utiliser les allowmBQS au lieu de fixmBQS
java -Xmx12g -jar /usr/local/gatk-3.3/GenomeAnalysisTK.jar -T RealignerTargetCreator -R ../test8_0_REFERENCE_with_faidx/MSU7.fa -I ../test6_sortsam/AX227_vsMSU7_bwamem_ara_nodup_sorted_idxed.bam -o ./AX227_vsMSU7_bwamem_ara_nodup_sorted_idxed_aMBQSsinceRTC.intervals -allowPotentiallyMisencodedQuals
  • prétraitement au calling : indelrealigner
java -Xmx12g -jar /usr/local/gatk-3.3/GenomeAnalysisTK.jar -T IndelRealigner -R ../test8_0_REFERENCE_with_faidx/MSU7.fa -I ../test6_sortsam/AX227_vsMSU7_bwamem_ara_nodup_sorted_idxed.bam -targetIntervals ../test8_1_RealignerTargetCreator/AX227_vsMSU7_bwamem_ara_nodup_sorted_idxed_aMBQSsinceRTC.intervals -o ./AX227_vsMSU7_bwamem_ara_nodup_sorted_idxed_realigned_aMBQSsinceRTC.bam -allowPotentiallyMisencodedQuals
  • calling : HC (done!)
java -Xmx12g -jar /usr/local/gatk-3.3/GenomeAnalysisTK.jar -T HaplotypeCaller -R ../test8_0_REFERENCE_with_faidx/MSU7.fa -I ../test8_2_IndelRealigner/AX227_vsMSU7_bwamem_ara_nodup_sorted_idxed_realigned_aMBQSsinceRTC.bam -o ./AX227_vsMSU7_bwamem_ara_nodup_sorted_idxed_realignedbam_aMBQSsinceRTC.bam_raw_snp_indels.vcf -rf BadCigar -allowPotentiallyMisencodedQuals

II-2- impact de garder les duplicats

  • ne pas restreindre les alignements aux reads parfaiment alignés en pair (pas d'option -f 0x2)
  • et garder les duplicats aussi (pas d'option -F 0x400)
[adiouf@node21 test5_sview]$ samtools view -b -S -h -o AX227_vsMSU7_bwamem_ara_withdup.bam ../test4_2_bwamem/AX227_vsMSU7_bwamem.sam
[adiouf@node21 test5_sview]$ ls
AX227_vsMSU7_bwamem_ara_withdup.bam

NB : ara : all reads aligned ; withdup : duplicats inclus

  • sort + index en une fois
java -Xmx12g -jar /usr/local/picard-tools-2.5.0/picard.jar SortSam VALIDATION_STRINGENCY=SILENT CREATE_INDEX=TRUE SORT_ORDER=coordinate INPUT=../test5_sview/AX227_vsMSU7_bwamem_ara_withdup.bam OUTPUT=./AX227_vsMSU7_bwamem_ara_withdup_sorted_idxed.bam

done in SortSam done. Elapsed time: 3,04 minutes.

  • stats sur lalignement
samtools flagstat ../test6_sortsam/AX227_vsMSU7_bwamem_ara_withdup_sorted_idxed.bam > AX227_vsMSU7_bwamem_ara_withdup_sorted_idxed.samtoolsFlagstat

rapport généré :

6823393 + 0 in total (QC-passed reads + QC-failed reads)
12507 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
6657909 + 0 mapped (97.57% : N/A)
6810886 + 0 paired in sequencing
3405443 + 0 read1
3405443 + 0 read2
5547448 + 0 properly paired (81.45% : N/A)
6592010 + 0 with itself and mate mapped
53392 + 0 singletons (0.78% : N/A)
444824 + 0 with mate mapped to a different chr
142771 + 0 with mate mapped to a different chr (mapQ>=5)
  • ce bam n'a pas de duplicats (on le voit en comparant ce flagstat avec celui precedent de ts les alignements sans les duplicats). conitnuer la comparaison n'aurait pas de sens car ce sont les memes contenus. Solution : construire un script qui va utiliser tous les fastq et faire les étapes jusqu'au flagstat puis on recupère les flagstats. Cependant cela a déjà été fait avec nos anciens travaux de generations de bam pour la selection d'outils de detection de SV. En examinant ces flagstat (on a verifié qu'aucune étape éliminant les duplicats n'a été faite avant les avoir obtenus), on voit que nous n'avons aucun cas de bam avec des duplicats.
  • To do : nous devons allons chercher des donnes avec des duplicats (2 fastq, et une ref) puis run le pipeline sans et avec les duplicats éliminés via les flags.

III - Problème 3 : impact de ne pas faire un indelrealigner

il s'agit de faire le HC directement avec le bam qui a été sorted et indexed avec picardtools sortsam.

  • calling : HC (done!)
java -Xmx12g -jar /usr/local/gatk-3.3/GenomeAnalysisTK.jar -T HaplotypeCaller -R ../test8_0_REFERENCE_with_faidx/MSU7.fa -I ../test6_sortsam/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed.bam -o ./AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed.bam_noindelrealign_raw_snp_indels.vcf -rf BadCigar -allowPotentiallyMisencodedQuals

IV - Problème 4 : impact de ne pas nettoyer les reads avant alignement

*4_2_bwa_mem :

bwa mem -t 2 -T 30 -M ../test4_1_bwaindex/MSU7.fa ../FASTQ-RAW-2MOD/AX227_1.fq.gz ../FASTQ-RAW-2MOD/AX227_2.fq.gz -R '@RG\tID:AX227\tSM:AX227\tPL:Illumina' > AX227_raw_vsMSU7_bwamem.sam
  • runtime : Real time: 1876.650 sec; CPU: 3777.946 sec *samtoolsview selection by tag + conversion in bam :
samtools view -b -S -h -f 0x2 -F 0x400 -o AX227_raw_vsMSU7_bwamem_pp_nodup.bam ../test4_2_bwamem/AX227_raw_vsMSU7_bwamem.sam

*picardtools sorting and bam indexing :

java -Xmx12g -jar /usr/local/picard-tools-2.5.0/picard.jar SortSam VALIDATION_STRINGENCY=SILENT CREATE_INDEX=TRUE SORT_ORDER=coordinate INPUT=../test5_sview/AX227_raw_vsMSU7_bwamem_pp_nodup.bam OUTPUT=./AX227_raw_vsMSU7_bwamem_pp_nodup_sorted_idxed.bam

-Elapsed time: 2,56 minutes.

  • samtools flagstats :
flagstat ../test6_sortsam/AX227_raw_vsMSU7_bwamem_pp_nodup_sorted_idxed.bam > AX227_raw_vsMSU7_bwamem_pp_nodup_sorted_idxed.samtoolsFlagstat

stats obtained:

5534363 + 0 in total (QC-passed reads + QC-failed reads)
5145 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
5534363 + 0 mapped (100.00% : N/A)
5529218 + 0 paired in sequencing
2764609 + 0 read1
2764609 + 0 read2
5529218 + 0 properly paired (100.00% : N/A)
5529218 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
  • realignertargetcreator :
java -Xmx12g -jar /usr/local/gatk-3.3/GenomeAnalysisTK.jar -T RealignerTargetCreator -R ../test8_0_REFERENCE_with_faidx/MSU7.fa -I ../test6_sortsam/AX227_raw_vsMSU7_bwamem_pp_nodup_sorted_idxed.bam -o ./AX227_raw_vsMSU7_bwamem_pp_nodup_sorted_idxed.intervals -allowPotentiallyMisencodedQuals
  • runtime : Total runtime 356.06 secs, 5.93 min, 0.10 hours
  • indelrealigner :
java -Xmx12g -jar /usr/local/gatk-3.3/GenomeAnalysisTK.jar -T IndelRealigner -R ../test8_0_REFERENCE_with_faidx/MSU7.fa -I ../test6_sortsam/AX227_raw_vsMSU7_bwamem_pp_nodup_sorted_idxed.bam -targetIntervals ../test8_1_RealignerTargetCreator/AX227_raw_vsMSU7_bwamem_pp_nodup_sorted_idxed.intervals -o ./AX227_raw_vsMSU7_bwamem_pp_nodup_sorted_idxed_realigned_aMBQSsinceRTC.bam -allowPotentiallyMisencodedQuals
  • runtime : Total runtime 170.82 secs, 2.85 min, 0.05 hours

(Your job 1377981 ("amad_rawfq_hc_pp_nodup") has been submitted)

  • runtime : Your job 1377982 ("amad_rawfq_hc_pp_nodup_bioinfq") has been submitted)
  • runtime :

Problème 5 : impact de faire un UnifiedGenotyper au lieu d'un HC

l'expérience des developpers de Toggle est que qu'en 7 jours HC a tourné et la sélection des variants n'était meme pas encore faite. Et qu'à coté Unifiedgenotyper donnait en 3 jours les résultats jusqu'à finir la sélection des variants. En comparant aussi, les résultats n'étaient pas si différents en terme de gain de variants détectés avec HC au point de devoir y mettre 4 jours+ encore. d'où UG a été choisi pour ces analyses. NB : ne pas oublier que HC fait une detection en suivant un cadre qu'il deplace et donc recorrige si certaines incohérences lui apparaissent sur un cadre ou alors cela donne naissance à de nouveau variants. Cela est utile en medecine humaine ou le moindre variant à fréquence aussi faible qu'il puisse etre peut une cible thérapeutique et donc dans ce cas HC ou des callers plus orientées low freaquencies variants vont etre demandés et les faibles différences entre sorties ne seront pas à occulter Explication correcte de la différence entre HC et UG: HC et UG utilisent tous deux un model de vraisemblance bayesien pour le genotypage. UG va faire ce calcul et émettre en meme temps les genotypes et frequences alléliques les plus vraissemblables et cela pour la population de N echantillons fournis en entrée, un genotype etant pour un echantillon. HC lui, avant de faire ce calcul va, dans la fenetre sur laquelle ce calcul se fait, , elire des active regions (regions où les preuves de variations sont significatives), construire des graphes DE Bruijn sur ces regions pour faire des assemblages. chacun de ces assemblages est un haplotype qu'il va aligner par l'algo Smith-Waterman à l'haplotype que represente la reference. Pour determiner la vraisemblance de chaque haplotype, chacun est aligné avec un read de manière pairwise par lalgo PairHMM : le resultat est une matrice dans laquelle des valeurs de vraisemblances sont affectées à des haplotypes selon selon les données liées aux reads (profondeur etc). Les genotypes pour chaque site potentiellement variant sont calculés en appliquant Bayes : utiliser les vraisemblances des alleles sachant que l'on dispose des données entières des reads pour calculer le genotype pour un échnatillon sachant les données de read liées à cet echantillon.

Pour voir l'impact en termes de pourcentage de variants en plus detectés par HC par rapport à UG, il faudrait comparer 2 detections HC et UG. Leurs options utilisées doivent être assez proche en terme de vérité biologique pour le calcul q'elles font. il faudrait donc un set de valeurs d'options de UG assez proche de celles utilisiées jusqu'ici par HC