Demande : un protocole de SNP calling pour 15 individus

Definition : Les SNPs sont des variations ponctuelles au niveau génomique menant à un variant allélique relativement important dans la population (comparé au trait de référence). La notion de fréquence est ici importante et requise pour qu'une variation soit considérée comme SNP.

Besoin : Dans le cadre d'analyses de diversité, nous pouvons nous poser des questions sur l'origine des variations au niveau phénotypique (phénotype = ensemble des traits observables chez un individu). Or,le phénotype peut être rapporté au génotype (la séquence génétique ayant codé pour ce caractère) chez un individu. Alors, dans une population, la diversité entre sous-groupes pour un caractère considéré peut etre mieux comprise en :

• reçensant pour une séquence d'intéret les sites ayant des variations au sein de la population
• énumérant,pour chaque site, les variations qui sont significatifs comparé au % du site majoritairemnt retrouvé jusque-là (SNP)
• estimant les genotypes présents dans cette population (genotypage)
• mettant en relation certains phénotypes à des genotypes

Différence entre SNP calling et genotype calling:

• La présence d'un SNP est conclut par la detection d'une variation d'un site par rapport à la référence pour ce même site. Pour un seul individu, une variation montre cela que l'individu soit homozygote ou heterozygote pour ce site; pour plusieurs individus, un SNP est montré chaque fois qu'un individu, homozygote ou heterozygote, porte une base différente de la référence à ce site.
• Le génotype : est l'estimation des bases modifiées ou non modifiées pour un site en considérant tous les allèles possédés par l'individu à cette localisation.

Schéma de travail : Nous voulons mettre sur pied un protocole de SNP Calling. Nous allons ici donner des outils pour les différentes étapes. Nous expliquerons, à chaque fois, pourquoi nous avons choisi cet outil. Nous choisissons surtout de nous focaliser sur ce qui est fait depuis les fichiers issus du séquençage à l'obtention de variants et genotypes

Obstacle :

Depending on the type of DNA sequence, genomic or transcript, and the availability of a reference genome this can be more or less challenging. The two main challenges in detecting reliable SNPs in next generation sequencing (NGS) data are to: Distinguish sequencing artefacts from true genetic variation. The very basic is just by doing RFLP but again, I guess you need to do it many times with different types of restriction enzyme. If the fragment is pure and small enough just sequence it. and with sequencing comes the story of the sequencing errors. They may be coming from the sample itself, the environnement that surronded the logistics of the sequencing and also the platform used. Whatever it is, it can only be guarranteed an error rate but not the certainty of a base.

==========================================================================

Ce sont :

• des reads issus du séquençage fait avec une technologie de 2nd génération. Ils sont le plus souvent dans des fichiers .fastq.
• la séquence de référence (génome de référence ou assemblage) : fichier auformat .fasta

Nous choissionns de nous concentrer sur les données issues de technologies de séquençage dites de 2nd génération. Leurs avantages sont :

• un débit 100 fois plus rapide que celui proposé par la Sanger sequencing
• plus rapide : sur une période de quelques, le genome de toute espèce peut être disponible
• sequencing de novo à moindre cout
• cout moindre a pour conséquence, reséquençage d'espèces connues ou d'échantillon représentation
• permet de detecter des variaions genomiqes : SNP, Indels et SV. combiné aux autres analyses omiques dans lesquelles les NGS sont utilisables (ex: metabolomique, transcriptomique etc), on a une image presque complète des états du génome

NB : technologie à haut débit => analyse sensible à des erreurs de la technologie Conséquence : utiliser les donées de 2nde génération NGS pour de bons résultats dépend aussi de faire une sélection d'outils bioinformatiques fiables et adaptés.

