I - Données

• folder FAST-RAW-2MOD : contains the 20 fastq files for the 10 samples of kasalath rice
• folder FAST-RAW-MINI : contains for 3 samples, the reduced versions of their fasq files to 10 000 reads each (20000 for a samples) and we have 3 samples of those. destined to smal and fast tests.
• folder test1_fastqc : contains the results of the FASTQC tool for quality control

II - Annalyse et correction eventuelle

1-Analyse qualité avant correction

fastqc FASTQ-RAW-2MOD/AX227_1.fq.gz -o test1_fastqc/
fastqc FASTQ-RAW-2MOD/AX227_2.fq.gz -o test1_fastqc/

result : from the sample AX227 we have 3405443 reads pairs of 83b.

• quality score across bases is above 28 for most reads. we can strip out the reads the bases with less than 30 quality
• per base sequence sequence content is quite shaky up until the 9th base: we can trim that portion 1st-9th base. their per base gc content is also shaky maybe their trimmingwill resolve it.

2-Correction : trimming and low quality removing

on pense que les 9 1ères bases sont des adaptateurs. or l'otption -u/-U enlève des bases memes si leurs qualités est bonne mais qu'elle sont juste en extrmité sur la distance de trimming désignée. nous n'allons alors pas utilser d'abord cette option. essayons d'nlever les adaptateurs puis voir ce que cela donne. si il y'an encore des positions en extremité qui sont mlgré toujours "shaky" en terme de per base content, on utilisera l'option -u/-U. model de la ligne de commande :

cutadapt -O 10 -b adapt1 -B adapt1 -m 35 -q 30,30 cutadapt_output_folder/R1_output.fq -p cutadapt_output_folder/R2_output.fq R1.fq R2.fq

notre vraie commande :

cutadapt -b ACACTCTTTCCCTACACGACGCTCTTCCGATCT -B ACACTCTTTCCCTACACGACGCTCTTCCGATCT -O 10 -m 35 -b TCTAGCCTTCTCGCAGCACATCCCTTTCTCACA -B TCTAGCCTTCTCGCAGCACATCCCTTTCTCACA -q 30,30 -o ./test2_cutadapt/AX227_1.CUTADAPT.fastq -p ./test2_cutadapt/AX227_2.CUTADAPT.fastq ./FASTQ-RAW-2MOD/AX227_1.fq.gz ./FASTQ-RAW-2MOD/AX227_2.fq.gz

sdtout de cutadapt (sortie 1)

3-Analyse qualité après correction

faisons une analyse qualité pour voir si la qualité est acceptable maintenant

fastqc ./test2_cutadapt/*.fq -o test3_fastqc

Resultat : les 9 bases en 5' ayant un per base content "shaky" sont toujours là. il faut donc ajouter les options -u/ -U avec 9 NB : transfert regulier des données vers local soit pour visualisation ou backup

en se plaçant dans le folder adiouf_div dans le node19: 
rsync -arv ../adiouf_div bioinfo-nas.ird.fr:/home/adiouf/
puis en étant en local dans le miroir de notre dossier /home/adiouf, faire : 
rsync -arv --exclude=FASTQ-RAW-2MOD --exclude=FASTQ-RAW-MINI [email protected]:/home/adiouf/adiouf_div .

**Nouveau cycle de correction : **

2bis-réessayons le cutadapt (placés ds l'output folder):

cutadapt -u 9 -U 9 -b ACACTCTTTCCCTACACGACGCTCTTCCGATCT -B ACACTCTTTCCCTACACGACGCTCTTCCGATCT -O 10 -m 35 -b TCTAGCCTTCTCGCAGCACATCCCTTTCTCACA -B TCTAGCCTTCTCGCAGCACATCCCTTTCTCACA -q 30,30 -o ./AX227_1.cleaned.fastq -p ./AX227_2.cleaned.fastq ../FASTQ-RAW-2MOD/AX227_1.fq.gz ../FASTQ-RAW-2MOD/AX227_2.fq.gz

sdtout de cutadapt (sortie 2)

3bis- analyse qualité de controle

(placés ds l'output folder)
fastqc ../test2_cutadapt/AX227_1.cleaned.fastq -o ./
fastqc ../test2_cutadapt/AX227_2.cleaned.fastq -o ./

• qui montre : des reads entre 59/60 et 74 b. juste deux warnings (pr sequence gc content et sequence length distribution) et per base sequence content est stable maintenant.
• pour le psGCc, c'est normal car il est assez etonnant d'avoir pile ce qu'on a en théorie.
• et pour la sld, c'est acceptable car la taille des sequences est ici donnée pour nous informer mais on doit de toute façon utiliser cs sequences quel que soit leur taille.

III - Alignement

4-1 Etape 1/2 : indexer la reference

copier la reference (MSU7.fa ici) dans un dossier et y faire un :
bwa index MSU7.fa

sdtout de bwa index (sortie 3)

5 fichiers ont été générés et constituent l'index utilisé par bwa mem pour l'alignement

4-2 Etape 2/2 : lancer bwa mem

bwa mem -t 2 -T 30 -M ../test4_1_bwaindex/MSU7.fa ../test2_cutadapt/AX227_1.cleaned.fastq ../test2_cutadapt/AX227_2.cleaned.fastq -R '@RG\tID:AX227\tSM:AX227\tPL:Illumina' > AX227_vsMSU7_bwamem.sam

sdtout de bwa mem (sortie 4) et et on a comme resultat un fichier sam

A revoir :

• le nombre de threads utilisés avec -t 2 (est ce pertinent si l'on a reservé qu'un seul coeur sur le noeud
• pour travailler sur les memes données pour des SV, refaire le bwamem mais sans l'option -M sinon les régions secondaires des splitreads seront justes des alignements secondaires et non supllémentaires
• comprendre le stdout de la commande bwamem dans les détails et ce qu'il fait succesivement

IV - Post-traitement de l'alignement

5 - restriction des alignement aux alignement necessaires

on veut faire du snp calling alors les paires de reads proprement alignés sont ce qui est nécessaire. on ne garde que cela et on enlève les duplicats en même temps. on en profite pour avoir deja un bam

samtools view -b -S -h -f 0x2 -F 0x400 -o AX227_vsMSU7_bwamem_pp_nodup.bam ../test4_2_bwamem/AX227_vsMSU7_bwamem.sam

6 - ordonner et avoir un index en une fois

java jvm-args -jar picard.jar SortSam VALIDATION_STRINGENCY=SILENT CREATE_INDEX=TRUE SORT_ORDER=coordinate INPUT=output_folder/samtools_output_file.bam OUTPUT=output_folder/ptss_output_file.ba

ne marche pas! commande corrigée :

/usr/bin/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_pp_nodup.bam OUTPUT=./AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed.bam

commande re-corrigée (handler java raccourci) :

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_pp_nodup.bam OUTPUT=./AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed.bam

• NB : -Xmx12g donne la taille maximale de memoire disponible pour la JVM. un fichier .bam est ordonné et un index .bai sont produits et voici le stdout : sdtout de picardtools sortsam (sortie 5)

• revoir la liste des outils picard tools avec java -jar /usr/local/picard-tools-2.5.0/picard.jar -h

• utiliser la commande de toggle pour picardtools sortsam mais reoir après dans le détail ce que c'est

7 - stats sur l'alignement

samtools flagstat ../test6_sortsam/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed.bam > AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed.samtoolsFlagstat

le rapport généré montre :

5550596 + 0 in total (QC-passed reads + QC-failed reads)
3148 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
5550596 + 0 mapped (100.00% : N/A)
5547448 + 0 paired in sequencing
2773724 + 0 read1
2773724 + 0 read2
5547448 + 0 properly paired (100.00% : N/A)
5547448 + 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)

V - Pré-traitement au calling

• cette partie est faite avec GATK. La version 3.3 est utilisée dans toggle. Si l'on veut tout repeter avec toggle en s'assurant d juste le refaire plus rapidement de manière automatisée sans rien changer dans les output, il faudrait tester avec cette version.
• cependant à la date actuelle du 26/06/2018, la dernière version de GA## K est la 4.0. il faudrait plus tard la comparer avec la 3.3 de toggle et savoir sur laquelle rester pour la suite. Dans toute la suite, GATK 3.3 est utlisé.

8-1-RealignerTargetCreator

cherchons les intervalles nécessitant un realignement

- voici la commande utilisée pour picardtools un outil utilisant java comme RealignerTargetCreator. sinspirer de celle-ci pour lancer celle de gatk RTC (dans toute la suite on utilisera cette méthode pour écrire vite les nouvelles commandes c'est à dire copier l'ancienne commande, s'en inspirer pour créer la nouvelle en la modifiant puis effacer plus tard l'ancienne commande.) :
/usr/bin/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_pp_nodup.bam OUTPUT=./AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed.bam
- on deduit la commande suivante pour gatk-RTC:
/usr/bin/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.intervals

• NB1 : le message important dans les erreurs avec gatk est en fin de la stderror.
• NB2 : Pour vérifier les options d'un outil, faire : /usr/bin/java -Xmx12g -jar /usr/local/gatk-3.3/GenomeAnalysisTK.jar -T RealignerTargetCreator --help
• NB3 : le repertoire test8_0_REFERENCE_with_faidx contient 3 fichiers nécessaire (la reference .fa, son index .fai et un .dict et ils doivent avoir tous les 3 le meme prefixe tel que [abc.fa abc.fa.fa abc.dict]). Explication ici (https://gatkforums.broadinstitute.org/gatk/discussion/1601/how-can-i-prepare-a-fasta-file-to-use-as-reference). En résumé, faire ces deux commandes suivantes pour creer 1-un dictionnaire (contenant liste des contigs) et 2-un index (permettant d'accéder à une position précise pour une recherche plus rapide) :

java -jar /usr/local/picard-tools-2.5.0/picard.jar CreateSequenceDictionary REFERENCE=MSU7.fa OUTPUT=MSU7.dict

samtools faidx MSU7.fa 

• NB4 : cette erreur peut être rencontrée : ERROR MESSAGE: SAM/BAM file SAMFileReader [your_bam_file] appears to be using the wrong encoding for quality scores: we encountered an extremely high quality score. il s'agit alors du fait que la qualité des bases sont encodés dans l'ancien encodage qui est déprécié et que certains outils prennent mais que cela risque d'offset leurs calculs alors ils le signalent. L format de qualité reconnu par la specification SAM est Q0 == ASCII 33, c'est il y'a un ajout de 31 sur ce que les valeurs doivent actuellement etre. Les outils doivent soustraire alors juste 31 de chaque qualité de base et calculer avec cela. Pour cela il faut ajouter le flag --fix_misencoded_quality_scores (-fixMisencodedQuals) geeralement. l'exception est quand on doit ajouter -allowPotentiallyMisencodedQualspour l'etape "PrintReads" de base recalibration (BaseRecalibrator sera fait avec "--fix_misencoded_quality_scores (-fixMisencodedQuals)") les liens suivants sont interessants pour la question : (https://gatkforums.broadinstitute.org/gatk/discussion/1991/version-highlights-for-gatk-version-2-3), (https://gatkforums.broadinstitute.org/gatk/discussion/2500/depthofcoverage-producing-extremely-high-quality-score-error), (https://gatkforums.broadinstitute.org/gatk/discussion/2335/realignment-with-high-base-qualities)
• Pour les score de qualité des bases, il existe 4 encodages : --If the quality scores contain characters in the range ASCII 33 - 58 -> can only be Sanger --If FastQ file is known to be from an Illumina/Solexa platform AND the quality scores contain characters in the range ASCII 59 - 63 -> can only be Solexa/Illumina 1.0 (can be very good score of sanger ) --If ASCII characters 64 or 65 are used in quality scores -> cannot be Illumina 1.5+ (see whole discussin on the subject here https://www.biostars.org/p/911/)
• cette image sur (https://wiki.bits.vib.be/index.php/Identify_the_Phred_scale_of_quality_scores_used_in_fastQ) resume ls diffrences entr encodages
• Trois options sont à connaitre pour faire travailler le caller avec le bon encodage de base quality : --fix_misencoded_quality_scores / -fixMisencodedQuals Fix mis-encoded base quality scores By default the GATK assumes that base quality scores start at Q0 == ASCII 33 according to the SAM specification. However, encoding in some datasets (especially older Illumina ones) starts at Q64. This argument will fix the encodings on the fly (as the data is read in) by subtracting 31 from every quality score. Note that this argument should NEVER be used by default; you should only use it when you have confirmed that the quality scores in your data are not in the correct encoding. --allow_potentially_misencoded_quality_scores / -allowPotentiallyMisencodedQuals Ignore warnings about base quality score encoding This flag tells GATK to ignore warnings when encountering base qualities that are too high and that seemingly indicate a problem with the base quality encoding of the BAM or CRAM file. You should only use this if you really know what you are doing; otherwise you could seriously mess up your data and ruin your analysis. --useOriginalQualities / -OQ Use the base quality scores from the OQ tag This flag tells GATK to use the original base qualities (that were in the data before BQSR/recalibration) which are stored in the OQ tag, if they are present, rather than use the post-recalibration quality scores. If no OQ tag is present for a read, the standard qual score will be used.

la commande finalement fonctionnant est donc (handler java raccourci aussi):

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.intervals -fixMisencodedQuals

le stdout (sdtout de RealignerTargetCreator (sortie 6)) montre ce qui suit et un fichier .intervals (chaque ligne est ChrX:posA/posA-posB) est produit:

8-2-IndealRealigner

Realignons maintenant

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.bam -fixMisencodedQuals

NB : on reconnait que si cette commande fonctionne tout vcf dans la sortie ne contient plus une information sur des indels valable et que les infos sur les indels seront à éliminer

• la commande a été allégée (/usr/bin/java devient java)
• se rappeller que l'index du bam est important ici car il est parcouru
• La sortie est un bam contenant les alignments rectifiés et son index .bai
• le stdout est :(sdtout de IndelRealigner (sortie 7))

Une certaine idée serait de générer un flagstat ici dans une partie 8-2-1-stats sur bam realigné. on passde dessus pour le moment car quoi que puisse etre le resultat on va travailler avec le contenu du bam realigné et ce serait juste pour voir ce qui a changé globalement dans les nombres.

8-3-BaseRecalibrator

Reestimer les scores de qualité d'alignement des bases on a :

java -Xmx12g -jar /usr/local/gatk-3.3/GenomeAnalysisTK.jar -T BaseRecalibrator -R ../test8_0_REFERENCE_with_faidx/MSU7.fa -I ../test8_2_IndelRealigner/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realigned.bam -o ./AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realignedbam_recaldata.table -fixMisencodedQuals

on a l'erreur suivante : 
##### ERROR MESSAGE: Invalid command line: This calculation is critically dependent on being able to skip over known variant sites. Please provide a VCF file containing known sites of genetic variation.

• Le problème est ici qu'il nous faut produire un vcf avec quelques variants pour que la commande fonctionne.

When working with a non-model organism for which little or no known variant information is available the GATK developers recommend that you "bootstrap" your own list of known variants following this procedure:

• from (http://seqanswers.com/forums/showthread.php?t=44765), we found a solution :

Quote:
I'm working on a genome that doesn't really have a good SNP database yet. I'm wondering if it still makes sense to run base quality score recalibration without known SNPs.

The base quality score recalibrator treats every reference mismatch as indicative of machine error. True polymorphisms are legitimate mismatches to the reference and shouldn't be counted against the quality of a base. We use a database of known polymorphisms to skip over most polymorphic sites. Unfortunately without this information the data becomes almost completely unusable since the quality of the bases will be inferred to be much much lower than it actually is as a result of the reference-mismatching SNP sites.

However, all is not lost if you are willing to experiment a bit. You can bootstrap a database of known SNPs. Here's how it works:

• First do an initial round of SNP calling on your original, unrecalibrated data.

• Then take the SNPs that you have the highest confidence in and use that set as the database of known SNPs by feeding it as a VCF file to the base quality score recalibrator.

• Finally, do a real round of SNP calling with the recalibrated data. These steps could be repeated several times until convergence.
lancer un 1er calling et utiliser son vcf resultat est une solution. il faudra alors lancer ce 1er calling directement sans recalibrer la qualité de base calling. Dans Toggle ceci est fait (car cette recalibration est sautée parce que celà n'apporte pas beaucoup vu le temps que cla prnd) (faire une comparaison avec et sans recalibration après pour voir).

Pour le 1er vcf, on propose : Exemple de fichier de configuration utilisé par Toggle pour le snp calling sur du paired end data (external help 1)

cete dettection est assez proche de la notre en donnant des indices pour finir restant assez concis sur les filtres et selection de variants. on peut donc utiliser la detection qu'ils ont faites avec des options pas trop éloignées pour produire notre 1er vcf.

8-3-1 Bootstrap d'un vcf de variants connus : 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.bam -o ./AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realignedbam_raw_snp_indels.vcf -rf BadCigar 

• on ajoute le -rf BadCigar car d'experience sans celui-ci cela bloque. Ceci serait du au fait que sans celui-ci les indels consecutifs sont maintenus dans les données et HC ne les supporte pas. Rappel sur BadCigar : This read filter will filter out the following cases:

different length and cigar length
Hard/Soft clips in the middle of the cigar
starting with deletions (with or without preceding clips)
ending in deletions (with or without follow-up clips)
fully hard or soft clipped
consecutive indels in the cigar (II, DD, ID or DI)

• L'output de HC est un vcf non filtré et un .idx
• voici la sortie stdout : sdtout de HaplotypeCaller avec -fixMisencodedQuals(sortie 8)

8-3-2 Bootstrap d'un vcf de variants connus : variant filtration

java -Xmx12g -jar /usr/local/gatk-3.3/GenomeAnalysisTK.jar -T VariantFiltration -R ../test8_0_REFERENCE_with_faidx/MSU7.fa -V ../test8_3_1_HaplotypeCaller/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realignedbam_raw_snp_indels.vcf -o ./AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realignedbam_raw_snp_indels_filtered.vcf -filterName 'LowDepth' -filter 'DP<10' -filterName 'LowQual' -filter 'QUAL<30' 

• sdtout de VariantFiltration (sortie 9)

8-3-3 Bootstrap d'un vcf de variants connus : Selection du type de variants

java -Xmx12g -jar /usr/local/gatk-3.3/GenomeAnalysisTK.jar -T SelectVariants -R ../test8_0_REFERENCE_with_faidx/MSU7.fa -V ../test8_3_2_VariantFiltration/AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realignedbam_raw_snp_indels_filtered.vcf -o ./AX227_vsMSU7_bwamem_pp_nodup_sorted_idxed_realignedbam_filtered_snp_only.vcf -selectType SNP 

-sdtout de SelectVariants (sortie 10)

• Problème 1 : Pour être plus stringeant et entrainer un bon model pour le recalibrage, il faut lui fournir un vcf qui ai pas trop de snp pour éviter de laisser trop peu de positions à recalibrer mais ceux qui sont dans ce vcf devront être lourdement supportées en termes de preuves. on va jouer sur le DP (nbre d reads suportant un variant). il faudrait étudier le nbre de reads globalement sur notre echantillon.

K : mettre ici les options de filtrages à utiliser pour un vcf super stringeant

8-3-4 Comparaison de plusieurs points critiques du pipeline

voir ici : - Impacts de certaines actions (Annexe 5) certaines actions dans le pipeline sont comprises an termes de "pourqoui on le fait?"; il est interessant de connaitre au moins une fois ce que cela fait qu'elles n'aient pas été faites ou si elles sont mal faites. Les résultats de telles comparaisons sont figurées dans ce tableau avec des commentaires liées (https://docs.google.com/spreadsheets/d/1QAOB-1CDk5rP7dtROi1Wl6qRok7b6LbTmoeKB-ZkFzI/edit?usp=sharing). Cette figure propose un version simplifiée du pipeline avec les choix faits.

8-4-PrintReads

renvoyons les alignements corrigés

K

prochaine commande :

K