• Un échantillonnage représentatif de la population (pour ne pas surestimer ni sous estimer la fréquence d'un variant) est fait.
• Extraction de l'ADN > Fragmentation de l'ADN > Creation de Librairie > Séquençage > fichiers .fastq

• L'échantillonnage va conditionner les données que l'on aura. Il faut s'assurer de bien prendre des individus représentatifs de la population et bien les marquer en conservant dans un fichier toutes les informations liées à leur existence, leur manifestation et leurs phénotypes.

• Le séquençage est fait en externe souvent. Cependant il est important de connaitre les étapes du séquençage pour comprendre le fond des choses. Nous choisissons de sauter cela pour le moment.

• La profondeur et la couverture : **La couverture : ** Elle indique l'étendue du génome sur laquelle chacune des bases qui s'y trouve a été séquencée un certain nombre de fois. Elle est donnée en %, pourcentage de bases sur cette ce génome ciblé qui ont été séquencée. La profondeur de couverture, attendue avant de faire le séquençage ou obtenue après qu'un séquençage soit fait, est le nombre moyen de fois qu'une base a a été séquencée, étant donné :

• un certain nombre de reads (N)
• une certaine longueur de chacun de ces reads (L)
• et en ce plaçant dans le cas d'une distribution aléatoire d'un génome dont la longueur de son haploide est G. NB : G est la taille de la partie du génome qui est ciblée et non la taille de tout le génome connue à ce jour. **Calculs : ** Et alors soit P cette profondeur de couverture, P = (L.N)/G. on a ici 4 paramètres l'un peut être déduit si les 3 autres sont disponibles en utilisant cette equation. Avantage d'une profondeur plus grande : Un séquençage de chaque base une seule fois dans une étude de diversité peut entrainer des faux positifs. En effet, une erreur de sequencage ne peut être dfférenciée d'une variation. Pour s'éloigner de cette ambiguité, plus on a de reads supportant une certaine à une position plus on est certain qu'il ne s'agit pas d'une erreur de sequençage. **NB : Pourquoi faire du paired-ends dans ce contexte d'éviter les inconvenances du séquencage? : ** la profondeur ne pourra lever toutes les ambiguités. Les gaps dans les assemblages dus aux séquences repetées de taille approchant ou même plus grande que celle des reads. Une paire de reads dans ce contexte permet de placer de manière non ambigue une séquence réppétée dont la taille n'excede pas la taille de l'insert (séquence dont les deux bouts ont été séquencés pour chacun 2 reads, chacun partant d'un bout).

**Lorsque ces reads sont alignés à une référence, nous devons faire resortir les sites où ils contiennent une variation d'une base à un site par rapport à cette référence. Cela se fait via une comparaison de chaque site avec son équivalent en position au niveau de la séquence de référence. Cependant il faut qu'au mieux possible les contenus que l'on compare soient comparables en étant non teintées par des variables apportées par la technique. C'est là qu'intervient l'analyse qualité. **

========================================================================== Nous allons parcourir un peu quelques travaux pour nous figurer le squelette de notre pipeline c'est à dire les grandes étapes qu'il faudra au moins

2 Pipelines sont rapportés par Ricke et al., 2018 : un de SAMtools et un autre de GATK.

aligner chaque individu à la référence en alignan les reads de son séquençage à celle-ci

convertir un fichier de sortie SAM en BAM

ordonner les fichiers BAM suivant les coordonnées

recréer un index au fichier BAM ordonné

• L'input de SAMtool mpileup : bam
• L'output de SAMtools mpileup : un fichier BCF (binary call format) ou vcf.

si un bcf est produit dans l'étape précédente, SAMtools calls (v1.3) est utilisé pour parser le BCF et faire la conversion en VCF.

output : un fichier VCF prduit à partir d'une liste de loci rsids En absence de rsid, les loci sont ajoutés manulmmenent

il n'appelle de variants que pour les loci dans le VCF d'entrée

• le nombre de fois qu'un locus est présent sur un brin

• Les runtimes des pipelines SAMtools et GATk sont comparées par Ricke et al., 2018 à un 3ème outil : SAMtools 7x plus rapide que GATK sous une puissance de calcul identique (3 threads) .
• GATK a arrété son support pour les plateformes Ion Torrent.
• Le pipeline SAMtools a de bonnes performances pour tous les SNP autosomaux. Cependant dans le cas où l'on aligne des reads de génomes humains à leur référence, ~50% des SNP du chr Y sont alignés au chr X pour ce qui est de la région pseudoautosomale (PAR).
• en tant que softs de calling, SAMtools (samtools mpileup) et GATK (UnifiedGenotyper) font ont les memes accomplissements (bam en entrée, vcf en sortie, realignement, recalibrage, singl comme multi-sample, calls=SNP+Indels,Linux). Cepednant, GATK supporte le multi-threading.

Dans les cas des deux pipelines proposés (celui avec SAMtools comme celui avec GATK), un parseur peut être construit pour renvoyer :

• la profondeur d couverture total
• les variatnts detectés (allèles) ==========================================================================

il inclut 5 parties jusqu'ici :

• analyse qualité des reads
• alignement des reads à un genome de reference
• post-traitemnt de l'alignement
• pré-traitement au calling
• calling ou la detection proprement dite des SNP

1 - FASTQC : FastQC is a tool to access the quality of reads (for example, to see where the quality scores start to significantly drop off and/or at what position that drop occurs). it does not access the quality of all reads though but only use a subsampling of the total. c'est globalement une analyse de la qualité des reads sur plusieurs tableaux. avantage = je le connais mieux et donc je peux mieux comprendre l'état de mes reads. fait ici pour voir la qualité des reads avant trimming.

cmd : fasqc filename.fq -o dir_to_put_the_output_files
options :
-o pour definir le dossier de sortie (nécessaire)

Inteprétation :

• how to spot the errors data? http://www.bioinformatics.babraham.ac.uk/training/Sequence_QC_Course/How%20sequencing%20experiments%20fail.pdf
• comment comprendr le gros de lanalyse : http://www.bioinformatics.babraham.ac.uk/training/Sequence_QC_Course/Developing%20and%20Implementing%20QC.pdf
• comprendre le niveau de duplication : https://www.biostars.org/p/107402/ 2 - CUTADAPT :
• pour le nettoyage de reads (avantage = je le connais mieux que d'autres comme trimmomatic) cutadapt enlève les bases respectant certaine caractéristiques fournies à l'outil. Un cas de ce travail est le trimming (raccourcir en coupant sur les extremités). les options -a (--adapter=ADAPTER) et -g (--front=ADAPTER) permettent par exemple de couper rspectivement en 3' et 5' des séquences d'adaptateurs et en même temps toutes les bases suivant vers la fin de cette extrémité. Avec $ et ^ (anchoring), il est possible de respectivement préciser si l'adaptateur est trimmé si et seulement si sa séquence finit ou débute celle du read. A ce point, des options plus puissantes existent et seront utilisées (voir ci dessous). voir (http://gensoft.pasteur.fr/docs/cutadapt/1.8.3/guide.html) pour une figure sur les positions 5', 3' sur les reads et comment apparait l'ancrage.

Cmd : cutadapt

options : 

Finding adapters:
-O MINLENGTH, --overlap=MINLENGTH
                        If the overlap between the read and the adapter is
                        shorter than MINLENGTH, the read is not modified.
                        Reduces the no. of bases trimmed due to random adapter
                        matches. Default: 3
-b ADAPTER, --anywhere=ADAPTER
                        Sequence of an adapter that may be ligated to the 5'
                        or 3' end (paired data: of the first read). Both types
                        of matches as described under -a und -g are allowed.
                        If the first base of the read is part of the match,
                        the behavior is as with -g, otherwise as with -a. This
                        option is mostly for rescuing failed library
                        preparations - do not use if you know which end your
                        adapter was ligated to!
Filtering of processed reads:
-m LENGTH, --minimum-length=LENGTH
                        Discard trimmed reads that are shorter than LENGTH.
                        Reads that are too short even before adapter removal
                        are also discarded. In colorspace, an initial primer
                        is not counted. Default: 0
Additional read modifications:
-q [5'CUTOFF,]3'CUTOFF, --quality-cutoff=[5'CUTOFF,]3'CUTOFF
                        Trim low-quality bases from 5' and/or 3' ends of each
                        read before adapter removal. Applied to both reads if
                        data is paired. If one value is given, only the 3' end
                        is trimmed. If two comma-separated cutoffs are given,
                        the 5' end is trimmed with the first cutoff, the 3'
                        end with the second.
see (http://gensoft.pasteur.fr/docs/cutadapt/1.8.3/guide.html) for the quality trimming algorithm works
For this to work correctly, the quality values must be encoded as ascii(phred quality + 33). If they are encoded as ascii(phred quality + 64), you need to add --quality-base=64 to the command line.
-u LENGTH, --cut=LENGTH
                        Remove bases from each read (first read only if
                        paired). If LENGTH is positive, remove bases from the
                        beginning. If LENGTH is negative, remove bases from
                        the end. Can be used twice if LENGTHs have different
                        signs.
Paired-end options:
    The -A/-G/-B/-U options work like their -a/-b/-g/-u counterparts, but
    are applied to the second read in each pair.

    -A ADAPTER          3' adapter to be removed from second read in a pair.
    -G ADAPTER          5' adapter to be removed from second read in a pair.
    -B ADAPTER          5'/3 adapter to be removed from second read in a pair.
    -U LENGTH           Remove LENGTH bases from second read in a pair (see
                        --cut).

Output:
    -o FILE, --output=FILE
                        Write trimmed reads to FILE. FASTQ or FASTA format is
                        chosen depending on input. The summary report is sent
                        to standard output. Use '{name}' in FILE to
                        demultiplex reads into multiple files. Default: write
                        to standard output
    -p FILE, --paired-output=FILE
                        Write second read in a pair to FILE. in the commande line, the name of each the 2 output fastq cleaned files have to be given before the name of each the 2 input fastq file that are to be cleaned. the -p is to be put between the name of 1st ouput and the 2nd output. also write the the 1st input file before the 2nd input file.

Essayer aussi Atropos :

3 - FASTQC : fait ici pour voir la qualité des reads après trimming

il s'agit de l'alignement des shorts reads obtenus au séquençage sur un génome dit référence. Il peut être un génome de référence comme il en existe actuellement pour beaucoup d'espèces ou un assemblage. L'objectif de l'alignement est de comparer un read donné à une séquence connue et d'enregistrer le resultat de cette comparaison. Pour la plupart des outils d'alignment, l'alignment commence par une 1ère étape d'indexation puis suit létape d'alignement proprement dit.

se fait pour la référence et/ou les reads à aligner. Les aligneurs les pls récents (ex : Bowtie2) indexent la référence car il est plus avantageux d'indexer une fois dans le cas de la réf que plusieurs fois , cas de reads. 2 types d'algorythmes d'indexation:

• ceux basés sur le tri des suffixes/prefixes : avantage = memory-efficient et plus rapide ; inconvenient = moins précis (ex : BWA utilise l'algo BWT, SOAP3 utilise la 2 BWT)
• ceux basés sur les tables de hashage : avantage = plus sensiblle; Inconvenient = plus lent (ex: BFAST utilise le spaced seed,SSAHA2 utilise le Q-gram filter and multiple seeds, Novoalign utilise le seed extension)
• NB : des exceptions. ex : Slider Software dédié aux sorties de séquenceur Illumina. il fait un merge-sorting des sous séquences de la ref et des reads.

2 algorithmes, indépendamment de la méthode d'indexation choisie :

• algo Smith-Waterman
• algo Needle-Wunsch

• Selon le soft utilisé, l'alignment pourra avoir des gaps ou non. Autoriser les gaps est préférable si par la suite l'alignement est utilisé pour detecter des réarrangments génomiques.
• Un bon soft pour l'alignement devrait supporter le paired-end mode pour un alignement d'une plus longue séquence et améliorer la qualité d'alignement.

=> le choix portera sur un algo basé sur tri des sufixes/préfixes pour privilégier memory efficient (x: BWA) avant sensibilité (Novoalign). BWA a 3 algo mais MEM est le plus récent et supporte les alignements splits en ayant à l'esprit l'idée de fournir aussi les fichiers bam à d'autres softs pour a detection de SV. NB : BWA MEM est rapporté comme le mieux adapté pour les reads entre 70 et 100 pb d'Illumina. Ceci est un avantage mais ne conaissant ni la taille ni l'origine de nos données encore cela reste une projection pour si de telles données sont utilisées.

• NB: il srait intéressant de voir ce que donne une comparaison de Bowtie (-D 20 -R 3 -N 1 -L 20) et bwa mem (apparemment 91% vs 85%)

=> on a choisit BWA MEM:

4 - BWA MEM : bon compromis entre vitesse et sensibilité (% reads alignés) (output = SAM). ne pas oublier d'ajouter les readgroups (-R str):

cmd : bwa mem
=> etape 1/2 : faire un bwa index reference.fasta pour créer un index 

=> etape 2/2 : les options de BWA MEM que l'on devrait utiliser : 
 - -t 2 : fixer le nombre de threads utilisé pour la tache, ici 2 mais par défaut 1
 - -T 30 : score minimal qu'un alignement doit avoir pour aller dans l'output, ici 30 mis par défaut 30 aussi.
 - -M : considérer tout alignement aditionnel à part l'alignement primaire comme un alignement secondaire . the thing about the -M option is that when you want the supplementary segments of a plit to be considered dont use the option (an interesting case would be a sv detection study). Picardtools mak duplicates used to have a problem with supplementary alignements so we use the option to make them all secondary (manually it's setting the secondary flag 0x100 for all the alignments that have the supplementary flag 0x800)
 - fournir la reference (un .fa le plus souvent). en tant que tel, bwa mem n'a pas besoin de la référence mais plutot des 5 fichirs d'index produits juste à coté de la référence par bwa index.
 - fournir le fichier contenant les reads 1 du paired-end
 - forunir le fichier contenant les reads 2 du paired-end
 - -R : -R '@RG\tID:AX227\tSM:AX227\tPL:Illumina' permet de tagger (marquer) nos readgroups d'un ID (le plus souvent le nom de l'individu), d'un SM (le nom de l'échantillon) et d'un PL (nom de la platforme de séquençage ayant été utilisé pour séquencer l'échantillon et avoir ces reads). 
 - avec > fichier de sortie.sam : donner le fichier de sortie qui sera un sam

4 étapes faisables avec SAMtools :

• 5 - Conversion du format SAM au format BAM : genère un BAM (binary alignment map avec une taille réduite) depuis le SAM d'entrée. avantage = reduction taille des fichiers, standardisation à un même format pour les variant callers. Nous en profitons également pour restreindre les alignements à ceux où les reads sont bien alignés et écarter les ducplicats PCR ou optiques

Cmd : samtools view 

options : 
 - -b: Convertir en fichier bam
 - -S: Le fichier d'entrée est un fichier sam
 - -h: Inclure le header (-H pour inclure tout sauf le header)
 - -f 0x2 : inclure les alignés n proper pair
 - -F 0x400 : ecarter les ducplicats (PCR comme optiques)
 - le fichier de sortie est défini avant le fichier d'entrée.
 - -o fichier_de_sortie.bam fichier_dentree.sam

• 6 - Ordonner le bam obtenu : ordonner les fichiers BAM suivant les coordonnées pour accélerer les étapes à suivre

Cmd : samtools sort fichier_obtenu.bam version_ordonnée_du_fichier.bam

• 7 - indexer le bam ordonné : dans la suite, des outils vont parcourir le bam pour y chercher des différences d'alignements à chaque position. cela peut être accélré en créant un index pour ce bam

Cmd : samtools index version_ordonnée_du_fichier.bam

• NB : on aurait pu faire les étapes 6 et 7 en une seule commande avec PicardtoolsSortSam (avoir un bam trié et l'index en meme temps). faire :

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

• Facultatif : Enlever les duplicats : les duplicats sont des paires de reads présentes plusieurs fois à travers les données au lieu d'une fois. Leur impact est qu'en ne les enlevant pas la profondeur de couverture est faussée. enlever les duplicats de la PCR sinon la couverture est faussée or celle-ci est utilisée par les softs de calling et la conséquence serait des genotypes surestimés (dans certains pipelines pour être sur de très tot enlever, on peut le faire dès l'étape de samtoolsview convertissant un sam en bam; si on veut les garder jusqu'à plus tard pt mark duplicates peut les enlever)

cmd : picardToolsMarkDuplicates REMOVE_DUPLICATES=true 

ou aussi

Cmd : samtools rmdup

NB : les duplicats peuvent aussi être enlevés avec PicardtoolsMarkDuplicates plus tard. cela nous permet de les avoir dans notre rapport de stats d'alignements et de les enlever à un moment plus proche de ce qu'il vont perturber (snp calling proprement dit et genotypage)

• 8 - Rapport de statistiques d'alignements (automatique avec certains aligners (MOSAIK, Bowtie, SOAP2), mais pour BWA faire un SAMtools flagstat). l'analyse porte sur un .bam (sortie de la conversion avec samtools view -b). ce rapport permet entre autres d surveiller le % de reads alignés (s'il est bas, la référence est soit incomplète ou alors qu'il y'a une grande distance génomique entre les reads alignés et la référence c'est à dire qu'on est entre de comparer 2 inds relativement distants). genere une fiche de stats sur l'alignement afin de vérifier sa qualité avant de continuer. surveiller le % de reads alignés : s'il est bas, sequençage de basse qualité pour les reads utilisés ou référence incomplete ou reférence et reads non concordants (essayer d'aligner des reads d'un ind A sur un génome référence d'un ind B très éloigné)

Cmd : samtools flagstat version_ordonnée_du_fichier.bam > stats_report.samtoolsFlagstat

NB : les commandes de GATK (soft utilisé dans la suite) sont introduites par une version raccourci du nom du soft dans gatk pour des souci de présentation. ils devraient être introduits par java -jar GenomeAnalysisTK.jar

• réalignement local : des artefacts au cours de l'alignement peut mener à un mismatch à certaines positions et donc de faux-positifs en SNP calling. Des outils permettnt des realignement locaux (inclus dans GATK). Explication : l'occurrence d'un indel chez l'individu séquencé, a pour conséquence des reads qui peuvent en début ou fin couvrir cet indel. alors a l'alignement, meme si ces reads sont correctement alignés, les autres bases du read peuvent etre en mismatch. ces mismatchs s'il ne sont pas corrigés en les minimisant au moins vont ressortir comme des SNPs alors qu'ils ne le sont pas t sont juste des bases décalées par rapport à celles contenues dans la référence. l'alignement local est fait pour cela et ce sera fait pour les intervalles sur les reads qui sont mal alignés à cause d'indels.

NB : Si dans la suite d'un tel workflow de detection de variants, il est utilisé un caller faisant une etape d'assemblage avant la detection (ex : HaplotypeCaller ou MuTect2), il n'est pas nécessaire de faire ce réalignement. Nous allons bien utiliser HC ici pour le calling de SNPs cependant nous traiterons cette étape car si notre caller est changé pour un caller de type legacy (ex : UnifiedGenotyper, Mutect 1ère version)

9 - gatkRealignerTargetCreator : realignement (etape 1/2).Entrée : un ou plusieurs bam. determiner les intervales suspectés d'avoir besoin de realignement puis les mettre dans Sortie : un fichier à fournir au realigner dans la prochaine etape. Limites : 1 - si plusieurs indels sont trouvés, l'outil en choisi un seul, celui qui est le plus plausible d'avoir causé le misalignement 2- les reads de la technologie 454 contient de manière inhérente de faux indels. le realigner ne fonctionnra pas sur les reads de tels outils. 3 - d'autres reads ignorés sont ceux ayant MQ0 et ceux ayant des operateurs d'indels consecutifs dans leur CIGAR string NB : Reads filters and multi-threading option can be given. Read filters are applied to the data before the processing by the RelignerTargetCreator. see the documentation (https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_indels_RealignerTargetCreator.php) for on overview on them

Cmd : 
java -jar GenomeAnalysisTK.jar \
   -T RealignerTargetCreator \
   -R reference.fasta \
   -I input.bam \
   -o forIndelRealigner.intervals
optional : 
--known indels.vcf : Input VCF file with known indels. Any number of VCF files representing known SNPs and/or indels. Could be e.g. dbSNP and/or official 1000 Genomes indel calls. SNPs in these files will be ignored unless the --mismatchFraction argument is used. This argument supports reference-ordered data (ROD) files in the following formats: BCF2, VCF, VCF3
-maxIntervalSize (500) : maximum interval size; any intervals larger than this value will be dropped. Because the realignment algorithm is N^2 (O(N2)), allowing too large an interval might take too long to completely realign.
-minReads (4) : minimum reads at a locus to enable using the entropy calculation
-mismatch (0.0) : fraction of base qualities needing to mismatch for a position to have high entropy. To disable this behavior, set this value to <= 0 or > 1. This feature is really only necessary when using an ungapped aligner (e.g. MAQ in the case of single-end read data) and should be used in conjunction with '--model USE_SW' in the IndelRealigner. 
-window (10) : window size for calculating entropy or SNP clusters. Any two SNP calls and/or high entropy positions are considered clustered when they occur no more than this many basepairs apart. Must be > 1.

• how does the target creator works : see at (https://gatkforums.broadinstitute.org/gatk/discussion/6284/realignertargetcreator-entropy-values-in-output)

looking for mismatches using this:
if ( p.getBase() != refBase )
    mismatchQualities += p.getQual();
    totalQualities += p.getQual();
and then
if ( pileup.getNumberOfElements() >= minReadsAtLocus && (double)mismatchQualities / (double)totalQualities >= mismatchThreshold )
    hasPointEvent = true;
and it separately looks for indels using:
if ( p.isDeletion() || p.isBeforeInsertion() ) {
    hasIndel = true;
    if ( p.isBeforeInsertion() )
        hasInsertion = true;
}
And finally, if any of **hasPointEvent**, **hasIndel** or **hasInsertion** is true, then the tool outputs the coordinates of the interval to the output file. So its understandable that it's not really set up to output an entropy score -- and indeed there is only a score for the mismatch case, not for the indel cases. So it would take some dev work to get this working to some extent.

• a possible way to find potentially off-target regions for CRISPR/Cas9 WGS data : 1 - using the RealignerTargetCreator to output the entropy scores would require dev time, instead use HC with a special option :

• A better option is to use a subset of the HaplotypeCaller machinery instead.
• HC has a hidden argument called --justDetermineActiveRegions which does just that: it runs only the entropy calculations (Percentage of mismatches at a locus ) and doesn't do any of the subsequent assembly, haplotype scoring or genotyping steps. There is also a documented argument called --activityProfileOut that allows you to output the activity profile, ie the per base probability that there is something going on at that site. from there it's a question of parsing the output to be useable in your work.

10 - gatkIndelRealigner : realignement (etape 2/2), lancer le realigneur sur ces intervales

java -jar GenomeAnalysisTK.jar -T IndelRealigner -R reference.fa -I sorted_bwamem_output_file.bam -targetIntervals forIndelRealigner.intervals -o realignedBam.bam 

Cmd : 
java -jar GenomeAnalysisTK.jar \
   -T IndelRealigner \
   -R reference.fasta \
   -I input.bam \
   -targetIntervals forIndelRealigner.intervals \
   -o realignedBam.bam
NB : 
1 - le(s) bam d'entrée, la ref, le indels.vcf pour les indels connus doivent être ceux utilisés précédemment à l'étape RealignerTargetCreator. Le fichier .intervals utilisé doit etre output de l'étape précédente. 
2 - les reads 454 ou d'autres technologies induisant de faux indels de manière inhérente ne sont pas supportés.
3- les reads à MQ0 et ceux avec des operateurs consecutifs d'indels dans leur CIGAR string ne sont pas supportés aussi

Options : 
-known indels_to_keep.vcf -model USE_READS -LOD 5.0 --nWayOut -entropy 0.15 --maxConsensuses 30 -maxsize 3000 -maxPosMove 200 -greedy 120 -maxReads 20000 -maxInMemory 150000 -notags false

--consensusDeterminationModel / -model : Determines how to compute the possible alternate consenses
We recommend that users run with USE_READS when trying to realign high quality longer read data mapped with a gapped aligner; Smith-Waterman is really only necessary when using an ungapped aligner (e.g. MAQ in the case of single-end read data). The --consensusDeterminationModel argument is an enumerated type (ConsensusDeterminationModel), which can have one of the following values:
KNOWNS_ONLY : Uses only indels from a provided ROD of known indels.
USE_READS : Additionally uses indels already present in the original alignments of the reads.
USE_SW : Additionally uses 'Smith-Waterman' to generate alternate consenses.

--LODThresholdForCleaning / -LOD
LOD threshold above which the cleaner will clean
This term is equivalent to "significance" - i.e. is the improvement significant enough to merit realignment? Note that this number should be adjusted based on your particular data set. For low coverage and/or when looking for indels with low allele frequency, this number should be smaller.

--nWayOut / -nWayOut
Generate one output file for each input (-I) bam file (not compatible with -output)
Reads from all input files will be realigned together, but then each read will be saved in the output file corresponding to the input file that the read came from. There are two ways to generate output bam file names: 1) if the value of this argument is a general string (e.g. '.cleaned.bam'), then extensions (".bam" or ".sam") will be stripped from the input file names and the provided string value will be pasted on instead; 2) if the value ends with a '.map' (e.g. input_output.map), then the two-column tab-separated file with the specified name must exist and list unique output file name (2nd column) for each input file name (1st column). Note that some GATK arguments do NOT work in conjunction with nWayOut (e.g. --disable_bam_indexing).

--entropyThreshold / -entropy
Percentage of mismatches at a locus to be considered having high entropy (0.0 < entropy <= 1.0)
For expert users only! This is similar to the argument in the RealignerTargetCreator walker. The point here is that the realigner will only proceed with the realignment (even above the given threshold) if it minimizes entropy among the reads (and doesn't simply push the mismatch column to another position). This parameter is just a heuristic and should be adjusted based on your particular data set.

--maxConsensuses / -maxConsensuses
Max alternate consensuses to try (necessary to improve performance in deep coverage)
For expert users only! If you need to find the optimal solution regardless of running time, use a higher number.

--maxIsizeForMovement / -maxIsize
maximum insert size of read pairs that we attempt to realign
For expert users only!

--maxPositionalMoveAllowed / -maxPosMove
Maximum positional move in basepairs that a read can be adjusted during realignment
For expert users only!

--maxReadsForConsensuses / -greedy
Max reads used for finding the alternate consensuses (necessary to improve performance in deep coverage)
For expert users only! If you need to find the optimal solution regardless of running time, use a higher number.

--maxReadsForRealignment / -maxReads
Max reads allowed at an interval for realignment
For expert users only! If this value is exceeded at a given interval, realignment is not attempted and the reads are passed to the output file(s) as-is. If you need to allow more reads (e.g. with very deep coverage) regardless of memory, use a higher number.

--maxReadsInMemory / -maxInMemory
max reads allowed to be kept in memory at a time by the SAMFileWriter
For expert users only! To minimize memory consumption you can lower this number (but then the tool may skip realignment on regions with too much coverage; and if the number is too low, it may generate errors during realignment). Just make sure to give Java enough memory! 4Gb should be enough with the default value.

--noOriginalAlignmentTags / -noTags
Don't output the original cigar or alignment start tags for each realigned read in the output bam

• recalibrage du score de qualité de base: Le Phred score es la forme sous laquelle le sequencur exprime la qualité de la base. Pour que ce Phred score soit une estimation plus précise de ce qu'est la nouvelle probabilité de mismatching après le(s) realignement(s) locaux (fait dans GATK avec cmd recalibrate)

11 - gatkBaseRecalibrator : recalibrage qualité des bases (etape 1/2) : Variant calling algorithms rely heavily on the quality scores assigned to the individual base calls in each sequence read. These scores are per-base estimates of error emitted by the sequencing machines. Unfortunately the scores produced by the machines are subject to various sources of systematic technical error, leading to over- or under-estimated base quality scores in the data. Base quality score recalibration (BQSR) is a process in which we apply machine learning to model these errors empirically and adjust the quality scores accordingly. first the program builds a model of covariation based on the data and a set of known variants (which you can bootstrap if there is none available for your organism), then it adjusts the base quality scores in the data based on the model. This tool performs the first step described above: it builds the model of covariation and produces the recalibration table. It operates only at sites that are not in the vcf supplied for positions to maintain as they are. we assume that all reference mismatches we see are therefore errors and indicative of poor base quality. This tool generates tables based on various user-specified covariates (such as read group, reported quality score, cycle, and context). Assuming we are working with a large amount of data, we can then calculate an empirical probability of error given the particular covariates seen at this site, where p(error) = num mismatches / num observations. The output file is a table (of the several covariate values, number of observations, number of mismatches, empirical quality score).

cmd :  java -jar GenomeAnalysisTK.jar -T baseRecalibrator -R reference.fa -I realignedBam.bam -o recal_data.table 
optional params : 
--no_standard_covs/-noStandard : Do not use the standard set of covariates, but rather just the ones listed 
using the -cov argument
--covariate/ -cov : One or more covariates to be used in the recalibration. Can be specified multiple times
--num_cpu_threads_per_data_thread / -nct : Number of CPU threads to allocate per data thread.
Each CPU thread operates the map cycle independently, but may run into earlier scaling problems with IO than data threads. Has the benefit of not requiring X times as much memory per thread as data threads do, but rather only a constant overhead. 
--lowMemoryMode : Reduce memory usage in multi-threaded code at the expense of threading efficiency
 
--bqsrBAQGapOpenPenalty/-bqsrBAQGOP : BQSR BAQ gap open penalty (Phred Scaled). Default value is 40. 30 is perhaps better for whole genome call sets

12 - gatkPrintReads : recalibrage qualité des bases (etape 2/2) PrintReads is a generic utility tool for manipulating sequencing data in SAM/BAM format. It can dynamically merge the contents of multiple input BAM files, resulting in merged output sorted in coordinate order. It can also optionally filter reads based on various read properties such as read group tags using the --read_filter/-rf command line argument (see documentation on read filters for more information).

Note that when PrintReads is used as part of the Base Quality Score Recalibration workflow, it takes the --BQSR engine argument, which is listed under Inherited Arguments > CommandLineGATK below. --BQSR / -BQSR Input covariates table file for on-the-fly base quality score recalibration Enables on-the-fly recalibrate of base qualities, intended primarily for use with BaseRecalibrator and PrintReads (see Best Practices workflow documentation). The covariates tables are produced by the BaseRecalibrator tool. Please be aware that you should only run recalibration with the covariates file created on the same input BAM(s) or CRAM(s).

Cmd : java -jar GenomeAnalysisTK.jar -T PrintReads -R reference.fa -I realignedBam.bam -o recalibrated.bam --BSQR recal_data.table
optionals params : 
--read_filter/ -rf : MappingQualityZero and others bad reads
-n : Print the first n reads from the file, discarding the rest
-dfrac DEC : dowsamples the bam file to DEC%
--platform : Exclude all reads with this platform from the output
--readgroup : Exclude all reads with this read group from the output

13 - picardToolsMarkDuplicates (REMOVE_DUPLICATES=true) :

marque et enlève les duplicats de paires de reads pour que dans la suite on puisse travailler avec des paires uniques. Avantage = cela rectifie la profondeur de couverture augmentée par les duplicats or celle-ci est utilisée pour l'estimation des genotypes. Expliquée d'une autre manière, les models statistiques utilisés par les callers se basent sur l'indépendance : si les duplicats proviennent d'artefacts de PCR alors cette notion d'indépendance s'écroule.

• Solutions de SAMtools pour realignement et recalibrage : SAMtools assigne à chaque alignement un base alignement quality (BAQ) score. Il est bas pour une base en mismatch dû à un artefact et sera plus grand en realignant localement

La detection proprement dite des SNP elle est faite par des SNP ou genotype callers. 2 types:

• ceux utilisant une approche heuristique : variant calling basé sur de multiple source d'informations concernant la structure et la qualité des données. Moins utilisés du fait de leur demande en calcul. Ex : VarScan2. Un caller basé sur un algo heuristique determine le genotype en se basant sur des seuils pour la profondeur de couv (minimum 33), la qualité de la base (minimum 20) et la fréquence de l'allèle du variant (minimum 0.08); cela est suivi par un test exact de Fisher sur le nombre de reads supportant chaque allèle en comparant avec la distribution attendue si l'on considère les erreurs de séqueçage seulement
• ceux utilisant une approche probabiliste : basée sur des mesures statistiques reflétant l'incertitude des genotypes estimés. Plus d'iformations telles les frequences allèliques peuvent être incluses dans l'analyse. le calling est ici basé sur un calcul de vraisemblance pavec le théorème de Bayes. Après l'étape de realignement et recalibration en prétraitement du calling, les vraisemblances de chaque genotype possible (homozygote de l'allèle de ref, homozygote d'un allèle alternatif ou heterozygote) sont calculées en se basant sur les scores de qualité et le decompte d'un allèle parmi les reads à ce site de SNP; en se basant sur le framework Bayesien, la vraismblance calculée est combinée à une probabilité antérieure au genotypage pour donner une probabilité postérieure au genotypage; le genotype ayant la pls haute probabilité postérieure au genotypage est choisi. Ex: SAMtools et GATK sont tous deux basés sur cette approche probabiliste utilsant le principe Bayesien. samtools mpileup les infos dans un BAM et calcule la vraisemblance des données pour chaque gentype possible; bcftools applique la probabilité postérieure au gentypage et fait le SNP calling

14 - gatkHaplotypeCaller (-T HaplotypeCaller): calling proprement dit. HC a de meilleures performances eentermes de nombre de detection de SNP et Indels que UnifiedGenotypeCaller. Sortie = un vcf non filtré contenant tous les variants trouvé sans restriction de type, de score de qualité sur sites ou de d'échantillons

• compréhension : (https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_haplotypecaller_HaplotypeCaller.php)

1. Define active regions
The program determines which regions of the genome it needs to operate on, based on the presence of significant evidence for variation.


2. Determine haplotypes by assembly of the active region
For each ActiveRegion, the program builds a De Bruijn-like graph to reassemble the ActiveRegion, and identifies what are the possible haplotypes present in the data. The program then realigns each haplotype against the reference haplotype using the Smith-Waterman algorithm in order to identify potentially variant sites.


3. Determine likelihoods of the haplotypes given the read data
For each ActiveRegion, the program performs a pairwise alignment of each read against each haplotype using the PairHMM algorithm. This produces a matrix of likelihoods of haplotypes given the read data. These likelihoods are then marginalized to obtain the likelihoods of alleles for each potentially variant site given the read data.


4. Assign sample genotypes
For each potentially variant site, the program applies Bayes' rule, using the likelihoods of alleles given the read data to calculate the likelihoods of each genotype per sample given the read data observed for that sample. The most likely genotype is then assigned to the sample.

cmd : 
optional params :
--sample_ploidy/-ploidy (2) : Ploidy per sample. For pooled data, set to (Number of samples in each pool * Sample Ploidy).
--maxReadsInRegionPerSample 10 000 : Maximum reads in an active region
--min_base_quality_score/ -mbq (10) : Minimum base quality required to consider a base for calling
--standard_min_confidence_threshold_for_calling/-stand_call_con (10.0): The minimum phred-scaled confidence threshold at which variants should be called
--activeRegionIn
 -AR : Use this interval list file as the active regions to process
--comp : Comparison VCF file
--bamOutput/-bamout : File to which assembled haplotypes should be written
--kmerSize (10,25) : Kmer size to use in the read threading assembler
--max_alternate_alleles/-maxAltAlleles (6) : Maximum number of alternate alleles to genotype
--max_genotype_count/-maxGT (1024) : Maximum number of genotypes to consider at any site
--maxNumHaplotypesInPopulation (128) : Maximum number of haplotypes to consider for your population
--maxReadsInMemoryPerSample (30000) : Maximum reads per sample given to traversal map() function
--maxTotalReadsInMemory (10000000) : Maximum total reads given to traversal map() function
--output_mode/-out_mode : EMIT_VARIANTS_ONLY : Which type of calls we should output
NB : see the assembly related options

• Restriction des résultats du calling*

15 - gatkVariantFiltration (-T VariantFiltration couples (--filterName " " --filterExpression " ")) : Filter variant calls based on INFO and FORMAT annotations. ajoute à la colonne FILTER un "pass" si l'entrée respecte certains critères tel que le score de qualité des sites, la profondeur (nbre de reads) au site, etc . (https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_filters_VariantFiltration.php). DP : nombre de reads à la position (depth) et AD : profondeur de allèle reference et alternatif (ex: DP>=20 on peut avoir AD=10,10.

Cmd :
 java -jar GenomeAnalysisTK.jar \
   -T VariantFiltration \
   -R reference.fasta \
   -o output.vcf \
   --filterExpression "AB < 0.2 || MQ0 > 50" \
   --filterName "SomeFilterName" 
 
optionals params : 
--variant / -V
Input VCF file
Variants from this VCF file are used by this tool as input. The file must at least contain the standard VCF header lines, but can be empty (i.e., no variants are contained in the file).

This argument supports reference-ordered data (ROD) files in the following formats: BCF2, VCF, VCF3
--clusterSize / -cluster (3)
The number of SNPs which make up a cluster
Works together with the {@code --clusterWindowSize} argument.
--clusterWindowSize / -window (0)
The window size (in bases) in which to evaluate clustered SNPs
Works together with the {@code --clusterWindowSize} argument. To disable the clustered SNP filter, set this value to less than 1.
--filterExpression / -filter
One or more expression used with INFO fields to filter
VariantFiltration accepts any number of JEXL expressions (so you can have two named filters by using {@code --filterName One --filterExpression "X < 1" --filterName Two --filterExpression "X > 2"}).
--filterName / -filterName 
Names to use for the list of filters
This name is put in the
FILTER
field for variants that get filtered. Note that there must be a 1-to-1 mapping between filter expressions and filter names.
--genotypeFilterExpression / -G_filter
One or more expression used with FORMAT (sample/genotype-level) fields to filter (see documentation guide for more info)
Similar to the

INFO
field based expressions, but used on the
FORMAT
(genotype) fields instead. {@link VariantFiltration} will add the sample-level
FT
tag to the
FORMAT
field of filtered samples (this does not affect the record's
FILTER
tag). One can filter normally based on most fields (e.g. {@code "GQ < 5.0"}), but the
GT
(genotype) field is an exception. We have put in convenience methods so that one can now filter out hets ({@code "isHet == 1"}), refs ({@code "isHomRef == 1"}), or homs ({@code "isHomVar == 1"}). Also available are expressions {@code isCalled}, {@code isNoCall}, {@code isMixed}, and {@code isAvailable}, in accordance with the methods of the {@link Genotype} object.
--genotypeFilterName / -G_filterName
Names to use for the list of sample/genotype filters (must be a 1-to-1 mapping); this name is put in the FILTER field for variants that get filtered
Similar to the

INFO
field based expressions, but used on the
FORMAT
(genotype) fields instead.
--mask / -mask
Input ROD mask
Any variant which overlaps entries from the provided mask rod will be filtered. If the user wants logic to be reversed, i.e. filter variants that do not overlap with provided mask, then argument {@code -filterNotInMask} can be used. Note that it is up to the user to adapt the name of the mask to make it clear that the reverse logic was used (e.g. if masking against Hapmap, use {@code -maskName=hapmap} for the normal masking and {@code -maskName=not_hapmap} for the reverse masking).

This argument supports reference-ordered data (ROD) files in the following formats: BCF2, BEAGLE, BED, BEDTABLE, EXAMPLEBINARY, RAWHAPMAP, REFSEQ, SAMPILEUP, SAMREAD, TABLE, VCF, VCF3
--filterNotInMask / -filterNotInMask
Filter records NOT in given input mask.
By default, if the {@code -mask} argument is used, any variant falling in a mask will be filtered. If this argument is used, logic is reversed, and variants falling outside a given mask will be filtered. Use case is, for example, if we have an interval list or BED file with "good" sites. Note that it is up to the user to adapt the name of the mask to make it clear that the reverse logic was used (e.g. if masking against Hapmap, use {@code -maskName=hapmap} for the normal masking and {@code -maskName=not_hapmap} for the reverse masking).
--maskName / -maskName
The text to put in the FILTER field if a 'mask' rod is provided and overlaps with a variant call
When using the {@code -mask} argument, the {@code maskName} will be annotated in the variant record. Note that when using the {@code -filterNotInMask} argument to reverse the masking logic, it is up to the user to adapt the name of the mask to make it clear that the reverse logic was used (e.g. if masking against Hapmap, use {@code -maskName=hapmap} for the normal masking and {@code -maskName=not_hapmap} for the reverse masking).
--maskExtension / -maskExtend (0)
How many bases beyond records from a provided 'mask' rod should variants be filtered
--invalidatePreviousFilters / NA (false)
Remove previous filters applied to the VCF
Invalidate previous filters applied to the {@link VariantContext}, applying only the filters here.
--invertFilterExpression / -invfilter (false)
Invert the selection criteria for --filterExpression
Invert the selection criteria for {@code --filterExpression}.
--invertGenotypeFilterExpression / -invG_filter (false)
Invert the selection criteria for --genotypeFilterExpression
Invert the selection criteria for {@code --genotypeFilterExpression}.
--missingValuesInExpressionsShouldEvaluateAsFailing / NA
When evaluating the JEXL expressions, missing values should be considered failing the expression
By default, if JEXL cannot evaluate your expression for a particular record because one of the annotations is not present, the whole expression evaluates as

PASS
ing. Use this argument to have it evaluate as failing filters instead for these cases.
--setFilteredGtToNocall / NA
Set filtered genotypes to no-call
If this argument is provided, set filtered genotypes to no-call (./.).

16 - gatkSelectVariants (-T SelectVariants -selectType SNP -selectType MNP) : generer un vcf ne contenant que les vrais SNP et les MNP (multi nucleotide polymorphisms) (https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_variantutils_SelectVariants.php)

-for more, see the tool page (https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_filters_VariantFiltration.php#--setFilteredGtToNocall)

Cmd : 
optionals params : 

• Allele depth (AD) and depth of coverage (DP). These are complementary fields that represent two important ways of thinking about the depth of the data for this sample at this site. AD is the unfiltered allele depth, i.e. the number of reads that support each of the reported alleles. All reads at the position (including reads that did not pass the variant caller’s filters) are included in this number, except reads that were considered uninformative. Reads are considered uninformative when they do not provide enough statistical evidence to support one allele over another.
• DP is the filtered depth, at the sample level. This gives you the number of filtered reads that support each of the reported alleles. You can check the variant caller’s documentation to see which filters are applied by default. Only reads that passed the variant caller’s filters are included in this number. However, unlike the AD calculation, uninformative reads are included in DP. See the Tool Documentation for more details on AD (DepthPerAlleleBySample) and DP (Coverage) for more details
• Reads that are not used for calling are not counted in the DP measure, but are included in AD
• ambiguously mapped reads (mapping quality = 0)indicative of high reference sequence homology, as such regions are known toproduce many false positives;
• shading of mismatched bases by read base quality, as clusters of bases of low quality can be indicative of sequencing errors;

GATK | Doc #4260 | Phred-scaled Quality Scores (https://software.broadinstitute.org/gatk/documentation/article.php?id=4260)

-conseils (https://www.researchgate.net/post/What_are_the_essential_SNP_filtering_criteria_one_should_follow_before_performing_GWAS_analysis)

• Evaluer la profondeur et l'uiliser pour filtrer les résukats du calling : http://www.gettinggeneticsdone.com/2014/03/visualize-coverage-exome-targeted-ngs-bedtools.html
• un pipeline éliminant les SNP proches car penant qu'ils sont issus de recombinaison...not understood and feared to be old resultats now outdated. to be reviewed. (http://userweb.eng.gla.ac.uk/cosmika.goswami/snp_calling/SNPCalling.html)
• essayer ce tuto (http://ddocent.com/filtering/)
• How can I filter a vcf filter a VCF file on minimum genotype depth and genotype quality for each sample. criterias sample.DP > 10 and sample.GQ > 15 :

bcftools view  -i  'MIN(FMT/DP)>10 & MIN(FMT/GQ)>15'   my.vcf.gz

(https://gatkforums.broadinstitute.org/gatk/discussion/4017/what-is-a-gvcf-and-how-is-it-different-from-a-regular-vcf) :

• les gVCF sont des VCF qui rapportent l'information liée à chaque position de la région séquencée et étudiée (intervalle d'intéret) et cela que cette position soit variante ou pas. Ils sont similaires au VCF obtenu en utilisant l'option --output_mode EMIT_ALL_SITES De UnifiedGenotyper sauf qu'on utilisera HC pour les avoir. Ils sont utilisés pour les analyses de cohortes. 2 types de gVCF can be obtained with HC :

• With BP_RESOLUTION, you get a gVCF with an individual record at every site: either a variant record, or a non-variant record.
• With GVCF, you get a gVCF with individual variant records for variant sites, but the non-variant sites are grouped together into non-variant block records that represent intervals of sites for which the genotype quality (GQ) is within a certain range or band. The GQ ranges are defined in the ##GVCFBlock line of the gVCF header. The purpose of the blocks (also called banding) is to keep file size down, and there is no downside for the downstream analysis, so we do recommend using the -GVCF option.

• use this workflow for cohort analysis : the principle is to analyse every bams of a sample until gvcf is obtained, join all sammple's bam and obtain one bam that you filter https://software.broadinstitute.org/gatk/documentation/article.php?id=3893

• vtnormalise and vcfeval in RTG tools sont les meilleures options actuelles pour cela.(https://doi.org/10.1016/j.csbj.2018.01.003). bcftools, gatk et vt ont été comparés pour la normalisation : vt is best . vt et bcftools ont la meme normalisation mais vt a eu l'algo en 1er. gatk ne normalise bien que les biallelic; les multiallelic ou autres variants mixés (des SNP adjacents à des indels) ne peuvent être alignés à gauche par gatk (de tels cas après réalignemnt avec vt par ex peuvent se rèveler n'etre que de simples indels d'où l'importance de la normalisation)
• vt : (https://github.com/atks/vt) et its wiki (https://genome.sph.umich.edu/wiki/Vt) voir ici Vt (outil 1)
• algorithmique de la normalisation : un vcf peut contenir plusieurs repésentations d'une meme variation. Ces variantions peuvent etre reduites en les normalisant pour atteindre une form canonique et puis enlevant par la suite les doublons. La normalisation effectue l'alignement à gauche et la parcimonie des variants ds le vcf. L'algorithme suivi par vt pour installer cet alignement à gauche et cet parsimonie a deux parties : ---Lines 1 to 8 performs the left alignment and ensures parsimonious representation on the right side. ---Lines 9 to 11 ensures parsimonious representation on the left side. voir figure ci dessous or explication ici (https://genome.sph.umich.edu/wiki/Variant_Normalization) :

*normaliser : voir les regles ici (https://genome.sph.umich.edu/wiki/Variant_Normalization)

• NB1 : revoir les variants communs comme SINE, LINE, ALU (https://fr.wikipedia.org/wiki/S%C3%A9quence_r%C3%A9p%C3%A9t%C3%A9e) et numt (https://en.wikipedia.org/wiki/NUMT)

• essayer ce tuto (https://wiki.bits.vib.be/index.php/NGS_Exercise.6)
• highlights changes and cases of things filtered by IGV: (https://www.mendeley.com/research-papers/variant-review-integrative-genomics-viewer/?utm_source=desktop&utm_medium=1.19.1&utm_campaign=open_catalog&userDocumentId=%7B3fef3af2-75ea-4f82-a0b5-b2517fa4d5db%7D) or more at (http://software.broadinstitute.org/software/igv/). you can get IGV here (https://www.igv.org/)

-essaer ce tuto (https://rdrr.io/github/jjellis/GenomicVis/man/vcf.venn.html)

• The genome is a large structure with localized regions of interest, frequently separated by large oceans of uninteresting sequence. To help visualize data in this context, Circos can create images with variable axis scaling, permitting local magnification of genomic regions to be controlled without cropping. Scale smoothing ensures that the magnification level changes smoothly. In combination with axis breaks and custom ideogram order, the final image can be easily tuned to offer the clearest illustration of your data. http://circos.ca/intro/genomic_data/img/circos-conde-nast-large.png
• The graphic shows the human genome annotated with data related to genes implicated in disease, regions of variation found in various populations, and regions of similarity between chromosomes.

Over the course of the 1000 Genomes Project, how variants were called from the samples changed quite dramatically. Two clear lessons from the project were, when considering low coverage data, calling from multiple samples at once produces more, higher quality variants and considering the sites discovered from multiple algorithms improved the discovery rate and accuracy of discovery. Many different programs and strategies were developed over the duration of the project. The publications referred to at the top of this page are the best place to get a description of what programs were used and how the 1000 Genomes variant calling pipeline was run (http://www.internationalgenome.org/analysis)

-atropos vs cutadapt -minimap2 vs bwa mem -SNPf -SnpSift : (http://snpeff.sourceforge.net/SnpSift.html)

• bedtools : http://bedtools.readthedocs.io/en/latest/
• comparaison de 2 vcfs : https://www.reddit.com/r/bioinformatics/comments/5jqqsj/comparing_variants_from_two_vcfs/
• comparer avec RTG tools (https://github.com/RealTimeGenomics/rtg-tools) :

• GATK : Genome Analysis Toolkit

Les softs utilisés pour un pipeline de detection de SNP sont nombreux et changent très rapidement. L'important est de surtout comprendre les différences entre algorithmes. Ceci permet de customiser un pipeline déja créé selon le but de l'analyse, la taille de nos reads après nettoyage, et le couple (taille de nos données, puissance disponible sur la plateform de calcul utilisée).

(1) Li H, Durbin R. Fast and accurate short read alignment with Burrows- Wheeler transform. Bioinformatics 2009;25(14):1754–60.

(2) Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics 2009;25 (16):2078–9.

(3) Ricke, Darrell O.Shcherbina, Anna Michaleas, Adam Fremont-Smith, Philip. 2018GrigoraSNPs: Optimized Analysis of SNPs for DNA Forensics,

(4) Mielczarek, M.Szyda, J.. 2016. Review of alignment and SNP calling algorithms for next-generation sequencing data

(5) https://www.biostars.org/p/788/

(6) https://www.biostars.org/p/18784/

(7) SNP et MNP. https://www.ncbi.nlm.nih.gov/books/NBK44488/

(8) [1] D. Sims, I. Sudbery, N. E. Ilott, A. Heger, and C. P. Ponting, “Sequencing depth and coverage: Key considerations in genomic analyses,” Nat. Rev. Genet., vol. 15, no. 2, pp. 121–132, 2014.