Manuscrit HDR Claire Lemaitre
Manuscrit HDR Claire Lemaitre
Manuscrit HDR Claire Lemaitre
présentée par
Claire Lemaitre
Préambule 4
1 Contexte et problématiques 6
1.1 Les objets biologiques : les variants de structure . . . . . . . . . . . . . . . . . 6
1.1.1 Le génome et les variations génétiques . . . . . . . . . . . . . . . . . . 6
1.1.2 Les différents types de variants de structure . . . . . . . . . . . . . . . 7
1.1.3 Pourquoi étudier les variants de structure ? . . . . . . . . . . . . . . . 7
1.2 Les données : données de séquençage et génomes de référence . . . . . . . . . 9
1.2.1 Le séquençage de l’ADN et ses différentes technologies . . . . . . . . . 9
1.2.2 Quelques exemples d’utilisation en génomique . . . . . . . . . . . . . . 10
1.2.3 Assemblage et utilisation d’un génome de référence . . . . . . . . . . . 11
1.3 Les méthodes : état de l’art des méthodes de détection et d’analyse des va-
riants de structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3.1 La base : l’alignement contre un génome de référence . . . . . . . . . . 13
1.3.2 Différents signaux de mapping utilisés pour la détection . . . . . . . . 13
1.3.3 Un problème difficile avec des lectures courtes . . . . . . . . . . . . . . 15
1.3.4 Plus de 70 logiciels de détection pour les lectures courtes . . . . . . . . 16
1.3.5 De nouvelles méthodes pour les lectures longues . . . . . . . . . . . . . 17
1.3.6 Après la découverte, les autres problèmes méthodologiques associés
aux variants de structure . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.4 Plan du manuscrit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2
3 Vers une meilleure compréhension des faibles performances des outils avec
des lectures courtes 32
3.1 Motivations et contexte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.1.1 Apport des lectures longues pour la détection des variants de structure 32
3.1.2 2019-2020 : les premiers catalogues exhaustifs et précis chez l’homme . 33
3.2 Caractéristiques des vraies insertions chez l’homme . . . . . . . . . . . . . . . 34
3.2.1 Quatre niveaux de caractérisation des insertions . . . . . . . . . . . . . 34
3.2.2 Résultats : la majorité des insertions sont difficiles . . . . . . . . . . . 35
3.3 Identification des causes de la perte de rappel des outils de découverte . . . . 37
3.3.1 Conception d’un benchmark à base de simulations . . . . . . . . . . . 37
3.3.2 Résultats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.3.3 Quelques leçons pour MindTheGap . . . . . . . . . . . . . . . . . . . . 39
3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
5 Discussion et perspectives 51
5.1 Est-ce la fin des lectures courtes pour étudier les variants de structure ? . . . . 51
5.2 Améliorer les méthodes de détection des variants de structure avec diverses
données de séquençage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5.2.1 Perspectives pour l’outil MindTheGap . . . . . . . . . . . . . . . . . . 53
5.2.2 Developper des outils pour les données linked-reads . . . . . . . . . . . 53
5.2.3 Affiner les points de cassure avec des données de lectures longues . . . 55
5.3 Représentation et quantification des Variants de Structure dans les graphes
de génome . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
5.4 Applications et questions biologiques . . . . . . . . . . . . . . . . . . . . . . . 57
Références bibliographiques 65
3
Préambule
Ce document présente une partie de mes travaux de recherche que j’ai effectués au sein
de l’équipe Genscale, équipe de bioinformatique commune au centre Inria Rennes Bretagne
Atlantique et au laboratoire d’informatique IRISA à Rennes. Mes travaux de recherche se
situent à l’interface entre plusieurs disciplines : l’informatique, les mathématiques et la bio-
logie. Cette pluri-disciplinarité est présente dès ma formation universitaire initiale et se
retrouve tout au long de mon parcours professionnel, notamment dans les différents envi-
ronnements de recherche dans lesquels j’ai travaillé : un laboratoire de biologie en thèse,
une plateforme de services de bioinformatique en post-doctorat et un laboratoire d’informa-
tique pour mon poste actuel de Chargée de Recherche. Mon activité consiste à développer et
utiliser des méthodes informatiques pour répondre à des questions biologiques. Depuis mes
premiers travaux en thèse, mes objets biologiques d’étude sont les génomes des organismes
vivants et plus généralement les séquences d’ADN. Les questions méthodologiques portent
sur leur analyse et leur comparaison.
Durant mes doctorat et post-doctorat, entre 2005 et 2010, les données de génomiques
étaient sous la forme de quelques génomes complets, avec une grande qualité de séquence
et d’annotation. Je me suis intéressée à l’évolution des génomes d’organismes divers, des
mammifères aux bactéries, j’ai développé des méthodes permettant de les comparer et j’ai
également analysé les résultats biologiques de ces comparaisons. Mon arrivée, en 2010, en
tant que Chargée de Recherche au centre Inria Rennes Bretagne Atlantique a coïncidé avec
l’apparition de nouvelles technologies de séquençage à haut débit qui ont posé des problèmes
algorithmiques nouveaux pour traiter et comparer une quantité beaucoup plus importante
de données génomiques. Mon activité de recherche s’est naturellement orientée vers ces nou-
veaux problèmes. Les questions biologiques sont similaires, il s’agit toujours d’identifier ce
qui est commun et ce qui diffère entre deux ou plusieurs génomes. Mais les méthodes pour
y répondre doivent s’adapter à deux différences majeures des données : leur qualité et leur
quantité. D’une part, l’information est beaucoup plus morcelée et plus bruitée, ce qui né-
cessite de développer des méthodes spécifiques de pré-traitement et d’analyse des données.
D’autre part, ces données sont beaucoup plus massives et nécessitent de développer de nou-
velles heuristiques plus rapides et/ou de nouvelles structures d’indexation plus légères en
mémoire. Ainsi, depuis 2010, mes travaux de recherche ont porté sur une grande diversité
de problèmes relatifs au traitement et l’analyse des données de séquençage à haut débit :
la correction des erreurs de séquençage, la compression des fichiers de séquençage, l’assem-
blage de génomes, la détection de variations génétiques, ponctuelles ou plus complexes, la
comparaison massive de séquençages métagénomiques, etc.
J’ai choisi de focaliser ce document sur une partie seulement de ces travaux, ceux relatifs
à une problématique biologique qui m’intéresse particulièrement : les variations de structure
dans les génomes. Mon intérêt pour ces variants en particulier remonte à mes tout premiers
travaux de recherche, en thèse. La question qui m’a passionnée durant ma thèse est celle de
4
l’organisation des génomes et l’évolution de leur structure : lorsqu’on compare des génomes
d’espèces proches, on observe un contenu en séquence relativement similaire mais l’ordre et
l’orientation de ces séquences le long des génomes peuvent être très différents. Ces modi-
fications d’organisation génomique sont générées par des réarrangements chromosomiques,
appelés également des variants de structure à l’échelle intra-spécifique. Comprendre les mé-
canismes et les impacts de ces évènements sont des questions biologiques fondamentales
mais qui nécessitent dans un premier temps d’être capable de détecter et caractériser ces
variations dans les génomes. Cette problématique méthodologique est devenu majoritaire
dans mes activités de recherche ces dernières années et constitue mon projet de recherche
pour les années à venir. C’est pourquoi, j’ai fait ce choix de centrer ce document sur cette
thématique. Cependant, mes autres travaux sur les données de séquençages ne sont pas com-
plètement absents de ce document, puisqu’ils se basent notamment sur des concepts et des
outils algorithmiques communs (par exemple, les représentations par ensembles de k-mers
et les graphes de séquences).
Ce document est organisé en 5 chapitres, le premier présente le contexte biologique et
les problématiques méthodologiques de mes travaux, les trois suivants offrent une synthèse
de mes principales contributions pour l’étude des variants de structure avec des données de
séquençage, et le dernier chapitre expose mes perspectives de recherche. Les trois publications
principales sur lesquelles se basent ce document sont ajoutées en annexe du document. Enfin,
un CV présentant les différents éléments nécessaires à l’établissement du dossier d’HDR
clôture ce document.
5
Chapitre 1
Contexte et problématiques
Ce chapitre présente le contexte scientifique dans lequel s’inscrivent les travaux de re-
cherche présentés dans ce document. Nous définissons dans un premier temps les objets
biologiques qui nous intéressent, les génomes et leurs variants de structure. Puis, nous dé-
crivons les données à disposition pour étudier ces objets et qui composent l’entrée de nos
méthodes, avant de présenter les problématiques méthodologiques et un état de l’art des
méthodes pour détecter et étudier les variants de structure dans les génomes.
6
Il existe différents types de variations génétiques. On les distingue principalement par
leur taille. Les plus petites sont les variations dites ponctuelles qui ne touchent qu’une seule
position du génome isolée : ce sont les substitutions d’un nucléotide par un autre (SNP ou
SNV), les délétions d’un nucléotide ou les insertions d’un nucléotide à une position donnée.
Les délétions et insertions d’un ou plusieurs nucléotides successifs sont également appelés
petits indels. Enfin, les plus grandes variations génétiques sont regroupées sous le terme de
Variations Structurales (ou variant de structure) et forment un groupe très hétérogène en
termes de type de mutation et de taille, allant généralement de 50 pb à des bras chromoso-
miques ou chromosomes entiers. En terme de fréquence dans les génomes, plus les variations
sont petites, plus elles sont fréquentes. Ainsi, chez l’homme, sur les 3 milliards de pb qui
composent son génome, on estime le nombre de SNPs de l’ordre du million, dix fois moins
d’indels et 100 fois moins de variants de structure (> 50 pb). Cependant, comme les variants
de structure portent sur des segments d’ADN plus longs, ils touchent globalement une pro-
portion plus importante du génome. Ainsi, on estime que la diversité génomique entre deux
individus humains est en moyenne de 0,1 % lorsqu’on ne considère que les mutations ponc-
tuelles, mais elle augmente à 1,5% si on tient compte des variants de structure (Pang et al.,
2010).
7
Figure 1.1 – Représentation schématique des principaux types de variants de structure, en
haut ceux qui modifient le nombre de copies (variants déséquilibrés), en bas ceux qui ne font
que déplacer des segments d’ADN (variants équilibrés).
8
homologue dans le segment inversé, peuvent générer une barrière reproductive et conduire
à la spéciation des espèces, être associées à de l’adaptation locale, ou encore favoriser la
divergence entre chromosomes sexuels (lire par exemple (Kirkpatrick, 2010) ou encore le
numéro spécial de la revue Molecular Ecology (Wellenreuther et al., 2019)).
Malgré ces nombreux impacts, les variants de structure ont été largement moins étudiés
que les variations ponctuelles et cela est principalement du à des difficultés techniques pour
les détecter dans les données génomiques.
9
Nanopore (ONT). Cette génération se démarque des NGS car elle permet le séquençage
direct de grandes molécules, produisant des lectures de tailles allant de 10 Kb à plusieurs
Mb. Leur principal inconvénient est leur plus fort taux d’erreurs, initialement autour de 20 à
30 % mais qui tombe actuellement en dessous de 13 %, jusqu’à 1 % pour les données PacBio
HiFi. Le type des erreurs est également différent par rapport aux autres générations : il
s’agit principalement d’erreurs d’insertions et de délétions. Enfin, actuellement, leur débit
est plus faible et leur coût plus élevé que la technologie Illumina (voir Table 1 1 de la revue
de Logsdon et al. (2020), pour plus de détails sur les technologies les plus récentes, leurs
caractéristiques et leur coût). Classiquement, on distingue ces deux dernières générations
par les termes lectures courtes versus lectures longues.
Un type intermédiaire de données de séquençage, entre les lectures courtes et les lectures
longues, sont les lectures liées ou lectures longues synthétiques (nous utiliserons le terme
anglais plus commun linked-reads). Ce type de données est apparu en même temps que
la troisième génération de séquençage et a été popularisé par l’entreprise 10X Chromium
Genomics. Ce n’est pas une technologie de séquençage à part entière, l’innovation réside
dans la préparation des librairies des fragments d’ADN avant le processus de séquençage. Le
principe est de marquer des longues molécules, typiquement de taille 50 Kb, par des barcodes
moléculaires (courtes séquences de 10 à 20 pb) avant de les séquencer avec la technologie
Illumina. Le résultat est un ensemble de lectures courtes, dont le début de la séquence
correspond au barcode. Les lectures issues d’une même molécule possèdent le même barcode
et on en tire une information longue distance. En 2020, l’entreprise 10X Chromium Genomics
a arrêté les ventes de ce produit pour des conflits de propriété intellectuelle, mais depuis
2019, d’autres technologies du même type ont été développées, telles que TELL-seq (Chen
et al., 2020), stLFR (Wang et al., 2019) et Haplotagging (Meier et al., 2021).
10
Ce projet compte actuellement plus de 3 000 génomes humains séquencés. D’autres initiatives
centrées sur des populations en particulier ont vu le jour par la suite, comme par exemple
UK10K au Royaume-Unis (10 000 génomes) ou encore GoNL aux Pays Bas (1 000 génomes).
11
complète et contigüe possible. L’annotation du génome consiste ensuite à localiser et décrire
ses éléments fonctionnels tels que les gènes, les séquences régulatrices, les éléments répétés.
Si on prend l’exemple du génome humain, suite à la publication en 2001 de la première
version du génome de référence, de nombreux efforts internationaux ont permis d’améliorer
à la fois sa séquence et son annotation. La version de 2001 ne représentait que la partie
euchromatique (sans les télémoères et centromères), soit environ 92 % du génome, et c’est
seulement en 2021 que l’intégralité des chromosomes a été assemblée sans aucun trou par le
consortium T2T (Telomere-to-telomere) (Nurk et al., 2021). Du côté de l’annotation, environ
20 000 gènes ont été identifiés dans la version de 2001. Cependant, les séquences codantes ne
représentent que 1,5 % des 3 milliards de paires de bases du génome et le projet ENCODE, en
particulier, a permis d’améliorer l’annotation des éléments fonctionnels non codants. Enfin,
il est important de noter que près de la moitié du génome est composé d’éléments répétés,
pour la plupart des éléments transposables répétés en un grand nombre de copies.
Lorsqu’on séquence plusieurs individus d’une même espèce et qu’on possède déjà un
génome de référence pour cette espèce, on parle alors de re-séquençage. L’objectif de tels
séquençages est souvent d’identifier les variations génomiques entre les différents individus
re-séquencés. Pour cela la stratégie usuelle est d’utiliser le génome de référence comme
base de comparaison, avec une première étape de mapping des lectures sur le génome de
référence. Le mapping consiste à identifier d’une part la position génomique d’où provient
une lecture donnée et d’autre part les différences s’il y en a entre la séquence de la lecture et
sa version dans le génome de référence. Cette tâche s’effectue par alignement de séquences, un
autre problème algorithmique fondamental de la bioinformatique. Là encore les algorithmes
développés dans les années 70-80 pour les premières données de séquences ont du être adaptés
pour faire face à la grande quantité de données à traiter. Dans le cas du mapping de lectures
de séquençage sur un génome de référence, la problématique majeure est de réduire le temps
de calcul et la consommation mémoire grâce à des structures de données d’indexation des
séquences.
L’étape suivant l’alignement de toutes les lectures sur le génome de référence s’appelle
l’appel des variants (variant calling en anglais). Elle diffère selon le type de variants re-
cherché. Pour les mutations ponctuelles comme les SNPs et les petits indels, cela consiste
principalement à parcourir le génome de référence et à identifier des positions pour les-
quelles un nombre significatif d’alignements indiquent une même différence. Pour les variants
de structure, nous verrons dans la Section suivante que les stratégies sont plus diverses et
complexes.
12
à une très forte résolution, jusqu’au nucléotide près.
13
Figure 1.2 – Les trois types de signaux de mapping aberrants sont représentés dans le cas
d’une délétion d’un segment (en vert) dans le génome re-séquencé par rapport au génome
de référence. Les lectures de séquençage sont représentées par des flèches, comme elles sont
échantillonnées sur le génome re-séquencé et comme elles sont mappées sur le génome de
référence, et leur appariement est représenté par un arc en pointillé. Les lectures grises sont
concordantes et n’indiquent aucune variation de structure. La paire de lectures rouges est
discordante car la distance de mapping sur le génome de référence est plus grande qu’attendu.
La lecture orange est mappée en deux morceaux de part et d’autre du segment supprimé
(split-read). En bas, la variation de la couverture de séquençage est représentée avec une
chute au niveau du segment vert indiquant l’absence de cette séquence dans le génome
re-séquencé.
Les deux premiers signaux sont caractéristiques de points de cassure, ce sont des ad-
jacences de régions génomiques différentes entre le génome re-séquencé et le génome de
référence. Dans le premier cas, d’une lecture alignée en deux morceaux, le début et la fin de
la lecture sont adjacents dans le génome re-séquencé mais distants ou orientés différemment
dans le génome de référence produisant deux alignements distincts pour ces deux parties
de la lecture. Par exemple, dans le cas d’une délétion d’un segment génomique, une lec-
ture séquencée à cheval sur le point de cassure de la délétion dans le génome re-séquencé,
sera alignée sur le génome de référence en deux alignements de part et d’autre du segment
supprimé.
Le deuxième type de signal suit le même principe mais à l’échelle du fragment de séquen-
çage, représenté par deux lectures appariées (voir Section 1.2.1 p. 9), plutôt que de la lecture
individuelle. La taille des fragments étant contrainte par la technologie, la distance sur le
génome re-séquencé entre deux lectures provenant d’un même fragment suit une distribution
particulière. Ainsi, lorsque deux lectures issues d’un même fragment ne mappent pas sur le
génome de référence de façon à former un fragment d’une taille raisonnable, c’est-à-dire que
leur distance de mapping ne suit pas la distribution attendue ou leur orientation relative
n’est pas compatible avec un fragment séquencé aux deux extrémités (de façon tête-bêche),
on dit qu’elles sont discordantes (voir une exemple dans la Figure 1.2). Cette discordance
de mapping indique là encore une adjacence de séquences dans le génome re-séquencé qui
n’existe pas dans le génome de référence. Si on reprend l’exemple de la délétion, une paire de
lectures avec une lecture de part et d’autre du point de cassure sur le génome re-séquencé,
donnera un mapping discordant car les deux lectures seront alignées de part et d’autre du
14
segment supprimé sur le génome de référence, donc à une distance plus grande qu’attendue.
On distingue plusieurs types de discordances en fonction de la distance et de l’orientation
relative des deux lectures, qui sont générées par les différents types de variants de structure
mais ne leur sont pas toujours spécifiques.
Enfin, le dernier signal n’est généré que par certains types de variants de structure,
les délétions et les duplications. Ces variants modifient le nombre de copies d’une région
génomique donnée. La couverture de séquençage, ou le nombre de lectures mappées couvrant
une position donnée du génome, est un proxy pour estimer ce nombre de copies. Les régions
du génomes ayant subi une diminution ou une augmentation de leur nombre de copies
présentent une couverture de séquençage qui ne suit pas la distribution observée sur le reste
du génome. Dans le cas d’une délétion, à l’état homozygote, le signal est encore plus fort,
puisqu’on observe l’absence de lectures mappées dans cette région.
15
paraître négligeable, et elle l’est dans d’autres contextes, comme la recherche de mutations
ponctuelles dans les portions codantes du génomes (fraction mappable de 97.4 %), mais
dans le cas des variations structurales, même une petite fraction de répétitions peut géné-
rer un grand nombre de fausses détections (faux positifs). Et si on décide de simplement
ne pas regarder ces régions, on réduit fortement la sensibilité de détection, car les variants
de structure ont tendance à être associés aux éléments répétés. D’une part, les répétitions
sont souvent impliquées dans les mécanismes de génération des variants de structure (par
exemple par recombinaison homologue non allélique), et d’autre part une grande partie
des variants de structure sont eux-mêmes des éléments répétés (les duplications, les trans-
positions d’éléments transposables). Plus la lecture est petite, plus la probabilité d’erreur
de mapping ou de mapping ambigu est grande. Cela explique pourquoi l’alignement par
morceaux (split-mapping) est complexe et limité pour les lectures courtes.
16
de logiciels (Kosugi et al., 2019; Cameron et al., 2019). Ces études confirment la supériorité
des dernières méthodes publiées utilisant l’assemblage local. Elles montrent que les outils
ont des performances très différentes en fonction des types de variants de structure, en par-
ticulier les variants de type insertions sont les plus difficiles. Enfin, même si les performances
sur données réelles restent faibles globalement à cause de la petite taille des lectures, elles
suggèrent qu’il reste de la place pour encore améliorer les méthodes.
17
de spécification très dense 3 . Par exemple, la position de fin d’une délétion peut être indiquée
dans un champ libre, soit en indiquant la position de fin sur le génome, soit en indiquant la
taille de la délétion. Une duplication peut être définie par l’un des trois types suivant : DUP,
INS ou CNV. Une inversion peut être décrite par un variant de type INV ou deux variants
de type BND (Break-end, équivalent à un point de cassure).
Il est important également de mentionner le problème de la normalisation de la représen-
tation qui est rencontré lorsqu’il existe des répétitions aussi petites qu’une paire de base au
niveau des points de cassure. Dans ce cas, plusieurs positions associées à différentes séquences
alternatives représentent de manière équivalente le même évènement (un exemple est donné
dans la Figure 1.3). Ce type de répétition s’appelle également de l’homologie jonctionnelle
et n’est pas rare. Ne serait-ce que par hasard, l’alphabet ne comportant que quatre lettres,
ce problème de normalisation est censé se poser pour près de la moitié des variants de type
insertion par exemple.
Figure 1.3 – Exemples des différentes représentations possibles au format VCF d’un variant
de type insertion dont les deux séquences alléliques sont représentées en haut. Cette insertion
possède une répétition de taille 2 (AC, en rouge) aux points de cassure. Dans ce cas, trois
évènements d’insertion sont possibles et aboutissent à la même paire de séquences alléliques.
Leurs représentations au format VCF diffèrent sur la position et les séquences REF et ALT.
La normalisation à gauche (resp. droite) consiste à choisir la représentation avec la position
la plus à gauche (resp. droite).
18
nière intelligente ces comparaisons, tels que SVanalyzer 4 , SURVIVOR 5 et Jasmine (Kirsche
et al., 2021). D’autres problématiques relèvent de la visualisation de ces variants ou de la
quantification des allèles dans de nouveaux échantillons. Ce dernier problème s’appelle le
génotypage et sera décrit plus en détails dans le Chapitre 4.
4. https://github.com/nhansen/SVanalyzer/
5. https://github.com/fritzsedlazeck/SURVIVOR
19
Chapitre 2
Dans ce chapitre, nous présentons une méthode que nous avons développée pour détecter
un type particulier de variant de structure, les insertions, dans un génome re-séquencé en
lectures courtes par rapport à un génome de référence. Ce travail a débuté dans les années
2012-2013, en 2014 la méthode a été publiée dans la revue Bioinformatics (Rizk et al.,
2014) (publication en Annexe A.1) et le logiciel associé, MindTheGap, a été diffusé à la
communauté. Puis, dans les années qui ont suivi et encore à l’heure actuelle, la méthode a
été améliorée et de nouvelles fonctionnalités ont été développées.
20
de référence dans le cas d’une insertion nouvelle, soit apparaissent correctement mappées
à un locus différent dans le cas d’une duplication. La Figure 2.1 montre deux exemples
d’insertions avec des signaux de mapping très différents et/ou similaires à d’autres types de
variants. La caractérisation de la séquence insérée nécessite donc une étape d’assemblage de
novo d’un ensemble de lectures courtes (à identifier).
Figure 2.1 – Types de signaux de mapping aberrants générés par deux types d’insertions
(segment inséré en vert), une insertion nouvelle (A) et une insertion duplicative (B). Contrai-
rement aux délétions qui génèrent les 3 types de signaux de mapping aberrants (Figure 1.2),
on en observe ici beaucoup moins. Pour une insertion nouvelle (A), la majorité des lectures
provenant du segment inséré ne sont pas mappées et il n’y a pas ou très peu de signal dans
la couverture de séquençage. Pour l’insertion duplicative (B), la couverture de séquençage
est modifiée mais elle ne permet pas d’identifier le site d’insertion et les signaux de map-
ping discordants peuvent être interprétés comme provenant d’autres types de variants de
structure (une délétion pour le signal rouge et une transposition pour le signal orange).
21
de calcul, notamment en termes de mémoire pour représenter le graphe d’assemblage. Les
outils génériques détectaient donc très peu de variant de structure de type insertion et les
quelques insertions renvoyées étaient définies uniquement par leur site d’insertion, jamais
avec la séquence insérée. Il existait à notre connaissance que deux outils incorporant un mo-
dule d’assemblage et spécialement dédiés à ce type de variants, SOAPindel (Li et al., 2013)
et NovelSeq (Hajirasouliha et al., 2010) mais qui ne permettaient de détecter chacun qu’un
sous-ensemble bien particulier d’insertions : les grandes insertions nouvelles pour NovelSeq
et les insertions de taille petite à moyenne (max 200 pb) pour SOAPindel. Ces limitations
sont le résultat de l’étape d’assemblage qui se limite à un sous-échantillon des lectures. Pour
SOAPindel, seules les lectures dont le mate est mappé dans la région du site d’insertion sont
utilisées pour l’assemblage, ce qui limite la taille de l’insertion à moins de deux fois la taille
des fragments de séquençage. Pour NovelSeq, seules les lectures non mappées sur le génome
de référence sont assemblées, limitant la détection aux insertions entièrement nouvelles.
Ces choix de sous-échantillonnage étaient en partie dictés par le problème de représenta-
tion en mémoire des graphes d’assemblage construits sur l’ensemble des lectures du génome
re-séquencé. En effet, en 2012, l’assembleur de novo ABYSS (Simpson et al., 2009) nécessi-
tait plus de 300 Go de mémoire pour assembler un génome humain (Chikhi and Rizk, 2013).
Si de telles ressources étaient envisageables pour une tâche ponctuelle d’assemblage d’un
génome de référence, ce n’était pas raisonnable pour la détection de variants dans de nom-
breux génomes re-séquencés. Or, en 2012, Rayan Chikhi et Guillaume Rizk, deux doctorants
de notre équipe, ont proposé de nouvelles méthodes et structures de données permettant de
calculer et représenter en très peu de mémoire le graphe de de Bruijn d’un ensemble de lec-
tures courtes. Le graphe de de Bruijn d’un ensemble de séquences est un graphe dirigé dont
les sommets sont l’ensemble des k-mers, ou mots de taille k, contenus dans les séquences
et les arcs sont les chevauchements exacts suffixe-préfixe de taille k − 1 (Figure 2.2). Si on
élimine les k-mers contenant des erreurs de séquençage et qu’on choisit une valeur de k suf-
fisamment grande (généralement > 30), le génome d’où les lectures ont été échantillonnées
est représenté par un ensemble de chemins dans ce graphe. Rayan Chikhi et Guillaume Rizk
ont proposé d’une part un algorithme de comptage des k-mers sur disque (DSK) qui permet
de filtrer les k-mers avec des erreurs de séquençage sur la base de leur abondance dans le jeu
de lectures (Rizk et al., 2013), et d’autre part une structure de données à base d’un filtre de
Bloom pour représenter en mémoire l’ensemble des sommets du graphe (Chikhi and Rizk,
2013). Ces deux innovations permettent de calculer, représenter et traverser le graphe de
de Bruijn de grands ensembles de lectures courtes avec très peu de mémoire et en temps
raisonnable. Par exemple, l’assemblage d’un génome humain nécessite moins de 6 Go de
mémoire vive avec Minia et peut donc être réalisé sur un ordinateur portable ou de bureau.
Ces travaux ont été ensuite ré-implémentés et diffusés à la communauté dans la librairie
C++ GATB qui facilite le développement d’outils basés sur cette implémentation du graphe
de de Bruijn (Drezen et al., 2014).
Ces innovations et leurs implémentations ont démocratisé l’utilisation du graphe de de
Bruijn pour bien d’autres tâches que celle à laquelle elle était cantonnée jusqu’à présent,
l’assemblage de novo d’un génome de référence. À partir de cette implémentation, nous
avons notamment développé un logiciel de correction de lectures courtes (Bloocoo) (Benoit
et al., 2014), un compresseur de fichier de séquençage (Leon) (Benoit et al., 2015), un outil
de découverte de petits variants sans génome de référence (discoSnp) (Uricaru et al., 2015),
et, ce qui nous occupe dans ce manuscrit, une méthode d’assemblage local pour la détection
des insertions (MindTheGap).
22
Figure 2.2 – Exemple d’un graphe de de Bruijn pour k=4, construit à partir de 3 lectures
(en haut). Après avoir décomposé les lectures en k-mers, ces derniers sont comptés (étape
non représentée ici) pour éliminer les k-mers rares. L’ensemble des k-mers distincts forment
l’ensemble des sommets du graphe, les arcs correspondent aux chevauchements exacts suffixe-
préfixe de taille exactement k −1. Utilisé pour l’assemblage de novo de séquences, on cherche
ensuite des chemins avec des propriétés particulières, en bas est représentée la séquence
assemblée issue du chemin le plus long.
23
Figure 2.3 – Les différentes étapes de la méthode MindTheGap.
En effet, contrairement à la quasi totalité des outils de détection des variants de structure
qui prennent en entrée le résultat d’une étape de mapping des lectures sur le génome de
référence, MindTheGap en est indépendant. La détection des sites d’insertion ou points de
cassure se fait là encore en exploitant cette représentation efficace des données par leur
graphe de de Bruijn, en comparant celui-ci avec le génome de référence, sur la base de leurs
k-mers communs et spécifiques.
Plus précisément, on recherche le long du génome de référence des motifs particuliers
d’absence de k-mers (insertions homozygotes) ou de k-mers ayant plusieurs voisins (insertions
hétérozygotes) qui sont générés par les insertions. Dans le cas d’une insertion homozygote,
la séquence du génome de référence au site d’insertion est absente dans le génome rere-
séquenséquencé. En termes de présence-absence de k-mers, cela se traduit dans la plupart des
cas par k −1 k-mers consécutifs sur le génome de référence qui sont absents dans le graphe de
de Bruijn du génome re-séquencé (Figure 2.3 à gauche dans l’encart du module Find). Nous
appelons ces stretchs de k-mers absents des gaps. D’autres variants homozygotes génèrent
également des gaps, mais pour la plupart leur taille est différente : elle est d’exactement k
pour un SNP et strictement supérieure à k pour une délétion. D’autres variants, comme les
inversions ou translocations, génèrent des gaps de mêmes caractéristiques que les insertions,
mais ils sont beaucoup moins abondants dans les génomes que les SNPs et les délétions.
Ainsi, ce motif de gap de taille k − 1 est relativement spécifique à des variants d’insertions et
est surtout très facile à énumérer en parcourant les k-mers du génome de référence. Il faut
noter cependant que la présence d’une séquence répétée entre une extrémité de la séquence
insérée et les séquences flanquant le site d’insertion peut altérer ce motif, réduisant la taille
du gap de la taille de la répétition (voir un exemple dans la Figure 1 du papier en Annexe
A.1). Un paramètre permet à l’utilisateur de prendre plus ou moins en compte ces répétitions
potentielles, qui sont appelées parfois micro-homologies ou homologies jonctionnelles et nous
verrons dans le chapitre 3 qu’il a de l’importance pour certaines données.
Dans le cas d’une insertion hétérozygote, le site de l’insertion est présent dans un des
haplotypes du génome re-séquencé, on n’observe donc pas d’absences de k-mers, mais on
24
observe un motif particulier en termes de k-mers branchants (des k-mers ayant plusieurs
voisins dans le graphe) : un k-mer branchant à droite suivi à une distance de k par un k-mer
branchant à gauche (Figure 2.3 à droite dans l’encart du module Find). De la même manière
que pour les gaps, c’est la distance de k entre les deux k-mers branchants qui différencie le
type de variant entre une insertion, un SNP ou une délétion, et qui peut être réduite en cas
de répétition entre la séquence insérée et le site d’insertion.
Ce choix d’éviter une étape de mapping a été guidé par plusieurs éléments : i) le gain
en temps de calcul puisque le mapping est une tâche longue et que le graphe de de Bruijn
est construit une seule fois pour les deux étapes, ii) contrairement aux signaux de mapping
anormaux, un unique motif (voire deux en fonction de la zygotie du variant) est recherché
quelque soit le type, la taille des insertions et les caractéristiques de taille de fragments des
données, iii) le besoin pour l’étape suivante de récupérer des k-mers voisins du site d’insertion
effectivement présents dans le graphe de de Bruijn et dont la position vis-à-vis de l’insertion
est maîtrisée.
Lors de l’étape Fill, le graphe de de Bruijn est parcouru en partant du k-mer à gauche de
chaque site d’insertion par un parcours en largeur. Un ensemble de contigs est construit avec
l’algorithme de l’assembleur de novo, Minia (Chikhi and Rizk, 2013), qui écrase les petites
bulles ou tips, jusqu’à atteindre des paramètres limites d’exploration locale du graphe en
termes de nombre de contigs construits ou nombre total de paires de bases assemblées. Puis
si le k-mer de droite est retrouvé dans un ou plusieurs de ces contigs, l’ensemble des chemins
reliant ces deux k-mers sont énumérés. Une autre originalité de l’approche est de renvoyer
plusieurs solutions si elles existent, plutôt que de s’arrêter à la première solution. Cependant
cela entraine une difficulté supplémentaire, car des petits polymorphismes ponctuels peuvent
générer un grand nombre de solutions peu différentes. Afin de limiter la redondance dans
les résultats due à du polymorphisme ponctuel, les différentes solutions sont alignées et les
plus similaires sont représentées par une unique séquence pour renvoyer un ensemble de
séquences présentant deux à deux au moins 10 % de divergence.
Pour les détails algorithmiques et d’implémentation nous renvoyons le lecteur vers la
publication de la méthode donnée en Annexe A.1.
25
dans le génome re-séquencé et dans le génome de référence. Cela empêche la détection de
sites d’insertions localisés à proximité (moins de k nucléotides) d’un autre variant, et en
particulier de variants souvent simples à détecter par ailleurs et très fréquents : les SNPs.
Nous avons donc étendu les définitions des motifs recherchés et implémenté des algorithmes
permettant de détecter des SNPs et délétions homozygotes, à proximité ou isolés des sites
d’insertions. Nous avons vu dans la section précédente que les variants homozygotes génèrent
des stretchs de k-mers le long du génome de référence qui sont absents dans le génome re-
séquencé (des gaps), dont la taille dépend du type de variant. Un site d’insertion isolé génère
un gap de taille inférieure à k − 1. Les SNPs isolés génèrent des gaps de taille exactement
k et les délétions des gaps de taille supérieure à k, dépendant principalement de la taille du
segment supprimé. Lorsque ces variants sont proches (à moins de k nucléotides), les motifs
se chevauchent sur le génome et génèrent un plus grand gap. Après avoir identifié ces trop
grands gaps, l’idée de la méthode est d’essayer les différentes façon de corriger le gap de
chaque côté par la recherche de k-mers proches dans le graphe de de Bruijn (ou micro-
assemblage). Cette méthode de correction des gaps dans le génome est notamment inspirée
d’une méthode de correction de novo des erreurs de séquençage dans les lectures courtes,
que nous avions proposée auparavant, appelée Bloocoo (Benoit et al., 2014).
Ces développements permettent de détecter d’autres types de variants (délétions et
SNPs) en plus des insertions, avec comme objectif de résoudre le maximum de gaps ou
différences de k-mers entre le génome re-séquencé et le génome de référence. Mais surtout,
cela a permis d’améliorer le rappel de la détection des variants de type insertion en présence
d’autres polymorphismes dans les données de séquençage.
26
le module Fill de MindTheGap intervient en deuxième étape de la méthode MinYS, pour
reconstruire les séquences manquantes ou trop divergentes du génome de référence.
La première étape construit un ensemble de contigs à partir des lectures alignées sur
le génome de référence. Dans la deuxième étape, les extrémités de ces contigs sont donnés
en entrée de MindTheGap, qui effectue des assemblages locaux avec le graphe de de Bruijn
construit à partir de la totalité des lectures de séquençage métagénomiques, et non plus les
seules lectures alignées sur le génome de référence. Afin d’être le moins possible influencé par
la structure du génome de référence, on effectue une recherche exhaustive entre toutes les
paires possibles d’extrémités de contigs, sans tenir compte de l’ordre ni de l’orientation de
ces contigs dans le génome de référence. C’est sur ce point, que de nouveaux développements
sur le module Fill ont été apportés. L’utilisation naïve de ce module impliquerait d’effectuer
2×N ×2×(N −1) explorations locales du graphe de de Bruijn, si N est le nombre de contigs
initial. En effet, pour chaque extrémité de contig (source), il y a 2 × (N − 1) cibles possibles.
L’algorithme a été modifié pour permettre de rechercher à partir d’une séquence source tous
les chemins possibles vers un ensemble de séquences cibles (et non plus une seule), en une
seule exploration locale du graphe, tout en évitant de trouver des solutions redondantes qui
ré-assembleraient, par exemple, des contigs initiaux.
Dans cette approche, nous avons ainsi tiré partie de deux spécificités de MindTheGap : la
capacité à passer à l’échelle de gros jeux de données (avec la structure de données compacte
du graphe de de Bruijn) et le fait que plusieurs solutions sont renvoyées si elles existent. Ce
dernier point permet notamment d’identifier des variants de structure ou des souches avec
des structures de génome différentes qui co-existent dans le microbiote séquencé. Sur ces deux
points MinYS se démarque nettement de l’état de l’art. Ce travail a fait l’objet d’une partie
de la thèse de Cervin Guyomar que j’ai co-encadré et d’un papier publié récemment dans la
revue NAR Genomics and Bioinformatics (Guyomar et al., 2020). Le logiciel, implémenté
en Python, est diffusé sur github 1 et comme un paquet Bioconda 2 .
27
pour des insertions de 1 Kb dans un chromosome humain). Les sites d’insertion sont bien
détectés quelque soit la taille de l’insertion, mais c’est bien l’étape d’assemblage local qui
est impactée par la présence de régions dans le graphe de de Bruijn plus complexes, très
denses en branchements, fréquemment générées par des séquences répétées en un grand
nombre de copies comme les éléments transposables. Dans ce cas, les limites d’exploration
locale du graphe sont atteintes avant d’avoir atteint la séquence cible. Les performances
sont également un peu moins bonnes pour les insertions qui sont à l’état hétérozygote dans
les données de séquençage par rapport aux insertions homozygotes (perte de 10 à 30 % de
sensibilité). Cela était attendu car, d’une part, le motif de k-mers généré par les insertions
hétérozygotes est moins spécifique que celui des insertions homozygotes (il peut notamment
être confondu avec des régions répétées), et d’autre part, la profondeur de séquençage est
deux fois plus faible pour assembler la séquence insérée.
Ces résultats obtenus sur des données simulées assez simples ont constitué une preuve de
concept que l’approche fonctionnait. Les tests sur données simulées ont permis d’analyser
finement les raisons des erreurs et/ou manques de détection (faux positifs et faux négatifs),
d’analyser l’influence des différents paramètres sur les résultats et d’optimiser la méthode et
ses paramètres. Les données simulées ont également permis de comparer notre approche à
d’autres logiciels. Cependant, au moment de la publication, seuls deux logiciels permettaient
d’assembler des variants d’insertion, SOAPindel (Li et al., 2013) et NovelSeq (Hajirasouliha
et al., 2010). Nous n’avons pas réussi à faire fonctionner Novelseq, qui n’est plus maintenu
depuis 2012 et SOAPindel, comme attendu, s’est révélé limité aux insertions de petites
tailles. Ainsi, cette comparaison a essentiellement mis en évidence que MindTheGap était le
seul outil en 2014 à pouvoir détecter et assembler des variants d’insertions plus grands que
la taille des fragments séquencés.
Si les données simulées présentent de nombreux avantages pour l’évaluation et la com-
paraison des méthodes, elles ont cependant une limite majeure de ne pas toujours bien
représenter la réalité biologique et de sous-estimer la difficulté du problème. Dans ce cas par-
ticulier, nos données simulées présentent deux limitations importantes : l’absence d’autres
polymorphismes dans les données de séquençage (SNPs et petits indels qui sont beaucoup
plus fréquents en pratique que les grandes insertions), et le fait que les insertions ont été
simulées à des localisations génomiques aléatoires (probabilité uniforme le long du génome),
ce qui ne reflète pas la réelle complexité des insertions et de leurs points de cassure. Or, en
2014, il existait très peu de données avec des longs variants d’insertions dont la séquence
insérée était résolue et validée (voir Chapitre 3, page 32). Ainsi, la validation de MindThe-
Gap sur des données réelles n’a été effectuée que sur une vingtaine d’insertions identifiées
au préalable chez un individu humain par du séquençage de fosmids. Cette expérimenta-
tion a notamment permis de montrer que MindTheGap était capable d’assembler avec une
grande précision certaines de ces insertions, notamment de très longues (> 4 Kb), avec des
ressources de calcul très raisonnables pour un jeu de données contenant plusieurs milliards
de lectures. Elle a également permis de mettre en évidence plusieurs pistes d’amélioration,
comme la gestion d’autres polymorphismes (SNPs et délétions) à proximité des points de
cassure qui a été résolue par la suite (voir Section 2.2.2).
28
étudiés sont très divers, allant de génomes viraux, de mitochondries, de bactéries à des
génomes plus complexes d’insectes, de plantes et des exomes humains. Ces résultats ont été
obtenus notamment lors de collaborations avec des biologistes, mais également sans notre
collaboration suite à la publication du logiciel comme l’attestent au moins cinq publications
(Burioli et al., 2017; Carrier et al., 2018; Daval et al., 2019; Feldman et al., 2019; Fuentes
et al., 2019). Parmi ces publications, on peut notamment citer celle de Fuentes et al. (2019)
parue dans la revue à fort impact Genome Research en 2019, qui étudie le polymorphisme
de structure de 3 000 génomes du riz, dans laquelle MindTheGap a été l’outil choisi pour
détecter les grandes insertions.
Dans le cadre d’une collaboration à long terme avec l’INRAE sur le génome du puceron
du pois, nous avons utilisé MindTheGap à plusieurs reprises : (i) pour rechercher des variants
d’insertion qui pourraient expliquer la variabilité du nombre de lectures non mappées sur le
génome de référence entre différentes populations et races d’hôtes de puceron (Gouin et al.,
2015), (ii) pour obtenir le génome de référence d’un symbiote encore peu caractérisé du
puceron du pois, Rickettsia sp (Guyomar et al., 2018) et (iii) pour répertorier les variants
structuraux d’un phage responsable de traits phénotypiques importants (phage APSE de
la bactérie Hamiltonella defensa, publication en cours de préparation). Dans le cadre de
collaborations encore en cours sur des génomes de papillons, nous utilisons MindTheGap
comme gap-filler pour l’assemblage de génomes mitochondriaux ou de régions d’intérêt.
C’est le cas par exemple, pour le locus Supergene chez le papillon mimétique Heliconius
numata, qui présente un polymorphisme d’inversion associé au patron de coloration des
ailes des papillons. Nous utilisons le pipeline MTG-link pour reconstruire cette région de
1,3 Mb dans une douzaine d’individus re-séquencés avec des lectures liées et pouvoir ainsi
étudier plus précisément leurs différences structurales.
Enfin, plus récemment, dans le cadre de la thèse de Wesley Delage, que j’ai co-encadré
avec Julien Thévenon, médecin-généticien spécialiste des maladies rares, MindTheGap a été
utilisé pour des applications de santé humaine, en collaboration avec l’Inserm de Grenoble
et le Centre Hospitalier Universitaire de Dijon. MindTheGap a notamment été appliqué sur
des données de séquençage d’exomes de patients atteints de maladies rares, pour lesquelles
certaines variations génomiques impliquées sont connues. Les essais réalisés sur une trentaine
de séquençages d’exomes n’ont pas permis de retrouver les variations attendues, mais ont
produit des ensembles de prédictions qui nécessitent d’être analysés plus en détails. En effet,
les quelques variations connues à l’heure actuelle étaient des CNV dont les localisations ne
sont pas précisemment connues, ce qui rend difficile une comparaison avec MindTheGap,
et les grandes insertions (> 50 pb) qui sont la cible principale de MindTheGap sont peu
fréquentes dans les exons et encore peu caractérisées dans le diagnostic clinique car aucune
méthode d’analyse standardisée n’existe. Cependant, une réussite de cette expérience a été
l’utilisation de MindTheGap sur des données réelles par des utilisateurs du domaine médical
et d’avoir des retours sur l’outil.
29
3 humain, MindTheGap met environ 15 minutes à détecter et assembler les 200 insertions
simulées. Sur des vraies données de séquençage du génome de C. elegans, on obtient un
temps de calcul très similaire. Mais les temps de calcul observés sur un génome humain
étaient d’un tout autre ordre de grandeur et montraient une augmentation non linéaire avec
la taille des données en entrée. En effet, au lieu d’une multiplication par 10 attendue du
fait de la taille du génome, soit quelques heures, on obtenait des temps déraisonnables de
plusieurs centaines d’heures, dus essentiellement au module Fill. Sans faire une analyse très
détaillée de la complexité des algorithmes, on peut facilement estimer que la complexité du
module Find de détection des sites d’insertion est linéaire avec la taille du génome, et celle
du module Fill est linéaire avec le nombre de sites d’insertions détectés à l’étape du Find.
On aurait pu penser que ce nombre de sites d’insertion augmenterait de manière linéaire
avec la taille du génome. Cela est sûrement vrai pour une séquence aléatoire ou un génome
simple. Mais dans le cas du génome humain, qui est un génome complexe, représenté à
plus de la moitié par des séquences répétées, on observe que ce nombre de sites n’augmente
pas linéairement : la proportion de sites faux positifs est plus importante dans un génome
complet qu’avec un seul chromosome. De plus, le temps de calcul moyen pour assembler
chaque site est plus de 20 fois plus grand avec le graphe construit sur le génome complet
qu’avec un seul chromosome (par exemple, près de 20 secondes versus moins d’1 seconde
par site), alors que les paramètres de limite d’exploration du graphe ne dépendent pas de la
taille de celui-ci.
Ainsi, la taille et la complexité du génome humain font augmenter le nombre de sites à
assembler et le temps moyen pour les assembler. Mais, il existe encore un autre facteur qui
aggrave encore le temps de calcul, et qui est souvent négligé dans les données simulées, c’est le
vrai polymorphisme contenu dans les données réelles de séquençage. La très grande majorité
de ce polymorphisme est sous la forme de SNPs et de petits indels (1 à 2 pb). En particulier,
ces petits indels (les insertions, mais aussi une partie des petites délétions) génèrent le même
motif de k-mers que les grandes insertions visées par MindTheGap. Ce qui semblait être un
avantage de MindTheGap de ne pas dépendre de la taille des insertions, devient alors un
énorme fardeau. Ainsi en appliquant le module Find à un vrai séquençage plein génome du
génome humain, on peut obtenir plus de 250 000 sites potentiels d’insertions à assembler,
soit une estimation du temps d’assemblage à plus d’un mois (sans parallélisation) !
Nous avons proposé plusieurs solutions, qui sans résoudre complètement ce problème sur
données réelles humaines, en diminuent l’impact (division par 2,5 voire 3 du nombre de sites
à assembler) :
— l’imputation des SNPs homozygotes de l’individu re-séquencé dans le génome de ré-
férence avant de lancer le module Find. Cette étape rajoute souvent peu de temps
de calcul, puisque dans de nombreuses analyses, les SNPs sont les premiers poly-
morphismes à être détectés et étudiés et ce type de données est donc très souvent
disponible.
— l’assemblage dès le module Find des petites insertions. L’approche est similaire à celle
utilisée pour les variants proches, des tentatives d’assemblage local de 1 à 2 pb sont
effectuées pour chaque site potentiel d’insertion. Cela évite d’explorer dans le module
Fill une partie du graphe beaucoup plus importante que nécessaire pour ces petites
insertions qui sont généralement triviales à assembler.
— un filtre sur les caractéristiques de branchement des k-mers voisins au site d’insertion
est appliqué pour limiter le nombre de sites faux positifs. Les sites faux positifs sont
le plus souvent détectés avec le motif hétérozygote, dans des régions répétées qui se
caractérisent par une grande densité de k-mers branchants.
30
2.4 Conclusion
Nous avons présenté dans ce chapitre une méthode originale pour détecter et assembler
des variants d’insertions avec des données de lectures courtes. Cette méthode était l’une
des premières quand elle a été publiée à faire intervenir une structure de données et des
algorithmes d’assemblage de séquence pour la détection de variants de structure. Cette
stratégie a montré son efficacité pour les variants de type insertion, mais elle est maintenant
de plus en plus courante dans la nouvelle génération d’outils de détection génériques pour
tous les types de variants de structure (Chen et al., 2016; Cameron et al., 2017; Wala et al.,
2018).
Le développement du logiciel MindTheGap a débuté dans les années 2012-2013, et a
impliqué, à des degrés d’implications variables, plusieurs personnes que j’ai encadrées ou
avec qui j’ai collaboré (Guillaume Rizk, Anaïs Gouin, Pierre Marijon, Cervin Guyomar, et
Wesley Delage). Il se poursuit encore à l’heure actuelle pour améliorer les performances de
la méthode ou ajouter de nouvelles fonctionnalités. Le logiciel, implémenté en C++, est
diffusé librement sur github 3 et comme un paquet bioconda 4 . Il dispose d’une communauté
d’utilisateurs comme en témoignent le nombre de téléchargements 5 et les retours d’utili-
sateurs que nous recevons. Ainsi, il a été appliqué avec succès pour produire des résultats
biologiques, principalement sur des petits génomes ou des exomes. Son intégration dans des
pipelines d’assemblage ou d’amélioration d’assemblage de génomes est également un succès
récent (outils MTG-link et MinYS (Guyomar et al., 2020)).
Un domaine d’application que nous avons commencé à étudier ces trois dernières années
mais qui mériterait d’être plus développé est celui de la santé humaine. Dans ce contexte,
nous avons rencontré deux difficultés principales : le passage à l’échelle de la méthode sur
des données plein génome du génome humain et l’évaluation des résultats obtenus dans un
contexte où peu de grandes insertions avec un intérêt clinique ont été caractérisées.
Enfin, nous verrons dans le chapitre suivant que l’arrivée des données de séquençage de
troisième génération nous permet d’évaluer différemment ces outils et remet en question
certains résultats et choix méthodologiques.
3. https://github.com/GATB/MindTheGap
4. https://anaconda.org/bioconda/mindthegap
5. Chiffres donnés par Bioconda : près de 4 000 sur les 2 dernières années.
31
Chapitre 3
32
nombreux variants de structure peuvent être entièrement inclus dans une même lecture
et que le mapping des lectures ou morceaux de lectures (split-mapping) sur le génome de
référence est beaucoup plus spécifique. Ainsi, en 2017, les premiers logiciels de découverte
des variants de structure avec des lectures longues ont été développés (voir Section 1.3.5
p. 17), mettant en évidence de nombreux nouveaux variants de structure par rapport aux
ensembles détectés seulement avec des lectures courtes (Sedlazeck et al., 2018b; Mahmoud
et al., 2019). L’assemblage de novo de génomes est également facilité avec ces technologies
et la stratégie de découverte de variants de structure par comparaison d’assemblages est
devenu possible, même si elle reste plus coûteuse en ressources de calcul que l’approche par
mapping.
33
sur ces jeux de données. Les résultats ont produit un choc et une grande déception : Mind-
TheGap n’était capable de détecter qu’une infime fraction de ces variants réels, moins de 5
%. Notre outil n’était pas le seul en cause, d’autres outils concurrents testés n’atteignaient
pas plus de 10 % et les publications du HGSV et de GiaB avaient montré le faible rappel
des détecteurs de variants de structure basés sur les lectures courtes et en particulier pour
les variants de type insertion. Cependant, l’écart de rappel avec nos résultats sur données
simulées (entre 75 et 95 % contre moins de 5%) était frappant et témoignait d’une faible
connaissance des caractéristiques des vraies insertions.
C’est pourquoi, dans ce travail, nous nous sommes attachés à caractériser le plus finement
possible l’ensemble des variants d’insertions catalogués par ces deux études. Puis, nous avons
cherché à identifier parmi les caractéristiques répertoriées, lesquelles étaient liées à une perte
de rappel des outils de détection d’insertions basés sur des lectures courtes. Ces analyses
ont été publiées dans le journal BMC genomics (Delage et al., 2020), nous présentons les
principaux résultats obtenus dans les deux sections suivantes et renvoyons la lectrice et le
lecteur à la publication produite en Annexe A.2 pour les détails des résultats, les figures et
les méthodes.
Annotation du type d’insertion Nous avons vu qu’il existait deux grands types : les
insertions nouvelles dont la séquence est absente du génome de référence et les insertions
duplicatives (ou duplications) dont la séquence existe dans le génome de référence. Parmi les
duplications, on peut distinguer des sous-types très différents, en fonction de la localisation
et du nombre de copies de la séquence dupliquée dans le génome. Par exemple, on distingue
l’insertion d’un élément transposable, qui possède plusieurs centaines voire milliers de copies
très similaires dans le génome, d’une grande séquence dupliquée en tandem, ou encore de
l’expansion d’un motif court répété en tandem un grand nombre de fois dans la séquence
insérée. On peut s’attendre à ce que ces différents types génèrent des signaux différents de
mapping des lectures et des difficultés différentes pour l’assemblage. Ainsi, pour analyser
34
les causes de perte de rappel des outils de détection, il est important de savoir annoter et
distinguer ces différents types. Ce travail d’annotation n’avait été fait que partiellement et de
manière hétérogène entre les deux catalogues de variants de structure du HGSV et du GiaB,
ne permettant pas de comparer les deux études, ni de représenter tous les types existants.
Nous avons alors proposé une annotation en cinq types, qui permet de couvrir l’ensemble
des différents types et sous-types de la base de données dbVar 3 :
— insertion de novo : absence de la séquence inserée dans le génome de référence ;
— insertion d’un élément mobile : séquence similaire à des éléments mobiles connus ;
— expansion d’une répétition en tandem : séquence composée d’une graine répétée en
tandem ;
— duplication en tandem : séquence dont une copie est présente dans le génome de
référence aux abords du site d’insertion ;
— duplication dispersée : séquence dont une copie est présente dans le génome de réfé-
rence à un locus différent du site d’insertion.
Nous avons alors proposé une méthodologie d’annotation qui permet d’assigner un type
de manière non ambigüe à une grande majorité des variants d’insertion. Cette méthodo-
logie est basée principalement sur l’alignement de la séquence insérée contre le génome de
référence, les régions adjacentes au site d’insertion et une banque d’éléments transposables.
Cependant, dans de très nombreux cas, les alignements obtenus ne couvrent pas entière-
ment la séquence insérée ou bien peuvent relever de différents types. Nous avons donc fait
intervenir un seuil de couverture de la séquence insérée qui est paramétrable et un arbre de
décision pour annoter les insertions qui pourraient relever de plusieurs types (voir la Figure
1 de l’article en Annexe A.2).
35
Chaque catalogue contient environ 14 000 insertions, avec plus de la moitié des insertions
d’un catalogue qui lui sont spécifiques. Malgré le faible nombre d’insertions partagées entre
catalogues, les distributions des différentes caractéristiques des insertions sont très similaires
d’un individu à l’autre. En particulier, la répartition dans les différents types d’insertion est
stable même entre les deux études qui ont utilisé des méthodologies différentes pour consti-
tuer leur catalogue. Les distributions des caractéristiques des insertions sont représentées
dans la Figure 3.1.
Figure 3.1 – Distributions des caractéristiques des variants d’insertion dans plusieurs ca-
talogues humains. Distributions (a) de la taille des insertions, (b) du type d’insertion, (c)
du contexte génomique d’insertion et (d) de la taille des homologies jonctionnelles. Abré-
viations : SimpleRep pour répétition simple, ME ou TE pour élément mobile/transposable,
TandemRep pour répétition en tandem, TandemDup pour duplication en tandem, Dispers-
Dup pour duplication dispersée. Figure extraite de (Delage et al., 2020).
Le point remarquable de cette analyse est que la grande majorité des insertions peuvent
être qualifiées de difficiles, c’est-à-dire qu’elles présentent des caractéristiques qui sont sus-
ceptibles de rendre leur détection difficile avec des lectures courtes. En quelques chiffres clés,
36
on observe que 63 % des insertions sont des expansions d’un court motif répété en tandem
(répétition en tandem), le deuxième type le plus représenté étant les insertions d’éléments
mobiles (16%). En terme de localisation génomique, la répartition des insertions n’est pas
uniforme le long du génome. On observe un fort biais vers les régions répétées et difficilement
mappables du génome, avec plus de 70% des insertions qui sont localisées dans des régions
répétées dites simples (court motif répété en tandem) alors que ces régions ne représentent
que 3 % du génome. Enfin, concernant la complexité des points de cassure, plus de 40 % des
insertions possèdent une homologie jonctionnelle de taille supérieure à 10 pb. Là encore, c’est
plus qu’attendu par chance, sur 2 000 insertions simulées uniformément le long du génome
et avec une séquence aléatoire, la taille médiane d’homologie jonctionnelle mesurée est de 0
et la plus grande observée est de 7 pb.
Les différents niveaux de caractérisation, c’est-à-dire le type, la taille, le contexte gé-
nomique et l’homologie jonctionnelle, ne sont pas indépendants les uns des autres. Ainsi,
on observe des distributions de taille, de localisation génomique et d’homologie jonction-
nelle différentes en fonction du type d’insertion. Par exemple, la très grande majorité des
insertions de type répétition en tandem sont localisées dans des répétitions en tandem. Ces
insertions en particulier cumulent donc plusieurs niveaux de difficulté pour les outils de dé-
tection avec des lectures courtes et on peut se demander si la perte de rappel est plutôt due
au type de la séquence insérée ou bien à son contexte génomique d’insertion.
37
des lectures de 2x250 pb sur le chromosome 3 humain qui est altéré par 200 insertions
dont les caractéristiques de taille, de type de séquence insérée, de localisation génomique et
d’homologie jonctionnelle au point de cassure sont contrôlées.
Un premier jeu de données est appelé référence et correspond aux insertions les plus
faciles théoriquement à détecter : des insertions de taille intermédiaire (250 pb, la taille
des lectures), de type insertion nouvelle (la séquence insérée est absente du chromosome 3
humain, elle contient pas ou peu de répétitions puisqu’elle est tirée d’une exon d’un autre
génome éloigné, Sacharomyces cerevisiae), localisée dans un contexte génomique non répété
(les exons) et il n’y a pas d’homologie jonctionnelle au site d’insertion. Ensuite, nous avons
décliné quatre scénarios de simulation où un seul des quatre niveaux de caractérisation est
modifié par rapport à la simulation de référence :
— scénario 1 : la taille des insertions. Il contient trois jeux de données simulées avec
trois tailles différentes d’insertion : 50, 500 et 1 000 pb. Pour chacun de ces jeux
de données, les séquences insérées de la simulation de référence ont simplement été
agrandies ou rétrécies sans changer leur localisation.
— scénario 2 : le type d’insertion. Cinq jeux de données ont été simulés avec les types
suivants : duplication dispersée, duplication en tandem, insertion d’éléments mobiles
(type SINE), répétitions en tandem avec deux tailles de motif (6 et 25 pb).
— scénario 3 : l’homologie jonctionnelle, avec cinq tailles testées de 10 à 150 pb. Pour
simuler l’homologie jonctionnelle d’une taille X donnée, le premier X-mer de chaque
séquence insérée de la simulation de référence est remplacé par le X-mer flanquant à
droite le site d’insertion correspondant.
— scénario 4 : le contexte génomique d’insertion. En se basant sur les annotations de
RepeatMasker des séquences répétées du chromosome 3, nous avons distingué cinq
contextes en plus de celui de la simulation de référence (les exons) : les régions
sans répétition, les répétitions simples courtes (<300 pb) et longues (>300 pb), les
éléments mobiles de type SINEs et LINEs.
Sur cette vingtaine de jeux de données simulées, nous avons appliqué quatre logiciels
permettant de détecter des insertions. Nous avons choisi trois logiciels de détection de va-
riants de structure génériques, parmi les plus récents et les plus utilisés par la communauté,
qui détectent tous types de variants de structure et pas seulement les insertions, mais qui
sont théoriquement les mieux adaptés aux variants d’insertions dans cette catégorie puis-
qu’ils utilisent tous des techniques d’assemblage local. Il s’agit de Manta (Chen et al., 2016),
SVaba (Wala et al., 2018) et GRIDSS (Cameron et al., 2017). Le quatrième outil testé
est naturellement l’outil que nous avons développé qui lui est dédié à ce type de variant,
MindTheGap (voir Chapitre 2). Pour chaque outil et sur chaque jeu de données, nous avons
mesuré le nombre d’insertions correctement prédites par rapport à la vérité simulée, soit le
rappel. Il faut noter que nous avons distingué deux types de rappel en fonction de la préci-
sion des prédictions : le rappel du site d’insertion où seule la localisation de l’insertion est
évaluée et le rappel avec séquence résolue où pour être considérée comme un vrai positif une
prédiction doit avoir une séquence assemblée correctement (avec au moins 90 % d’identité
de séquence avec la séquence simulée).
3.3.2 Résultats
Les résultats détaillés de ce benchmark sont présentés dans les Tableaux 1 et 2 de l’article
en Annexe A.2, le deuxième tableau est repris de manière simplifiée ci-après dans la Table
3.1. Nous résumons ici quelques points remarquables. Tout d’abord, comme attendu, les ou-
38
tils n’ont aucun problème, avec des rappels de 100 %, si les insertions simulées sont simples
comme dans la simulation de référence (de taille intermédiaire (250 pb) sans répétition dans
la séquence insérée, ni dans le contexte génomique, ni au point de cassure). Cependant,
dès qu’on introduit des caractéristiques plus difficiles, les rappels des outils peuvent chuter
jusqu’à moins de 5 % et les valeurs sont extrêmement variables en fonction de la caracté-
ristique simulée mais également de l’outil en question. Parmi les caractéristiques qui ont le
plus d’impact, on notera que le type de la séquence insérée est plus problématique que le
contexte génomique dans lequel elle s’insère, même lorsqu’on ne s’intéresse qu’au site d’in-
sertion. Ainsi, pour répondre à la question donnée en exemple, c’est le caractère de motif
répété en tandem présent dans la séquence à assembler qui semble limiter le plus la détection
des répétitions en tandem plutôt que leur contexte d’insertion lui aussi composé de motifs
répétés en tandem. Cependant, vraisemblablement que la combinaison de ces deux difficultés
rend leur détection encore plus difficile avec des lectures courtes.
Un résultat surprenant concerne l’homologie jonctionnelle. Le rappel de tous les outils
chute fortement en présence d’homologie de grande taille (>50 pb), mais même une petite
taille de 10 à 20 pb peut impacter le rappel de certains outils. Cela s’explique par le fait que
ces petites répétitions altèrent les signatures de mapping, comme elles altèrent les motifs
de k-mers dans MindTheGap (voir Chapitre 2 page 23). Or, notre étude montre que ces
répétitions ne sont pas si anecdotiques qu’on pouvait le penser, puisque près de 40 % des
insertions de ces catalogues possèdent de telles répétitions. Les algorithmes de détection
de variants de structure bénéficieraient donc de tenir compte de cette caractéristique des
variants.
Le résultat le plus marquant de ce benchmark est l’absence de résolution de la séquence
insérée pour la plupart des outils, même quand le site d’insertion est correctement prédit (voir
Table 3.1). Ainsi, à l’exception de MindTheGap, les outils ne fournissent la séquence insérée
que si elle est entièrement nouvelle et que sa taille est limitée à la taille des lectures. Dans
les autres cas, le rappel obtenu est très souvent de 0 %. MindTheGap, quant à lui, n’est pas
limité par la taille de la séquence et peut assembler certaines insertions duplicatives. Enfin,
aucun des quatre outils testés n’est capable d’assembler les insertions de type répétitions en
tandem, qui sont pourtant le type majoritaire dans les jeux de données réelles.
39
Rappel avec séquence résolue
GRIDSS Manta SvABA MindTheGap
Simulation de référence 81 100 96 100
Scenario 1 50 pb 56 100 100 100
Taille 500 pb 0 0 0 99
d’insertion 1 000 pb 0 0 0 98
Dup. dispersée 0 0 16 96
Scenario 2 Dup. tandem 0 0 0 0
Type Element transp. 0 0 61 58
d’insertion Rep. tandem (motif de 6 pb) 0 0 1 0
Rep. tandem (motif de 25 pb) 0 0 0 0
10 pb 99 100 92 0
Scenario 3 20 pb 100 100 78 0
Homologie 50 pb 6 46 10 0
jonct. 100 pb 0 11 0 0
150 pb 0 0 0 0
Non répété 77 97 93 83
Scenario 4 Rep. tandem (<300 pb) 77 98 97 73
Localisation Rep. tandem (>300 pb) 77 93 90 58
génomique SINE 77 99 94 53
LINE 76 97 95 89
Table 3.1 – Rappel avec séquence résolue de plusieurs logiciels de détection d’insertions
avec des lectures courtes en fonction des différents scénarios de simulation. Les cellules du
tableau sont colorées en fonction de la variation de la valeur de rappel de l’outil donné par
rapport au rappel obtenu avec la simulation de référence (première ligne, colorée en bleu) :
les cellules en rouge montrent une perte de rappel >10 %, les cellules en gris montrent peu
de différence ou une amélioration du rappel. Tableau repris de (Delage et al., 2020)
bien évidemment à la valeur de k et entraine une perte de spécificité du motif, ce qui peut
générer un grand nombre de faux positifs. Enfin, une duplication en tandem d’une séquence
à l’identique est un cas extrême d’une très grande homologie jonctionnelle, ce qui explique
l’absence du motif de k-mers pour ce type d’insertion. Ainsi, ces résultats suggèrent que
l’approche par k-mers du module Find de MindTheGap atteint ses limites pour détecter
certains types de sites d’insertions, et ceux-ci sont malheureusement les plus fréquents dans
les catalogues d’insertions humaines.
3.4 Conclusion
Ce travail représente la première caractérisation des variants de type insertion à grande
échelle chez l’homme. Il a été permis grâce à l’exhaustivité, la précision et la résolution de
séquence de catalogues de variants de structure obtenus seulement récemment grâce aux
nouvelles technologies de séquençage en lectures longues. Le résultat marquant de cette
caractérisation est la très grande diversité des insertions, tant du point de vue de leurs
caractéristiques génomiques, que du point de vue de leur capacité à être détectées avec les
lectures courtes.
40
Si l’inadéquation des lectures courtes pour la découverte de variants de structure en
général était déjà bien reconnue, notre étude met en lumière l’ampleur du problème pour les
variants de type insertion en particulier et plaide pour une meilleure évaluation des outils.
Les faibles valeurs de rappel obtenues dans notre benchmark contrastent avec les bonnes
performances montrées dans les publications des outils ou dans les précédents benchmarks
publiés pour évaluer les outils de détection des variants de structure. Cela provient d’un
manque de réalisme des données simulées et/ou d’un biais de représentativité des vrais
variants utilisés pour la validation. Ainsi, en ne considérant que certains types d’insertions,
et en particulier les plus faciles, les évaluations précédentes ont largement sur-estimé les
capacités des outils actuels.
Le tableau n’est pourtant pas si noir, car les performances des outils ne sont pas si
catastrophiques pour tous les types d’insertion et les insertions les plus difficiles ne sont pas
forcément celles qui sont le plus recherchées dans les applications biologiques. Cependant,
une évaluation plus réaliste, qui tient compte des différents types d’insertion notamment,
permettra à l’utilisateur de faire un choix d’outil plus éclairé en fonction de ses besoins.
Cette étude a permis également d’identifier plus précisément les facteurs responsables
d’une telle perte de rappel et de suggérer des pistes pour améliorer les algorithmes de détec-
tion des insertions avec des lectures courtes, et en particulier pour notre outil MindTheGap.
De plus, la variabilité des rappels entre les outils montre que bien souvent le signal du variant
est présent même avec des lectures courtes et qu’il peut être détecté, mais que les différents
outils n’appliquent pas les mêmes filtres. Cela donne donc de l’espoir de pouvoir faire mieux,
soit en améliorant les algorithmes grâce à l’analyse fine de ces résultats, soit en cherchant
une combinaison des différents outils qui prendrait le meilleur de chacun et maximiserait le
rappel. À l’heure actuelle, les méthodes dites meta ou de consensus qui combinent plusieurs
outils ont pour principal objectif de réduire le nombre de faux positifs et donc effectuent
simplement des intersections entre les ensembles de prédictions. Effectuer des unions per-
mettrait d’améliorer la sensibilité mais augmenterait considérablement le nombre de faux
positifs. Une étude, similaire à celle effectuée ici, non pas sur le rappel mais sur la précision
des outils permettrait probablement d’identifier des combinaisons d’outils plus intelligentes
et de réduire le nombre de faux positifs. Cependant, pour cela, il faudrait caractériser les
prédictions et non plus les vraies insertions et cela n’est possible que si les séquences insérées
ont été assemblées. Il faut donc travailler dans un premier temps à améliorer la résolution de
séquence des prédictions. Pour cela, notre benchmark montre que l’outil MindTheGap, et
en particulier son module d’assemblage Fill, est bien placé dans la compétition et pourrait
avantageusement s’intégrer dans une méthode qui combine plusieurs outils.
Toutefois, compte tenu de ces résultats, on peut se demander s’il est rentable de consacrer
des efforts sur ces méthodes alors que les lectures longues se généralisent. C’est une question
que nous discuterons dans le Chapitre 5 Section 5.1.
Ce travail représente une partie du travail de thèse de Wesley Delage que j’ai co-encadré.
Les résultats ont été publiés récemment dans BMC Genomics (Delage et al., 2020) et les
outils développés pour l’annotation des insertions et l’évaluation des outils sont diffusés à la
communauté sur github 4 .
41
Chapitre 4
Dans ce chapitre, nous présentons une méthode de génotypage des variants de structure
avec des lectures longues, qui a été développée dans le cadre de la thèse de Lolita Lecompte
(2017-2020) et qui a été publiée en 2020 dans la revue Bioinformatics (Lecompte et al.,
2020) (publication en Annexe A.3).
42
Les deux problèmes sont étroitement liés et sont parfois confondus. La confusion peut
provenir du fait que les méthodes de découvertes possèdent souvent un module de géno-
typage permettant d’enrichir l’information des variants découverts avec le génotype ou la
quantification des allèles dans les échantillons ayant permis leur découverte. Cependant, dans
cette situation, seuls les variants découverts sont génotypés. Dans l’hypothèse où les outils
de découvertes seraient extrêmement sensibles, on pourrait supposer que tous les variants
non détectés sont absents et déduire ainsi leur génotype. Cependant, nous avons vu que
cette hypothèse est loin d’être valide, il est donc important d’avoir une méthode capable
de vérifier la présence mais aussi l’absence d’un variant donné. De même, les méthodes
de découvertes produisant de nombreux faux positifs, elles demandent une profondeur de
séquençage plus importante pour valider la découverte d’un variant qu’elles n’en auraient
besoin pour simplement quantifier la présence d’un variant déjà validé par ailleurs. En effet,
dans de nombreuses applications, on effectue la découverte des variants avec un petit en-
semble d’individus fortement couverts, puis on les génotype dans des plus grands ensembles
d’individus plus faiblement couverts. Enfin, la comparaison de variants de structure prédits
indépendamment chez des individus différents est une tâche complexe, car il n’est pas rare
qu’un même variant soit prédit et décrit différemment par un même outil pour des jeux de
lectures différents. Le génotypage permet de s’affranchir de cette étape de comparaison et
d’identification des variants communs entre échantillons et d’obtenir des informations quan-
titatives pour chaque échantillon basées sur la même représentation des variants pour tous
les individus.
43
quement utilisés pour la découverte des variants (paires de reads discordantes, split-reads,
profondeur de séquençage). En conséquence, la séquence de l’allèle alternatif des variants
n’est pas représentée, ni utilisée pour la quantification. Cela induit un biais en faveur de
l’allèle de référence, et empêche de génotyper les variants de type insertions nouvelles pour
lesquels la séquence du variant alternatif est complètement absente du génome de référence.
La plupart des méthodes plus récentes ont pour objectif d’éliminer ce biais et ont bénéficié
de l’essor des représentations de génomes par des graphes. Dans un graphe de génomes ou de
variations, les sommets sont des séquences et les arêtes des chevauchements ou adjacences
de séquences observés dans au moins un allèle ou génome. Chaque allèle ou haplotype est
représenté par un chemin dans le graphe. Ainsi, contrairement à un génome de référence qui
ne représente qu’un seul allèle, avec ce type de graphe on dispose d’une structure pouvant
représenter tous les allèles connus d’un ensemble de variants. Lors du mapping des lectures
sur le graphe, il n’y a pas a priori de biais de mapping vers la référence. Les génotypeurs
les plus récents, tels que Paragraph (Chen et al., 2019) et graphTyper2 (Eggertsson et al.,
2019), exploitent ces structures de données. Cependant, la tâche de mapping des lectures
est plus complexe et plus coûteuse sur un graphe. Pour ne pas avoir à mapper la totalité
des lectures sur le graphe du génome complet, une pré-sélection des lectures d’intérêt est
effectuée à partir du mapping sur le génome de référence avant d’effectuer le mapping sur le
graphe. Le biais vers la référence n’est donc pas complètement éliminé.
L’ensemble des méthodes présentées précédemment ont été développées spécifiquement
pour les lectures courtes. En 2018, lorsque nous avons commencé ce travail, il existait plu-
sieurs méthodes de découvertes des variants de structure dédiées aux lectures longues, mais
aucune méthode de génotypage n’avait été décrite dans la littérature. Deux outils permet-
taient cependant le génotypage de variants de structure avec des lectures longues, mais ils
n’ont pas été dédiés à cette tâche : Sniffles (Sedlazeck et al., 2018b) est un outil de décou-
verte de variants de structure qui possède une option de génotypage (-Ivcf) non décrite dans
la littérature, et svviz2 (Spies et al., 2015) est un outil de visualisation d’alignements qui
estime des génotypes en sous-produit de sortie. Pour ces deux outils, là encore, il existe un
biais vers la référence puisque les lectures ne sont alignées que sur le génome de référence
pour le premier, et sont pré-sélectionnées avec la référence dans le deuxième.
Dans ce travail, nous avons donc proposé la première méthode dédiée au génotypage
des variants de structure avec des lectures longues. Cette méthode est implémentée dans le
logiciel SVJedi et a été publiée dans le journal Bioinformatics en 2020 (article en Annexe
A.3).
44
Figure 4.1 – Les différentes étapes de la méthode SVJedi. Figure extraite de (Lecompte
et al., 2020).
L’originalité de cette méthode est double : d’une part la représentation des allèles qui
donne le même poids à chacun des deux allèles de chaque variant de structure et d’autre part
l’alignement des lectures sur une petite sous partie du génome au lieu du génome complet. La
représentation des allèles permet de s’affranchir complètement du biais vers la référence et la
limitation de la séquence de référence aux séquences contenant des différences de structure
à génotyper permet de gagner en temps de calcul. Ces choix contrastent notamment avec
les méthodes existantes pour les lectures courtes, ils sont en effet largement influencés par
le type de lecture de séquençage à aligner : les lectures étant longues de plusieurs Kb,
elles s’alignent de manière plus spécifique que les lectures courtes de type Illumina. Ainsi,
elles sont moins sujettes aux problèmes de multi-mapping le long du génome de référence
(plusieurs positions génomiques équivalentes en terme de score d’alignement) et de sur-
mapping. Le sur-mapping est un problème qui apparait si la séquence de référence sur
laquelle les lectures sont mappées est incomplète : le mapper cherche à maximiser le nombre
de lectures alignées et force l’alignement de lectures provenant de ces régions manquantes à
d’autres loci représentés dans la séquence de référence. Ces alignements non légitimes issus
de multi-mapping ou de sur-mapping sont la principale source d’erreurs de génotypage et
doivent donc être évités.
En ne mappant que sur les séquences représentatives des allèles des variants de structure,
on gagne en temps de calcul car les lectures longues de part leur taux et leur type d’erreurs
de séquençage peuvent être coûteuses à aligner sur le génome entier. En contre partie, on
s’expose au problème de sur-mapping, c’est-à-dire d’aligner des lectures sur des séquences
représentatives de variant de structure alors qu’elles proviennent d’une autre région du
génome qu’on ne souhaite pas génotyper. Une étape importante de la méthode consiste
donc à identifier et éliminer ces alignements non légitimes parmi les résultats d’alignement.
Pour cela, nous avons implémenté des filtres basés sur le score de qualité de mapping et sur
l’analyse des coordonnées des alignements sur les lectures et les séquences de référence. En
particulier, un alignement, pour être considéré comme légitime et informatif, doit couvrir au
moins un point de cassure et doit être semi-global, c’est-à-dire qu’il implique les extrémités
soit de la lecture soit de la séquence de référence (voir les détails dans l’article présenté en
annexe).
4.2.2 Implémentation
La principale difficulté algorithmique dans cette méthode est l’alignement des lectures
longues sur les séquences représentatives des variants de structure. Les alignements recher-
chés sont des alignements locaux ou semi-globaux permettant de détecter des similarités
faibles jusqu’à 70 % d’identité environ avec un nombre important d’insertions et délétions
45
qui correspondent aux profil des erreurs de séquençage des données de lectures longues, telles
que celles obtenues avec les technologies PacificBioscience et Nanopore. Ce problème a fait
l’objet de recherche et développement dès l’apparition de ces lectures et est très bien résolu
par les outils de mapping actuels dédiés aux lectures longues. Parmi eux, minimap2 (Li,
2018) est l’un des plus populaires et probablement le plus rapide. Il est également facile
d’utilisation et d’intégration dans un pipeline. Nous avons donc estimé qu’il n’était pas
nécessaire de ré-implémenter un mapper dédié.
SVJedi a ainsi été implémenté en Python comme un pipeline faisant intervenir séquen-
tiellement 3 outils ou modules :
1. l’interprétation du fichier de variants de structure en entrée et la génération des
séquences représentatives,
2. le mapping des lectures longues sur les séquences représentatives avec minimap2,
3. l’analyse des alignements pour sélectionner les alignements informatifs et l’estimation
des génotypes.
L’implémentation modulaire permet de ne pouvoir effectuer que certaines parties du
pipeline et de ré-utiliser des résultats intermédiaires pour d’autres runs. Par exemple, la
représentation des variants de structure peut ainsi être ré-utilisée pour différents individus
à génotyper. Elle permet également de pouvoir interchanger le mapper de lectures longues.
4.3 Résultats
4.3.1 Données pour l’évaluation de la méthode
Au moment du développement de la méthode, il n’existait pas de jeux de données réelles,
du moins pour un génome complexe tel que le génome humain, permettant d’évaluer préci-
sément la qualité de nos estimations des génotypes. En effet, nous avons besoin de données
de séquençage en lectures longues pour au moins un individu pour lequel on connait les
génotypes d’un ensemble représentatif et bien caractérisé de variants de structure. Comme
nous l’avons vu précédemment, des ensembles exhaustifs, ou au moins représentatifs, et bien
caractérisés de variants de structure sont rares ou très récents. Chez l’homme, la base de
données dbVar est biaisée vers les variants de structure de type délétions, identifiées prin-
cipalement avec des lectures courtes, et ce n’est qu’en 2019 qu’ont été rendus disponibles
les callsets du HGSV et de GiaB pour quatre individus. Parmi ces différents jeux de don-
nées, seul le jeu de données du GiaB dispose de l’information de génotype pour un individu
(HG002) et nous l’avons utilisé dans un second temps pour évaluer notre méthode.
Ainsi, dans un premier temps, nous avons donc généré des données simulées pour déve-
lopper, optimiser et valider la méthode. Le principe de la génération d’un jeu de données
simulées est le suivant : on définit un ensemble de variants de structure à génotyper sur
un génome donné, on génère deux haplotypes à partir de ce génome en incorporant dans
chacun un sous-ensemble de l’ensemble des variants de structure à génotyper et on simule
un séquençage de lectures longues sur ce génome diploïde synthétique. Pour un variant de
structure donné, son absence ou sa présence dans un ou les deux haplotypes définit son
génotype dans l’individu simulé. Nous avons effectué ces simulations sur le chromosome 1
humain et nous avons sélectionné des variants de structure déjà caractérisés dans ce chro-
mosome grâce à la base de données de dbVar lorsque c’était possible, c’est-à-dire pour les
variants de structure de type délétion. La majorité des simulations a été effectué avec ce
46
type de variant de structure qui est également le plus fréquent et le mieux représenté dans
les bases de données.
Cette approche de simulation permet non seulement d’évaluer précisément la qualité
des estimations des génotypes puisque la vérité est connue, mais également d’évaluer la
robustesse de la méthode vis-à-vis de différentes caractéristiques des données que l’ont peut
facilement faire varier dans le processus de simulation. Ainsi, nous avons fait varier dans nos
simulations le type de variant de structure, la technologie de séquençage, la profondeur de
séquençage, le taux d’erreurs de séquençage et la précision de caractérisation des variants.
L’évaluation se base sur deux métriques : le taux de génotypage et la précision de génotypage.
Le premier est le pourcentage de variants pour lesquels un génotype a été assigné (dans
SVJedi, la raison pour laque un variant n’est pas génotypés est un nombre insuffisant de
lectures mappées de manière informative sur ses allèles). La précision de génotypage est le
pourcentage de variants, parmi les variants génotypés, dont l’estimation du génotype est
correcte.
47
4.3.3 Comparaison avec d’autres approches
Nous avons comparé ces résultats avec ceux obtenus par d’autres approches, et dans un
premier temps les deux outils directement concurrents utilisant des lectures longues : l’outil
de découverte de variants de structure Sniffles (Sedlazeck et al., 2018b) avec l’option -Ivcf,
et l’outil de visualisation des variants de structure svviz2 (Spies et al., 2015). Sur le jeu de
données réelles du GiaB, ces deux outils estiment un génotype pour presque tous les variants
(taux de génotypage de 100 % ou presque), mais au prix d’une précision bien plus faible
que SVJedi : 82 et 66 % de précision pour Sniffles et svviz2 respectivement, comparé à 92
% pour SVJedi.
SVJedi se démarque aussi nettement de ses concurrents concernant le temps de calcul :
pour génotyper les 12 745 variants de structure de GiaB avec un re-séquençage PacBio à
30X d’un génome humain, le temps de calcul est seulement de 2h25 avec 40 cœurs. La
très grande majorité du temps est pris par l’alignement des lectures avec minimap2 (2h15).
Surtout, SVjedi est 7 fois plus rapide que Sniffles (17h15) et 50 fois plus rapide que svviz2
(plus de 5 jours, mais il ne dispose pas d’une version parallélisée).
Pour ce jeu de données, il existe également des données de lectures courtes Illumina. La
comparaison de SVJedi avec une approche de génotypage avec des lectures courtes montre
clairement l’apport des lectures longues pour la qualité du génotypage, puisque l’un des
meilleurs génotypeurs de lectures courtes, SVtyper (Chiang et al., 2015), a assigné un mau-
vais génotype à plus de la moitié des variants (précision de 46 %).
Enfin, la dernière comparaison a pour objectif d’évaluer l’apport d’une approche dédiée
de génotypage par rapport à une approche de découverte qui génotype les variants détectés.
Dans cette approche, seuls les variants correctement découverts peuvent être génotypés et
le taux de génotypage correspond donc au rappel de détection des méthodes. Pour ce jeu de
données, en moyenne seulement la moitié des variants de structure ont pu être correctement
prédits par Sniffles ou PBsv. Plus surprenant, parmi les variants détectés la précision de
génotypage obtenue est plutôt faible, 44 et 78 % respectivement pour Sniffles et PBsv.
Cette expérience montre donc que la tâche de génotypage est bien différente de celle de la
découverte et qu’elle mérite d’avoir ses méthodes dédiées.
L’ensemble de ces résultats sont détaillés par type de variant de structure dans la Table
4.1.
48
Délétions Insertions
Outil précision taux précision taux temps
SVJedi 91,7 85,8 92,5 93,6 2h25m
Sniffles–Ivcf 82,5 99,9 81,7 99,8 17h16m
svviz2 72,5 100 61,0 100 5 jours∗
SVtyper (Illumina) 46,5 99,2 - - 5h32m
Sniffles (découverte) 48,7 52,4 39,8 44,8 18h04m
pbsv 90,1 72,7 68,8 59,8 5h29m
Table 4.1 – Comparaison de plusieurs outils et approches pour le génotypage des 12 745
délétions et insertions du catalogue du GiaB chez l’individu HG002. Trois approches sont
comparées : des outils de génotypage pour lectures longues (trois premiers outils), un outil
de génotypage pour lectures courtes (SVTyper) et des outils de découverte pour lectures
longues (deux derniers outils). À l’exception de SVtyper qui utilise un jeu de données de
séquençage Illumina 30X, tous les autres outils ont été exécutés avec un jeu de données de
lectures longues PacBio 30X. Les temps d’exécution ont été mesurés sur un nœud de calcul
de 40 CPU. ∗ svviz2 n’est pas parallélisé. Table reprise de (Lecompte et al., 2020).
4.5 Conclusion
Alors que la découverte des variants de structure est nettement améliorée par les techno-
logies de séquençage en lectures longues, les catalogues de variants de structure s’enrichissent
pour de nombreux organismes. Le génotypage de ces variants dans des nouveaux séquençages
individuels est devenu une problématique importante et nous avons assisté au développe-
49
ment de nombreuses méthodes bioinformatiques ces dernières années pour effectuer cette
tâche. Cependant, aucune méthode utilisant des données de lectures longues n’existait et
nous avons ainsi proposé la première, appelée SVJedi.
Cette méthode est dédiée aux lectures longues et exploite leur spécificité : la grande
longueur des lectures rend leur mapping bien plus spécifique que les lectures courtes et
nous permet de s’affranchir du mapping sur le génome complet. À la place, les différents
allèles de chaque variant sont explicitement représentés, ce qui permet d’éviter le biais vers la
référence qui est un défaut classique des méthodes utilisant les lectures courtes. Les résultats
obtenus sur des données simulées et réelles humaines montrent que cette nouvelle approche de
génotypage est efficace et rapide. Ce travail constitue le travail de thèse de Lolita Lecompte,
il a été publié récemment dans la revue Bioinformatics (Lecompte et al., 2020) (voir Annexe
A.3). Le logiciel est diffusé sous licence libre sur github 1 et distribué dans Bioconda 2 . Les
améliorations avec la représentation en graphes de séquences sont en cours de validation et
de diffusion avec l’outil SVJedi. Une des perspectives de ce travail que j’envisage à court
terme exploitera encore la représentation par graphe de séquences pour estimer non plus les
génotypes, mais les haplotypes de variants de structure proches, c’est-à-dire comment les
différents allèles des variants proches sont liés sur les chromosomes homologues.
1. https://github.com/llecompte/SVJedi
2. https://anaconda.org/bioconda/svjedi
50
Chapitre 5
Discussion et perspectives
Dans ce document, j’ai choisi de présenter parmi mes différentes contributions en bio-
informatique des séquences, celles qui portent sur une thématique particulière : l’étude des
variants de structure et ses problématiques méthodologiques. C’est une thématique qui me
motive particulièrement depuis de nombreuses années, et sur laquelle je veux poursuivre mes
recherches futures.
C’est un domaine de recherche qui est en plein essor en ce moment car, d’une part
les technologies de séquençage sont plus adaptées pour détecter les variants de structure,
et d’autre part, après avoir étudié en profondeur les variations ponctuelles, les biologistes
se tournent désormais vers des variants plus complexes pour expliquer la variabilité phé-
notypique restante. Ainsi, de plus en plus d’études montrent l’impact de ces variants sur
des traits phénotypiques. Un exemple récent d’une telle étude a été publié récemment sur
la tomate (Alonge et al., 2020), cette étude identifie notamment de nombreux variants de
structure qui modifient des traits d’intérêt agronomique tels que la saveur du fruit, sa taille
ou la productivité de la plante. La recherche méthodologique pour l’étude de ces variants
est par conséquent un domaine très compétitif, et qui évolue très vite du fait de l’évolution
des technologies de séquençage. Ainsi, dans ce chapitre, je discuterai dans un premier temps
l’évolution de l’utilisation des différentes technologies de séquençage pour ces questions avant
de présenter les perspectives de mes travaux et mes axes de recherche actuels et futurs.
5.1 Est-ce la fin des lectures courtes pour étudier les variants
de structure ?
La détection des variants de structure avec des données de séquençage est un problème
complexe et non résolu pour tous les types de variants et tous les types de données de
séquençage. Depuis plusieurs années, nous avons observé et quantifié l’inadéquation des lec-
tures courtes pour détecter et analyser ce type de variations génomiques. Le chapitre 3 en
particulier, montre bien la faible sensibilité de ces données (et de leurs méthodes associées)
pour une majorité des variants de type insertion. L’amélioration des technologies de séquen-
çage permettant d’obtenir des lectures longues de plusieurs dizaines de Kilobases voire des
Mégabases a permis de considérablement augmenter la sensibilité de détection et réduire
le nombre de fausses prédictions. La réduction actuelle des taux d’erreur de ces lectures
promet d’améliorer encore les performances de détection des variants de structure avec ces
lectures. Cependant, ces technologies restent encore très coûteuses et pour de nombreuses
applications, les technologies moins chères seront encore probablement largement utilisées.
51
Par exemple, les études de recherche d’association avec des phénotypes et les études
de génomique des populations nécessitent d’analyser des grands nombres d’individus, typi-
quement plusieurs centaines. À part pour des projets sur quelques espèces modèles ou des
espèces à fort intérêt économique, le coût des technologies ONT et PacBio (et encore plus
pour les données PacBio HiFi) reste pour le moment prohibitif pour de telles études (Coster
et al., 2021). Ainsi, une stratégie moins coûteuse pour des projets à grands nombres d’échan-
tillons est de partitionner les échantillons en deux ensembles. Le premier est composé d’un
petit effectif d’individus, le plus représentatif possible de la diversité génomique à étudier,
ces individus seront séquencés avec des lectures longues à forte couverture pour découvrir
de nouveaux variants de structure et constituer un catalogue de qualité. Dans le second
ensemble, de taille plus grande, les individus sont séquencés avec une technologie moins
chère, les lectures courtes. Ces lectures permettent de quantifier dans des grandes popula-
tions le sous-ensemble de variants de structure découverts dans le premier ensemble. Ainsi,
de grandes quantités de lectures courtes vont continuer à être générées. Dans un tel plan
expérimental, a priori les lectures courtes ne sont pas utilisées pour découvrir de nouveaux
variants de structure, seulement pour le génotypage. Il est donc important d’améliorer au
moins les méthodes de génotypage avec ces données, qui comme nous l’avons montré dans
le Chapitre 4 sont encore peu efficaces. Il est également probable que pour certaines appli-
cations la découverte de variants de structure avec ces lectures s’avère utile. Par exemple,
une étape de découverte avec toutes les lectures courtes en amont du séquençage en lectures
longues peut être envisagée pour optimiser la sélection les individus du premier ensemble
(c’est le cas dans l’étude effectuée sur 900 échantillons de tomates de Alonge et al. (2020),
dont 100 échantillons ont été sélectionnés grâce à l’outil SVcollector (Sedlazeck et al., 2018a)
pour maximiser la diversité structurale).
Une autre alternative pour ce type de projet est d’utiliser des données de séquençage
"linked-reads", qui permettent d’associer une information longue distance à des données de
type lectures courtes. Les premières technologies à produire ce type de données restaient en-
core trop chères (par exemple 10X chromium genomics) pour des projets à grands nombres
d’échantillons. Mais récemment, d’autres technologies ont été développées, telles que TELL-
seq (Chen et al., 2020), stLFR (Wang et al., 2019) et Haplotagging (Meier et al., 2021). En
particulier, Haplotagging est un nouveau protocole libre qui permet de séquencer avec un
sur-coût très faible par rapport à un séquençage Illumina classique des centaines d’individus
avec une information longue distance. Dans la publication, ils démontrent qu’avec des cen-
taines d’individus séquencés à seulement 2 à 3 x chacun, ils peuvent reconstruire de grands
haplotypes et identifier des variants de structure pour des analyses poussées de génomique
des populations. C’est une alternative très intéressante dans le domaine de la génomique
environnementale où l’intérêt économique est plus faible et les quantités d’ADN sont plus
limitées.
Dans le domaine médical, le diagnostic de nombreuses maladies par séquençage plein
génome est en développement notamment en France grâce au Plan France Médecine Géno-
mique, mais la question du coût est cruciale pour adopter ces tests en routine. Ainsi, c’est
actuellement et probablement pour encore quelques années la technologie Illumina de lec-
tures courtes qui est préconisée. Ainsi, dans ce domaine, malgré les limitations intrinsèques
des courtes lectures, il y a un vrai besoin d’améliorer les méthodes de détection avec des
lectures courtes.
Même si les lectures longues ont un avantage certain par rapport aux lectures courtes
(et cet avantage ne va faire que croître avec l’amélioration de leur qualité et la diminution
du coût), je pense que les méthodes développées pour les lectures courtes seront encore
52
beaucoup utilisées et méritent d’être encore améliorées. Ainsi, je propose par la suite quelques
pistes d’amélioration et développements sur ces données et en particulier les données linked-
reads. Concernant les lectures longues, elles vont bien entendu prendre une part de plus en
plus importante dans les études des variants de structure et donc naturellement dans mes
problématiques de recherche méthodologique. Les méthodes actuelles sont encore jeunes et
ont encore des limitations sur lesquelles du travail reste à faire.
53
Genomics pour des problèmes de propriété intellectuelle, l’engouement pour cette technologie
a largement décru et ces outils ont arrêté d’être maintenus. Ce n’est que récemment que
d’autres technologies ont repris le flambeau et suscitent de nouveau l’intérêt notamment
pour des projets de génomique des populations de génomes non modèles. Or, les outils
actuels sont peu adaptés à des génomes non modèles, qui peuvent présenter une diversité
génétique et un taux d’hétérozygotie accrus, et ont probablement été sur-paramétrés pour
des données humaines de type 10X. D’autre part, ils sont très gourmands en ressources de
calcul et notamment en utilisation de la mémoire vive, ce qui les rend parfois inutilisables.
Ainsi un premier axe de recherche porte sur ces problèmes de performances en temps et
mémoire. Nous travaillons sur la représentation de ces données en mémoire et notamment des
informations spécifiques à ces données : les barcodes. En effet, de nombreuses applications
demandent d’extraire ou de comparer des ensembles de barcodes en fonction des régions
génomiques où sont mappées leurs lectures, ou bien de récupérer l’ensemble de lectures
contenant un barcode donné. Nous développons LRez, une librairie C++ et une suite d’outils,
qui permet d’indexer les lectures par barcodes et de faire des opérations simples mais rapides
sur les lectures en fonction de leurs barcodes. LRez est disponible sur github et bioconda
et est en cours de publication (Morisse et al., 2021b). C’est une brique de base essentielle
pour ensuite développer d’autres méthodes utilisant efficacement ses données.
Parmi ces méthodes, je travaille actuellement sur deux applications : l’assemblage local
pour boucher les trous d’un assemblage non fini ou pour reconstruire la séquence de points
de cassure, et la découverte de variants de structure. Dans le premier cas, l’information
des barcodes est utilisée pour sélectionner, pour chaque séquence à assembler, un sous-
ensemble de lectures susceptibles de provenir du locus génomique en question, en fonction
des ensembles de barcodes observés aux extrémités de la séquence. La réduction du jeu de
lectures permet de réduire la complexité du graphe d’assemblage. Le choix de la structure de
données pour le graphe d’assemblage se pose alors : un graphe de de Bruijn en ré-utilisant
directement MindTheGap, ou bien un graphe de chevauchements de lectures qui est plus
classiquement utilisé pour des lectures longues puisque la taille des données le permet ici. La
première solution pose des problèmes de paramétrage du graphe et notamment de la taille des
k-mers, car les données de séquençage en entrée ont une profondeur qui dépend de l’efficacité
de la sélection des barcodes. Alors que la deuxième solution permet des chevauchements de
lectures de taille variable, ce qui s’adapte mieux aux variabilités de couverture. Nous étudions
et comparons ces deux approches dans l’outil MTG-link 1 qui est actuellement en cours de
développement.
Dans le problème de la découverte des variants de structure avec ces données, l’informa-
tion longue distance portée par les barcodes fournit un signal spécifique et une valeur ajoutée
par rapport aux lectures courtes seules notamment pour les grands variants de structure,
et en particulier les grandes inversions. Le signal de barcodes classiquement recherché est
une paire de régions distantes sur le génome qui partagent un nombre de barcodes com-
muns plus grand qu’attendu. L’identification de ce signal soulève alors deux problèmes, l’un
algorithmique sur le comptage efficace des barcodes partagés, et l’autre statistique sur la
modélisation et l’évaluation statistique de ces comptages. Le premier se résout grâce aux
structures d’indexation développées dans la librairie LRez (actuellement implémenté dans
le prototype LEVIATHAN (Morisse et al., 2021a)), et le second est encore à explorer. Des
développements sont également encore nécessaires pour l’étape suivante, qui est la caracté-
risation des variants de structure et de leurs points de cassure dans les régions identifiées
1. https://github.com/anne-gcd/MTG-Link
54
puisque le signal de barcodes seul donne une information fiable mais peu précise, et égale-
ment pour la découverte des variants de structure qui ne possèdent pas exactement deux
points de cassure, comme les insertions (un seul point de cassure) et les transpositions (trois
points de cassure).
Enfin, les données issues de la technologie Haplotagging comportent une difficulté sup-
plémentaire. Chaque individu étant séquencé à faible couverture (2-3 X), la découverte de
variants de structure se fait sur l’ensemble poolé des individus. Les méthodes devront certai-
nement être adaptées pour gérer ce polymorphisme dans les données qui peut gêner notam-
ment les assemblages locaux. L’étape suivante est alors de revenir à l’information individuelle
pour génotyper les variants de structure découverts pour chaque individu. Cette probléma-
tique rejoint celle adressée dans le chapitre 4 et pourra bénéficier des travaux effectués sur
la représentation des variants de structure et de leurs allèles. Dans ce cas, l’utilisation des
informations de barcode sera capitale car la couverture physique en molécules est bien plus
importante que la couverture en lectures courtes.
5.2.3 Affiner les points de cassure avec des données de lectures longues
Concernant l’utilisation des lectures longues, certains types de variants restent encore
mal prédits par les outils actuels ou bien leur description n’est pas suffisamment précise.
Notre benchmark sur divers types d’insertions, présenté dans le Chapitre 3, montre par
exemple que l’outil Sniffles renvoie des séquences insérées de mauvaise qualité avec plus de
20 % de divergence avec la séquence simulée et n’est pas capable de détecter les homologies
jonctionnelles de grande taille ou les duplications en tandem (voir le matériel supplémen-
taire 2 de (Delage et al., 2020)). De même, les inversions bordées par des répétitions inversées
(ce qui est très fréquent car cela fait partie du mécanisme de génération de ces inversions)
sont encore très mal prédites avec des lectures longues (Mahmoud et al., 2019). Dans ces
deux cas, le simple mapping des lectures longues sur le génome de référence ne suffit pas pour
correctement prédire ces variants. Plusieurs axes d’amélioration sont possibles comme faire
de l’assemblage local des lectures ou construire des séquences consensus pour les séquences
insérées et les points de cassure des variants, ou encore rechercher de manière explicite la pré-
sence d’homologie jonctionnelle puisque nous avons vu que c’est une caractéristique présente
dans plus de 40 % des insertions.
Le problème de détection des inversions m’intéresse particulièrement, car ce sont des
variants qui jouent un rôle important dans l’adaptation, l’évolution et la spéciation des es-
pèces, notamment dans des espèces d’insectes sur lesquelles j’entretiens des collaborations
depuis plusieurs années. Dans le cas des inversions avec des répétitions inversées, la détection
précise des divergences de séquence si elles existent entre les deux copies répétées et avec
les lectures pourrait permettre de mieux assigner les lectures et d’identifier précisément le
point de recombinaison dans la séquence répétée. Cela fait notamment écho à mes travaux
de thèse de génomique comparée sur l’affinement des points de cassure de réarrangements
évolutifs. J’avais développé une méthode appelée Cassis, basée sur la segmentation du si-
gnal des alignements des séquences flanquant les points de cassure (Lemaitre et al., 2008;
Baudet et al., 2010). Par la suite, j’ai ré-utilisé et adapté cette méthode pour des données
de séquençage dans le cadre d’un projet d’analyse de la biologie de virus intégrés dans le
génome d’une guêpe parasitoïde. C’est un système où des répétitions directes du génome
viral permettent à celui-ci de s’exciser du génome hôte pour former des cercles viraux in-
dépendants qui sont ensuite exportés pour défendre les œufs de la guêpe contre le système
2. https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-020-07125-5#Sec34
55
immunitaire de l’organisme qu’elle parasite. L’utilisation de cette approche de segmentation
sur le signal des mutations ponctuelles a permis l’identification précise des points d’excision
et l’analyse de leur répartition pour mieux comprendre ce mécanisme moléculaire (Legeai
et al., 2020). Dans le cas de lectures longues, la problématique algorithmique sera d’obtenir
des alignements de séquences précis malgré les forts taux d’erreurs de séquençage de ces
données et de distinguer ces dernières du polymorphisme des génomes re-séquencés.
Ces pistes d’amélioration visent à préciser la description des variants et affiner la posi-
tion du ou des points de cassure. C’est une problématique cruciale lorsqu’on veut ensuite
comparer des prédictions obtenues avec des outils différents ou chez des individus différents
ou encore pour le génotypage des variants dans des populations.
56
des défis méthodologiques. Les délétions et insertions sont les variants de structure les plus
faciles à représenter, ils forment des bulles avec deux chemins alternatifs dans le graphe. Les
variants équilibrés, comme les inversions, les translocations ou les transpositions sont plus
difficiles à implémenter, notamment car ils peuvent créer des cycles. Par exemple, la repré-
sentation GRAF de l’entreprise SevenBridges est par nature un graphe acyclique (DAG) et
les segments inversés sont dupliqués dans le graphe pour éviter les cycles (Rakocevic et al.,
2019). L’outil de génotypage de VG-toolkit montre également des performances nettement
moins bonnes pour les inversions que pour les insertions ou délétions (Hickey et al., 2020). À
cela s’ajoutent également des problématiques d’explosion combinatoire du nombre de che-
mins et de passage à l’échelle lorsque le nombre d’individus ou de variants augmentent dans
le graphe. Ainsi, cette thématique offre de nombreux axes de recherche pour l’étude des
variants de structure, avec des défis méthodologiques importants, que je voudrais explorer
sur le plus long terme.
57
d’hybridations entre espèces parentales différentes, rendant le choix d’un génome de référence
difficile et risqué. L’approche que nous privilégierons sera celle de la représentation des
génomes par un graphe de séquences. Des questions méthodologiques sur la découverte et le
génotypage des variants de structure avec des lectures longues et des données linked-reads
font naturellement partie de ce projet.
De manière plus générale, les données et les méthodes commencent tout juste à être
efficaces pour détecter les variants de structures dans les génomes et de nombreuses ques-
tions biologiques peuvent enfin être abordées. Parmi elles, les mécanismes d’apparition et
de maintien de tels variants, leurs impacts fonctionnels sur la biologie des organismes, ainsi
que leurs impacts sur l’évolution des génomes et des populations, sont des problématiques
biologiques qui m’intéressent. J’aimerais les aborder par des recherches méthodologiques (et
les problèmes ne manquent pas comme nous l’avons vu précédemment) mais également par
l’analyse des données qui sont et seront produites en toujours plus grande quantité.
58
Bibliographie
C. Baudet, C. Lemaitre, Z. Dias, C. Gautier, E. Tannier, and M.-F. Sagot. Cassis : Detection
of genomic rearrangement breakpoints. Bioinformatics, 26(15) :1897–1898, Jun 2010. doi :
10.1093/bioinformatics/btq301. [lien].
G. Benoit, D. Lavenier, C. Lemaitre, and G. Rizk. Bloocoo, a memory efficient read corrector.
European Conference on Computational Biology (ECCB), Sept. 2014. [lien]. Poster.
59
positional de bruijn graph assembly. Genome Research, 2017. doi : 10.1101/gr.222109.117.
[lien].
Z. Chen, L. Pham, T.-C. Wu, G. Mo, Y. Xia, P. L. Chang, D. Porter, T. Phan, H. Che,
H. Tran, V. Bansal, J. Shaffer, P. Belda-Ferre, G. Humphrey, R. Knight, P. Pevzner,
S. Pham, Y. Wang, and M. Lei. Ultralow-input single-tube linked-read library method
enables short-read second-generation sequencing systems to routinely generate highly ac-
curate and economical long-range sequencing information. Genome Research, 30(6) :
898–909, jun 2020. doi : 10.1101/gr.260380.119.
R. Chikhi and G. Rizk. Space-efficient and exact de bruijn graph representation based on a
bloom filter. Algorithms Mol Biol, 8(1) :22, 2013. doi : 10.1186/1748-7188-8-22. [lien].
60
S. Daval, A. Belcour, K. Gazengel, L. Legrand, J. Gouzy, L. Cottret, L. Lebreton, Y. Aigu,
C. Mougel, and M. J. Manzanares-Dauleux. Computational analysis of the plasmodio-
phora brassicae genome : mitochondrial sequence description and metabolic pathway da-
tabase design. Genomics, 111(6) :1629–1640, dec 2019. doi : 10.1016/j.ygeno.2018.11.013.
W. J. Delage, J. Thevenon, and C. Lemaitre. Towards a better understanding of the low
recall of insertion variants with short-read based variant callers. BMC Genomics, 21(1),
nov 2020. doi : 10.1186/s12864-020-07125-5.
E. Drezen, G. Rizk, R. Chikhi, C. Deltel, C. Lemaitre, P. Peterlongo, and D. Lavenier. Gatb :
Genome assembly & analysis tool box. Bioinformatics, 30(20) :2959–2961, Oct 2014. doi :
10.1093/bioinformatics/btu406. [lien].
H. P. Eggertsson, S. Kristmundsdottir, D. Beyter, H. Jonsson, A. Skuladottir, M. T. Hardar-
son, D. F. Gudbjartsson, K. Stefansson, B. V. Halldorsson, and P. Melsted. GraphTyper2
enables population-scale genotyping of structural variation using pangenome graphs. Na-
ture Communications, 10(1), nov 2019. doi : 10.1038/s41467-019-13341-9.
A. C. English, W. J. Salerno, O. A. Hampton, C. Gonzaga-Jauregui, S. Ambreth, D. I.
Ritter, C. R. Beck, C. F. Davis, M. Dahdouli, et al. Assessing structural variation in a
personal genome-towards a human reference diploid genome. BMC Genomics, 16(1) :286,
Apr 2015. doi : 10.1186/s12864-015-1479-3. [lien].
D. Feldman, D. J. Kowbel, A. Cohen, N. L. Glass, Y. Hadar, and O. Yarden. Identifi-
cation and manipulation of neurospora crassa genes involved in sensitivity to furfural.
Biotechnology for Biofuels, 12(1), sep 2019. doi : 10.1186/s13068-019-1550-4.
R. R. Fuentes, D. Chebotarov, J. Duitama, S. Smith, J. F. D. la Hoz, M. Mohiyuddin,
R. A. Wing, K. L. McNally, T. Tatarinova, A. Grigoriev, R. Mauleon, and N. Alexandrov.
Structural variants in 3000 rice genomes. Genome Research, 29(5) :870–880, apr 2019.
doi : 10.1101/gr.241240.118.
I. Gabur, H. S. Chawla, R. J. Snowdon, and I. A. P. Parkin. Connecting genome structural
variation with complex traits in crop plants. Theoretical and Applied Genetics, 132(3) :
733–750, Mar. 2019. doi : 10.1007/s00122-018-3233-0.
E. Garrison, J. Sirén, A. M. Novak, G. Hickey, J. M. Eizenga, E. T. Dawson, W. Jones,
S. Garg, C. Markello, M. F. Lin, B. Paten, and R. Durbin. Variation graph toolkit improves
read mapping by representing genetic variation in the reference. Nature biotechnology, 36 :
875–879, Oct. 2018. ISSN 1546-1696. doi : 10.1038/nbt.4227.
A. Gouin, F. Legeai, P. Nouhaud, A. Whibley, J.-C. Simon, and C. Lemaitre. Whole-genome
re-sequencing of non-model organisms : lessons from unmapped reads. Heredity, 114(5) :
494–501, May 2015. doi : 10.1038/hdy.2014.85. [lien].
C. Guyomar, F. Legeai, E. Jousselin, C. Mougel, C. Lemaitre, and J.-C. Simon. Multi-
scale characterization of symbiont diversity in the pea aphid complex through me-
tagenomic approaches. Microbiome, 6(1) :181, Oct 2018. ISSN 2049-2618. doi :
10.1186/s40168-018-0562-9. [lien].
C. Guyomar, W. Delage, F. Legeai, C. Mougel, J.-C. Simon, and C. Lemaitre. MinYS : mine
your symbiont by targeted genome assembly in symbiotic communities. NAR Genomics
and Bioinformatics, 2(3), jul 2020. doi : 10.1093/nargab/lqaa047.
61
I. Hajirasouliha, F. Hormozdiari, C. Alkan, J. M. Kidd, I. Birol, E. E. Eichler, and S. C.
Sahinalp. Detection and characterization of novel sequence insertions using paired-end
next-generation sequencing. Bioinformatics, 26(10) :1277–1283, May 2010. doi : 10.1093/
bioinformatics/btq152. [lien].
D. Heller and M. Vingron. Svim : Structural variant identification using mapped long
reads. Bioinformatics (Oxford, England), Jan. 2019. ISSN 1367-4811. doi : 10.1093/
bioinformatics/btz041.
M. Kirkpatrick. How and why chromosome inversions evolve. PLoS Biology, 8(9) :e1000501,
sep 2010. doi : 10.1371/journal.pbio.1000501.
62
H. Lee and M. C. Schatz. Genomic dark matter : the reliability of short read mapping
illustrated by the genome mappability score. Bioinformatics, 28(16) :2097–2105, jul 2012.
doi : 10.1093/bioinformatics/bts330.
F. Legeai, B. F. Santos, S. Robin, A. Bretaudeau, R. B. Dikow, C. Lemaitre, V. Jouan,
M. Ravallec, J.-M. Drezen, D. Tagu, F. Baudat, G. Gyapay, X. Zhou, S. Liu, B. A. Webb,
S. G. Brady, and A.-N. Volkoff. Genomic architecture of endogenous ichnoviruses reveals
distinct evolutionary pathways leading to virus domestication in parasitic wasps. BMC
Biology, 18(1), jul 2020. doi : 10.1186/s12915-020-00822-3. [lien].
C. Lemaitre, E. Tannier, C. Gautier, and M.-F. Sagot. Precise detection of rearrangement
breakpoints in mammalian chromosomes. BMC Bioinformatics, 9(1) :286, Jun 2008. doi :
10.1186/1471-2105-9-286. [lien].
H. Li. Minimap2 : pairwise alignment for nucleotide sequences. Bioinformatics, 34(18) :
3094–3100, may 2018. doi : 10.1093/bioinformatics/bty191.
H. Li, X. Feng, and C. Chu. The design and construction of reference pangenome graphs
with minigraph. Genome Biology, 21(1), oct 2020. doi : 10.1186/s13059-020-02168-z.
S. Li, R. Li, H. Li, J. Lu, Y. Li, L. Bolund, M. H. Schierup, and J. Wang. Soapindel : efficient
identification of indels from short paired reads. Genome Res, 23(1) :195–200, Jan 2013.
doi : 10.1101/gr.132480.111. [lien].
G. A. Logsdon, M. R. Vollger, and E. E. Eichler. Long-read human genome sequencing
and its applications. Nature Reviews Genetics, 21(10) :597–614, jun 2020. doi : 10.1038/
s41576-020-0236-x.
M. Mahmoud, N. Gobet, D. I. Cruz-Dávalos, N. Mounier, C. Dessimoz, and F. J. Sedlazeck.
Structural variant calling : the long and the short of it. Genome Biology, 20(1), nov 2019.
doi : 10.1186/s13059-019-1828-7.
J. I. Meier, P. A. Salazar, M. Kučka, R. W. Davies, A. Dréau, I. Aldás, O. B. Power, N. J.
Nadeau, J. R. Bridle, C. Rolian, N. H. Barton, W. O. McMillan, C. D. Jiggins, and
Y. F. Chan. Haplotype tagging reveals parallel formation of hybrid races in two butterfly
species. Proceedings of the National Academy of Sciences, 118(25) :e2015005118, jun 2021.
doi : 10.1073/pnas.2015005118.
M. Mohiyuddin, J. C. Mu, J. Li, N. B. Asadi, M. B. Gerstein, A. Abyzov, W. H.
Wong, and H. Y. Lam. MetaSV : an accurate and integrative structural-variant cal-
ler for next generation sequencing. Bioinformatics, 31(16) :2741–2744, apr 2015. doi :
10.1093/bioinformatics/btv204.
P. Morisse, F. Legeai, and C. Lemaitre. LEVIATHAN : efficient discovery of large structural
variants by leveraging long-range information from linked-reads data. bioRxiv, mar 2021a.
doi : 10.1101/2021.03.25.437002.
P. Morisse, C. Lemaitre, and F. Legeai. Lrez : C++ api and toolkit for analyzing and
managing linked-reads data. arXiv, Mar. 2021b. [lien].
S. Nurk, S. Koren, A. Rhie, M. Rautiainen, A. V. Bzikadze, A. Mikheenko, M. R. Vollger,
N. Altemose, L. Uralsky, A. Gershman, et al. The complete sequence of a human genome.
bioRxiv, may 2021. doi : 10.1101/2021.05.26.445798.
63
D. Ottaviani, M. LeCain, and D. Sheer. The role of microhomology in genomic structural
variation. Trends in Genetics, 30(3) :85–94, mar 2014. doi : 10.1016/j.tig.2014.01.001.
G. Rizk, D. Lavenier, and R. Chikhi. Dsk : k-mer counting with very low memory usage.
Bioinformatics, 29(5) :652–653, Feb 2013. doi : 10.1093/bioinformatics/btt020. [lien].
R. M. Sherman and S. L. Salzberg. Pan-genomics in the human genome era. Nature Reviews
Genetics, 21(4) :243–254, feb 2020. doi : 10.1038/s41576-020-0210-7.
N. Spies, J. M. Zook, M. Salit, and A. Sidow. svviz : a read viewer for validating structural
variants. Bioinformatics, page btv478, aug 2015. doi : 10.1093/bioinformatics/btv478.
64
phasing of structural variation in patient genomes using nanopore sequencing. Nature
Communications, 8(1), nov 2017. doi : 10.1038/s41467-017-01343-4.
M. Wellenreuther, C. Mérot, E. Berdan, and L. Bernatchez. Going beyond SNPs : The role
of structural genomic variants in adaptive evolution and species diversification. Molecular
Ecology, 28(6) :1203–1209, mar 2019. doi : 10.1111/mec.15066.
K. Ye, M. H. Schulz, Q. Long, R. Apweiler, and Z. Ning. Pindel : a pattern growth approach
to detect break points of large deletions and medium sized insertions from paired-end short
reads. Bioinformatics, 25(21) :2865–2871, Nov 2009. doi : 10.1093/bioinformatics/btp394.
[lien].
65
Annexe A
A.1 Publication 1
MindTheGap : integrated detection and assembly of short and long inser-
tions
Guillaume Rizk, Anaïs Gouin, Rayan Chikhi, Claire Lemaitre.
Bioinformatics 2014 30(24) :3451-3457.
66
Vol. 30 no. 24 2014, pages 3451–3457
BIOINFORMATICS ORIGINAL PAPER doi:10.1093/bioinformatics/btu545
insertions by combining paired-end mapping with read depth MINDTHEGAP on an actual whole-genome human dataset,
analysis. Finally, several methods detect sites of mobile element which required only 14 GB of memory.
insertions using collections of known transposable element se-
quences, by searching for read pairs where one mate is mapped
to a known element and the other to a unique part of the refer- 2 METHODS
ence genome (Ewing and Kazazian, 2011; Hormozdiari et al., The input of MINDTHEGAP is a set of reads and a reference genome. The
2010; Stewart et al., 2011). software performs three steps: (i) construction of the de Bruijn graph of
the reads, (ii) detection of insertion breakpoints on the reference genome
(find module) and (iii) local assembly of inserted sequences (fill module).
1.3 Reconstruction of inserted sequences Both the detection step and the assembly step rely solely on the con-
While short insertions are easy to reconstruct (as seen in Section structed graph.
1.2), to the best of our knowledge, only a few methods are cap- The output of the second step is a set of putative insertion positions on
able of handling medium or long insertions. They are based on the reference genome, whereas the output of the last step is, for each
global or local de novo assembly of reads that are potentially insertion site, one or several assembled sequences.
involved in an insertion.
SOAPindel (Li et al., 2013), Scalpel (Narzisi et al., 2013) and 2.1 de Bruijn graph construction
3452
MindTheGap
Heterozygous SNPs and deletions yield similar patterns, but the left
and right flanking k-mers are further separated from each other (k + 1
nucleotides apart for SNPs and 1-bp deletions, k+d 1 nucleotides
apart for deletions of size d). However, heterozygous inversions and
translocations do exhibit identical patterns. Also, inexact repetitions in
the reference genome create branching k-mers, which may yield by chance
the same pattern as a heterozygous insertion. To reduce this effect, we
apply an additional filter: the k – 1 suffix (resp. prefix) of the right (resp.
left) flanking k-mer must have less than h occurrences in the reference
genome. When h is set to 1, this prevents the detection of patterns that
may be generated by repetitions in the reference genome alone, in absence
of any sequence variants.
3453
G.Rizk et al.
3 RESULTS
3.1 Results on simulated datasets
MINDTHEGAP was applied on several simulated datasets to pre-
cisely estimate its recall and precision. This enabled to quantify
the impact of different levels of genome complexity, to independ-
ently evaluate each module and modes (detection versus assem-
bly, homozygous versus heterozygous) and to analyze the range
of insertion sizes MINDTHEGAP is able to detect and assemble.
Fig. 2. Fill module. A graph of contig is constructed from the left flank-
3.1.1 High recall and precision in homozygous mode For inser-
ing kmer L, in a breadth-first search order. Construction stops when a
tions of 1 kb, MINDTHEGAP recovered between 65 and 98.4% of
maximum number of nodes is reached, or when a branch becomes too
deep. The right flanking kmer R is searched within all nodes, finally all the simulated insertions, depending mainly on the complexity of
paths (in blue) between L and R are outputted as putative insertions the studied genome (Table 1). Almost all predicted homozygous
insertions are true-positive results, resulting in high precision
(consistently above 97%). Table 1 shows that almost all insertion
3454
MindTheGap
Table 1. Precision and recall results for MINDTHEGAP in homozygous mode on simulated and real datasets
Dataset Recall (%) Precision (%) N sim. Find module Fill module
TP FP TP FP
Note. Simulated insertions of size 1000 (homozygous). The number of deletions simulated in the reference genome appears in the column ‘N sim.’
size of detectable insertions, depending on the insert size of the 3.3 Application on real insertions of human
reads: given our simulation parameters, we observed that individual NA12878
SOAPindel recall decreased for insertions larger than 175 bp,
To evaluate the ability of MINDTHEGAP to recover real insertions
and the largest insertion detected was of length 189 bp.
in real data, we executed it on a human individual NA12878
Noticeably, the performance of SOAPindel was independent of
dataset containing 2.8 G 100 bp reads. As the coverage was
the genotype of insertions.
high, parameter k was empirically set to 63 and t to 5 (k-mers
with less than five occurrences were discarded). Predictions were
3.2 Evaluation on a real sequencing dataset of C.elegans then compared with a set of 74 large insertions predicted by
To evaluate the impact of real reads and a real donor genome alignment of fosmid sequences to the reference hg18 genome
with some degree of polymorphism in the reference genome, (see Section 2.4).
MINDTHEGAP was run on a C.elegans strain N2 read dataset 20 insertion sites were recovered by the find module. No het-
against the reference genome containing simulated deletions. erozygous insertions were predicted. We set r = 15, which
This is to simulate homozygous insertion variants in the donor enabled to find twice more sites than with r = 5. This suggests
genome. Additional insertions variants are likely to exist that real insertions contain longer repeated sequences at their
C.elegans strain. Thus, the number of FP could not be evaluated, breakpoints than expected in a random simulation. By analyzing
as the true set of insertions present in these reads is unknown. paired-end reads that mapped near each fosmid-predicted break-
For 1 kb insertion variants, 81.1% were correctly predicted point, we could infer the genotypes: only 23 breakpoints could be
and assembled by MINDTHEGAP (Table 1). Compared with the confidently assigned to a homozygous genotype (i.e. with less
fully simulated dataset on the same simulated insertions, the find than five read pairs spanning the breakpoint). The find module
module missed more insertion sites, whereas the fill module had a recovered 11 of them. Of the remaining 12 likely homozygous
better recall of inserted sequences. The first observation could be sites, the breakpoints of 8 of them were included in a large gap
explained by small polymorphism near the insertion breakpoints (k) in the reference binary string. This suggests that these sites
that generated longer gaps (see Section 2), whereas the second by were close to other form of polymorphism, which would explain
a higher read coverage in this dataset. why MINDTHEGAP did not detect them.
Additionally, we compared MINDTHEGAP and SOAPindel on Among the 20 detected insertions by the find module, the fill
this dataset with 1–100 bp simulated insertions. Recall values module succeeded in reconstructing correctly two inserted
were similar for both tools: 89% and 91%, respectively. sequences of sizes 4137 bp and 6729 bp, with respectively
3455
G.Rizk et al.
Table 2. Precision and recall results for MINDTHEGAP in heterozygous mode on simulated datasets, containing each 1000 simulated heterozygous
insertions of size 1000 bp
Dataset Recall (%) Precision (%) N sim. Find module Fill module k = 51 Fill module k = 31
TP FP TP FP TP FP
Note. Parameter r was set to 2, and assembled insertions smaller than 5 bp were filtered out.
4 DISCUSSION
MINDTHEGAP is the first integrated method to detect and assem-
3456
MindTheGap
for mobile element insertions. However, there exist methods tai- Chen,K. et al. (2014) Tigra: a targeted iterative graph routing assembler for break-
point assembly. Genome Res., 24, 310–317.
lored to the assembly of MEI, based on local assembly with
Chikhi,R. and Rizk,G. (2013) Space-efficient and exact de bruijn graph representa-
recruitment of mate reads. tion based on a bloom filter. Algorithms Mol. Biol., 8, 22.
Our tests on the NA12878 dataset showed there is room for DePristo,M.A. et al. (2011) A framework for variation discovery and genotyping
improvement: only two long homozygous insertions were suc- using next-generation dna sequencing data. Nat. Genet., 43, 491–498.
cessfully assembled out of 23 predicted ones. We postulate that Ewing,A.D. and Kazazian,H.H. Jr. (2011) Whole-genome resequencing allows
detection of many rare line-1 insertion alleles in humans. Genome Res., 21,
(i) polymorphism or repetitions near the insertion sites hinder
985–990.
detection by the find module, and (ii) the complexity of the Hajirasouliha,I. et al. (2010) Detection and characterization of novel sequence
human genome makes de novo assembly of large contigs difficult. insertions using paired-end next-generation sequencing. Bioinformatics, 26,
As no other tool was able to assemble long insertions, we could 1277–1283.
not assess whether our results were owing to weaknesses in our Hormozdiari,F. et al. (2010) Next-generation variationhunter: combinatorial algo-
rithms for transposon insertion discovery. Bioinformatics, 26, i350–i357.
method, or to specificities of this particular dataset (complex
Iqbal,Z. et al. (2012) De novo assembly and genotyping of variants using colored de
insertion sequences or mispredicted insertions). bruijn graphs. Nat. Genet., 44, 226–232.
Kidd,J.M. et al. (2010) A human genome structural variation sequencing resource
reveals insights into mutational mechanisms. Cell, 143, 837–847.
ACKNOWLEDGEMENTS Kim,S. et al. (2013) Reprever: resolving low-copy duplicated sequences using tem-
3457
A.2 Publication 2
Towards a better understanding of the low recall of insertion variants with
short-read based variant callers
Wesley Delage, Julien Thevenon, Claire Lemaitre.
BMC Genomics 2020, 21(1) :762.
74
Delage et al. BMC Genomics (2020) 21:762
https://doi.org/10.1186/s12864-020-07125-5
Abstract
Background: Since 2009, numerous tools have been developed to detect structural variants using short read
technologies. Insertions >50 bp are one of the hardest type to discover and are drastically underrepresented in gold
standard variant callsets. The advent of long read technologies has completely changed the situation. In 2019, two
independent cross technologies studies have published the most complete variant callsets with sequence resolved
insertions in human individuals. Among the reported insertions, only 17 to 28% could be discovered with short-read
based tools.
Results: In this work, we performed an in-depth analysis of these unprecedented insertion callsets in order to
investigate the causes of such failures. We have first established a precise classification of insertion variants according
to four layers of characterization: the nature and size of the inserted sequence, the genomic context of the insertion
site and the breakpoint junction complexity. Because these levels are intertwined, we then used simulations to
characterize the impact of each complexity factor on the recall of several structural variant callers. We showed that
most reported insertions exhibited characteristics that may interfere with their discovery: 63% were tandem repeat
expansions, 38% contained homology larger than 10 bp within their breakpoint junctions and 70% were located in
simple repeats. Consequently, the recall of short-read based variant callers was significantly lower for such insertions
(6% for tandem repeats vs 56% for mobile element insertions). Simulations showed that the most impacting factor
was the insertion type rather than the genomic context, with various difficulties being handled differently among the
tested structural variant callers, and they highlighted the lack of sequence resolution for most insertion calls.
Conclusions: Our results explain the low recall by pointing out several difficulty factors among the observed
insertion features and provide avenues for improving SV caller algorithms and their combinations.
Keywords: Short reads, Variant calling, Structural variants; Insertions
© The Author(s). 2020 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License,
which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate
credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were
made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless
indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your
intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly
from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative
Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made
available in this article, unless otherwise stated in a credit line to the data.
Delage et al. BMC Genomics (2020) 21:762 Page 2 of 17
The classical approach to discover SVs from Whole (MMBIR), generate such homologies whose size depend
Genome Sequencing (WGS) with short reads relies on precisely on the type of the involved mechanism. These
a first step consisting in mapping the reads to a refer- homologies can have an impact on insertion calling per-
ence genome. Then SV callers look for atypical mapping formance, since the concerned region at the inserted site is
signals, such as discordant read pairs, clipped reads or no longer specific to the reference allele and it is no longer
abnormal read depth, to identify putative SV breakpoints possible to identify the exact location of the insertion site.
along the reference genome [5, 6]. More than 70 SV callers However, little is known at present about the prevalence
have been developed up to date and several benchmarks of these homologies and their sizes for human insertion
have revealed great variability between results obtained variants due to their poor referencing in databases.
by different methods, demonstrating that SV detection More recently, novel long reads sequencing technolo-
using short read sequencing remains challenging [7, 8]. gies have overcome these limitations and allowed the
The challenge is to resolve two issues: a technical and generation of more accurate datasets, finally referenc-
a methodological one. The technical issue concerns the ing sequence-resolved insertion variants in the human
sequencing technology: insert size, read size and sequenc- genome [8, 17]. Thanks to several international efforts,
ing coverage have been shown to impact SV discovery. some gold standard callsets have been produced in 2019,
The second issue concerns SV caller algorithms and their referencing tens of thousands of insertions in several
ability to decipher and translate the biological signal from human individuals [18, 19]. Among the reported inser-
the alignments. Thus, SVs located in repeated regions or tions by Chaisson et al, a great majority (83%) could
containing repeats larger than the read size are difficult to not be discovered by any of the tested short-read based
detect [9]. SV callers. This result of recall below 17% is drastically
In particular, insertions are one of the most difficult different from the announced performances of insertion
SV types to call [7, 8]. Because the inserted sequence callers when evaluated on simulated datasets [20]. Indeed,
is absent from the reference genome, or at least at Chaisson et al showed that 59% of insertion variants
the given locus of insertion, calling such variants and were found in a tandem repeat context, suggesting that
resolving the exact inserted sequence requires finely most of the real insertion variants in human individuals
tuned approaches such as de novo or local assembly are probably occurring in complex regions and involving
[10, 11]. This increased difficulty is well exemplified by complex sequences. So far, such complexity factors were
the dramatic under-representation of such SV type in rarely included nor analysed in method benchmarks and
usual reference databases or standard variant callsets. to do so, actual insertion variants require to be better
For instance, dbVar at present references only 28% of characterized.
insertions or duplications among the SVs larger than Numerous countries are developing genomic medicine
50 bp. On the opposite, deletions represent more than programs, based on short-read sequencing. Although
70% of the database, although both types are expected third generation sequencing offers an unprecedented
to be roughly equally abundant in human populations technique for exploring the complexity of individual struc-
[12]. Moreover, only 1.5% of the reported insertions are tural variants, most of the genomic sequencing facilities
sequence-resolved, that is with an inserted sequence fully will still use short-read based sequencing in coming years
characterized. for its reduced cost. Hence, there is a critical need to mea-
One explanation is that the size of the reads is small sure and control the caveats of standard procedures for
compared to the target event size and the detection is detecting SVs with short-read sequencing data.
mainly based on alignments which may produce artefacts In this work, we performed an in-depth analysis of these
[13]. Another source of difficulty for insertion detec- unprecedented insertion callsets, in order to investigate
tion is the presence of repeated patterns at the precise the causes of short read based caller failures. We have
rearrangement breakpoints. Several molecular mecha- first established a precise classification of insertion vari-
nisms involved in rearrangement genesis are known to ants according to four different layers of characterization:
produce such repeated sequences, referred as junctional the nature and size of the inserted sequence, the genomic
homology [14–16]. Junctional homology is defined as a context of the insertion site and the breakpoint junc-
DNA sequence that has two identical or nearly identi- tion complexity. Because these levels are intertwined, we
cal copies at the junctions of the two genomic segments then used simulations to characterize the impact of each
involved in the rearrangement, when the sequence is complexity factor on the recall of several SV callers.
short (<70 bp) this is often called a micro-homology
[16]. The repair of DNA double strand breaks by diverse Results
mechanisms, such as Non-Allelic Homologous Recombi- In-depth analysis of an exhaustive insertion variant callset
nation (NAHR), Non-Homologous End Joining (NHEJ) In this work, we first aimed at precisely characterizing
or Microhomology-Mediated Break-Induced Replication an exhaustive set of insertion variants present in a given
Delage et al. BMC Genomics (2020) 21:762 Page 3 of 17
human individual. We based our study on a recently pub- else in the genome. Among tandem duplications, tandem
lished SV callset published by Chaisson and colleagues in repeats are characterized by multiple tandem repetitions
2019 [18]. Using extensive sequencing datasets, combin- of a seed motif within the inserted sequence. Mobile
ing different sequencing technologies and methodologi- element insertions are a very specific sub-type whose
cal approaches (short, linked and long reads, mapping- sequences are known and referenced in families. They
based and assembly-based SV calling), three human trios are notably characterized by very high copy numbers in
were thoroughly analysed to establish exhaustive and gold the genome (typically greater than 500). Other dispersed
standard SV callsets (Supplementary Table S1). We first duplication types were then required to have a copy num-
focused our study on the individual NA19240, son of the ber lower than 50, in order not to be confounded with
so-called Yoruban (YRI) Nigerian trio, whose SV callset potential mobile element insertions. We did not define
contained 15,693 insertions greater than 50 bp. segmental duplications and CNVs as additional sub-types
of dispersed duplications, as they are defined in the lit-
Nature and size of the inserted sequences erature by their size (above 1 Kb), the size being another
Insertion variants can be classified in different sub- independant level of characterization.
types according to the nature of the inserted sequence. In order to classify the insertion callset, all inserted
Three insertion categories were distinguished in the orig- sequences were aligned against the human reference
inal publication, namely tandem repeats, mobile element genome, a mobile element database and were scanned
insertions and complex ones for all the other types. We for tandem repeats (see Methods). We used a minimal
proposed to refine this classification in five insertion sub- sequence coverage threshold to annotate each insertion
types, illustrated in Fig. 1. A classical subdivision con- to an insertion sub-type according to the decision tree
sisted in distinguishing novel sequence insertions from described in Fig. 1. Insertions that did not meet any
insertions of exiting sequences, namely duplicative inser- requirement to be annotated as one of the previous sub-
tions. Several sub-types of duplicative insertions were types were qualified as unassigned insertions.
then defined according to the location or amount of the We set the threshold to 80% for our analysis to ensure
inserted sequence copies in the reference genome. Among a compromise between specificity and quantity of anno-
duplicative insertions, we proposed to stratify (i) tan- tated insertions in all sub-types. With such threshold, 88%
dem duplications, with at least one copy of the inserted of insertions could be assigned to a given type. Among the
sequence being adjacent to the insertion site, (ii) dispersed 13,850 annotated insertions, 8,735 (63%) were annotated
duplications, with copies that can be located anywhere as tandem repeats, 2,473 (17%) as mobile elements, 1,000
Fig. 1 Decision tree used to classify insertion variants. Five insertion sub-types are defined according to the nature of the inserted sequence : novel
sequence, tandem repeat and mobile element insertions and tandem and dispersed duplications. Unassigned insertions refer to insertions which
do not meet the requirements to be assigned to at least one sub-type
Delage et al. BMC Genomics (2020) 21:762 Page 4 of 17
(7%) as tandem duplications, 869 (6%) as novel sequences distributions differed between insertion types (Fig. 3a).
and 773 (5%) as dispersed duplications (Fig. 2b and Mobile elements showed the most contrasting size distri-
Supplementary Table S2 for results obtained with other bution with a strong over-representation of the 250-500
coverage thresholds). 46% of tandem repeats had a repeat bp size class (61%). This can be explained by the most
seed smaller than 10 bp and 93% smaller than 50 bp. Com- frequent and active mobile element class in the human
pared to the classification of Chaisson et al, the propor- genome being the SINE elements of size around 300
tions of tandem repeats (57% vs 56%) and mobile elements bp. Notably, the novel sequence insertion type carried a
(23% vs 16%) were close. The difference in mobile ele- greater proportion of large insertions than other types,
ment proportions mainly represented insertions that were with 164 (19%) of the 869 novel sequences larger than
unassigned in our annotation. The 1,843 (12%) unassigned 1,000 bp.
insertions at 80% threshold showed partial annotations
of mobile element (57%), tandem repeats (22%), tandem Characterization of insertion locations in the genome
duplications (15%) or dispersed duplications (5%). We then characterized the insertions based on the
Concerning the size of the insertions, 67% of the inser- genomic context of their insertion site. We investi-
tions were smaller than 250 bp and only 8% had a gated in particular genomic features that might make
size greater than 1 Kb (Fig. 2a). Interestingly, the size read mapping and SV calling difficult, such as the
Fig. 2 Distributions of insertion variant features across several callsets. Distributions of a insertion size, b insertion type, c repeated context of
insertion and d homology size at the breakpoint for NA19240, HG00514, HG00733 and HG002 insertion variant callsets. Abbreviations: SimpleRep for
simple repeat, ME for mobile element, TandemRep for tandem repeat, TandemDup for tandem duplication, DispersDup for dispersed duplication
Delage et al. BMC Genomics (2020) 21:762 Page 5 of 17
Fig. 3 Proportions of insertion variant features according to the type of insertion. Proportions of classes of a insertion size, b insertion location, c
homology size at the breakpoint according to the type of insertion in the NA19240 callset. Abbreviations: SimpleRep for simple repeat, ME for
mobile element, TandemRep for tandem repeat, TandemDup for tandem duplication, DispersDup for dispersed duplication
repetitive content. A strong over-representation was is an amplification of a seed in the reference genome.
found in regions annotated as simple repeats, with 9,675 Thus the largest homology size corresponded to the seed
(70%) of the annotated insertions located in these regions size presents at the right breakpoint (in case of left
that only represent 1.2% of the genome (Fig. 2c). The pre- normalization). As for tandem duplications, the discor-
ferred genomic context of insertions varied between inser- dance between their annotation as tandem duplication
tion types (Fig. 3b). 8,047 (92%) tandem repeats, 723 (73%) and the smaller size of the detected junctional homol-
tandem duplications and 519 (63%) dispersed duplica- ogy is related to the difference in the methods used to
tions were found in simple repeat regions. Conversely, 580 define the homology, where small distances (<10 bp) to
(67%) novel sequence insertions and 1,383 (56%) mobile the insertion site and to the inserted sequence extremity
element insertions were located in other regions. We did were required in the junctional homology case, whereas
not find a higher rate of insertions within exonic, intronic in the tandem duplication annotation case, the homolo-
or intergenic regions compared to a uniform distribution gous segment had only to cover at least 80% of the inserted
along the genome. sequence.
different variants. Using a rough estimation of shared variant call as annotated in the callsets. For the individual
variants, we identified only 1,169 insertion sites com- NA19240, 2,363 (17%) insertion variants were comforted
mon to the four callsets within a 1 kb size window. by SR-based SV callers. As shown in Fig. 4, this SR-
On average 3,344 insertions were shared between two based recall was highly heterogeneous with respect to the
given callsets, and overall, more than 55% of the studied previously described insertion features. Each described
insertions were specific to a given callset. The distribu- feature in this work (ie. nature and size of the inserted
tions of insertion types, sizes, locations and junctional sequence, insertion site genomic context and junctional
homology sizes were similar between the three individ- homologies) impacted the SR-based recall. As shown in
uals of the Chaisson et al study and the GiaB callset Fig. 4a, insertions larger than 500 bp were poorly discov-
(Fig. 2). ered by SR-based methods (<3%). An increased SR-based
recall for the 250-500 bp insertion size class corresponded
Short-read-based recall to mobile element insertions. The greatest difference
In order to investigate whether the previously described in SR-based recall was observed among the insertion
insertion features impacted the recall of short-read- types: 1,410 (56%) mobile elements and 342 (40%) novel
based (SR-based) SV callers, we reproduced our previous sequence insertions could be detected with SR-based SV
analysis according to the technology involved in the callers compared to only 87 (9%) tandem duplications, 484
Fig. 4 Proportions of SR-based insertion discoveries according to insertion features. Proportions of insertions that were called using short read
technology data according to a insertion size, b insertion type, c insertion location and d homology size at the breakpoint, in the NA19240 and
HG002 callsets. These callsets were provided by two different studies using different discovery tools and methodologies
Delage et al. BMC Genomics (2020) 21:762 Page 7 of 17
(6%) tandem repeats and 40 (5%) dispersed duplications or inserted sequence. As a more stringent evaluation, the
(Fig. 4b). sequence-resolved recall considered as true positives only
The variations of the SR-based recall with respect those insertion calls having a correct genomic position
to insertion features were very similar between the and whose inserted sequence was very similar to the sim-
three studied individuals from the Chaisson et al. study ulated one (>90% sequence identity and +/- 10% insertion
(Supplementary Figure S2). However, the same compar- size).
ison across two different studies with different method-
ologies was much more contrasted. Firstly, overall around Factors impacting insertion site recall
1.6 times more insertions in proportion could be detected Recalls of insertion sites for all four methods are pre-
by SR-based methods in the GiaB study compared to the sented for the different simulated datasets in Table 1.
Chaisson et al study (SR-based recalls of 28% and 17% On the baseline simulation, all tools succeeded to detect
for HG002 and NA19240 callsets respectively). Secondly, 100% of simulated insertions, except for GRIDSS with
the SR-based recall was more homogeneous with respect 81% of recall. The size of the inserted sequence impacted
to insertion features in the GiaB callset (Fig. 4). The fea- the recall of the insertion sites for most tools, except
ture showing the most impact was the insertion size with MindTheGap. GRIDSS was challenged by small insertions
a decrease of the SR-based recall with the insertion size, (50 bp) whereas Manta and SvABA had more issues with
reaching below 5% for insertions larger than 500 pb for large insertions. The most extreme behavior was observed
both studies (Fig. 4a). Similarly to the NA19240 callset, for SvABA which was not able to find the insertion sites of
tandem repeats appeared more difficult to discover with any of the simulated novel sequences larger than 500 bp.
SR-based methods, but to a lesser extent in the GiaB When simulating various insertion types, GRIDSS was
callset (Fig. 4b). Insertions located in simple repeats were the only tool whose recall was not negatively impacted.
less discovered using SR-based methods but this SR-based Manta could not find any type of dispersed duplica-
recall of 25% remained higher than for NA19240 where tions and showed a lower recall to detect tandem repeats
it only reached 5% on these locations (Fig. 4c). Junctional with 25 bp size seeds. MindTheGap was unable to detect
homology of the insertions of individual HG002 did influ- any type of tandem duplications and found only 58%
ence its SR-based recall, but in a different manner than in of mobile element insertions. SvABA was not able to
the Chaisson et al study (Fig. 4d). detect any tandem repeat insertion but was able to
detect all dispersed and tandem duplications and mobile
Using simulations to investigate the factors impacting the elements.
insertion calling recall Concerning junctional homology, the tools showed con-
In real insertion callsets, most of the previously identified trasting behaviors. GRIDSS was the only tool unaffected
factors impacting SV discovery are intertwined. In order by the presence and size of repeated sequence at the inser-
to quantify the impact of each factor independently, we tion junctions. On the contrary, MindTheGap was the
produced various simulated datasets of 2x150 bp reads most impacted by junctional homology, being unable to
at 40x coverage, containing each 200 homozygous inser- detect insertions with homology at any tested size. This
tion variants on the human chromosome 3. As a baseline, feature is actually controlled by a parameter of MindThe-
we simulated 250 bp novel sequences taken from Sac- Gap, increasing the max-repeat parameter value to 15 bp
charomyces cerevisiae exonic sequences inserted inside (default : 5bp), MindTheGap discovered 99% of the inser-
human exons. This is meant to represent the easiest type tion sites with 10 bp junctional homomolgies. Manta’s
of insertions to detect. Then, we considered four scenar- recall decreased with the size of junctional homologies,
ios of simulations, where only one of the four previously whereas SvABA handled small (less than 20 bp) or very
studied factors is changed at a time with respect to the large (150 bp) junctional homologies but was affected by
baseline simulation. medium sizes.
Four insertion variant callers were evaluated on these Concerning the impact of the genomic context of inser-
datasets. They were chosen according to their good per- tions, no loss of recall was observed in non repeated
formances in recent benchmarks [7] and to maximise the locations. Alignment-based SV callers showed no change
methodological diversity. GRIDSS [11], Manta [20] and in recall in small simple repeat (<300 bp), SINE and LINE
SvABA [6] are based on a first mapping step to the ref- locations. Manta and SvABA recalls lost 5 to 6% of recall
erence genome, contrary to MindTheGap [10] which uses in simple repeat regions larger than the insert size (>300
solely an assembly data structure (the De Bruijn graph). bp). MindTheGap lost 42 and 47% of recall in large simple
Two types of recall were computed depending on the pre- repeat and SINE location simulations. Simulating inser-
cision and information given for each call: insertion-site tions close to each other on the genome, at less than 150
only recall only evaluated if an insertion was called at an bp, reduced the recall of SvABA (-98%), MindTheGap
expected genomic position regardless of the predicted size (-33%) and Manta (-15%).
Delage et al. BMC Genomics (2020) 21:762 Page 8 of 17
Table 1 Insertion site recall of several short-read insertion callers according to different simulation scenarios. Cells of the table are
colored according to the variation of the recall value of the given tool with respect to the recall obtained with the baseline simulation
(first line, colored in blue): cells in red show a loss of recall >10%, cells in grey show no difference compared to baseline recall at +/- 10%
Insertion site only recall (%)
GRIDSS Manta SvABA MindTheGap
Baseline simulation: 250 bp novel seq. in exons 83 100 100 100
50 bp 56 100 100 100
Scenario 1
Insertion 500 bp 100 86 0 99
size
1,000 bp 100 88 0 98
Dispersed duplication 100 1 100 96
Tandem duplication 100 100 100 0
Scenario 2
Insertion Mobile element 100 2 100 58
type
Tandem repeat (6 bp pattern) 100 90 1 0
Tandem repeat (25 bp pattern) 99 66 0 2
10 bp 100 100 96 0
20 bp 100 100 85 0
Scenario 3
Junctional 50 bp 77 68 12 0
homology
100 bp 100 22 49 0
150 bp 100 0 100 0
Non repeat 83 100 99 96
Simple repeat (<300 bp) 82 100 100 73
Simple repeat (>300 bp) 87 94 95 58
Scenario 4
Genomic SINE 90 100 99 53
location
LINE 80 100 97 90
Clustered insertions (<150 bp) 85 85 2 77
Novel sequences at real locations 84 80 71 38
Scenario 5
Real insertions in exonic regions 84 74 57 24
Real insertions
Real insertions at real locations 39 35 44 6
Finally, when simulating the 889 insertions of NA19240 Removing this quality filtering allowed to increase the
callset located on chromosome 3, with their reported recall mainly for GRIDSS and SvABA (see Supplemen-
inserted sequence at their real locations as described in tary Table S3). Remarkably, GRIDSS reached a 100% recall
the variant calling file (scenario 5), the recall of all tools on almost every scenario, except the scenario simulat-
dropped to less than 44%, reaching for many tools their ing the real insertions where still a 35% loss of recall
lowest values among the different simulated datasets. This was observed (Supplementary Table S3). These differ-
was particularly marked for GRIDSS whose recall was ences indicated that a substantial amount of true posi-
greater than 77% in all simulated scenarios, but achieved tive insertions were detected but reported as low quality
only 39% on this simulation. When relaxing one complex- calls.
ity factor, the type or the location, ie. simulating either
Sequence-resolution of predicted insertions
novel sequences at the real locations or the real types in
We then investigated whether the SV callers were also
exonic regions, the drop of recall is much smaller for all
able to recover the full inserted sequences in the dif-
tools, indicating that there is a synergetic effect of combin-
ferent simulation scenarios (Table 2). On the base-
ing in a single insertion event these two factors, insertion
line simulation with 250 bp novel sequence insertions,
type and insertion location.
every tools reported for almost all detected insertion
Impact of quality filtering sites a resolved and correct inserted sequence. However,
Previous results were computed using only the calls these high sequence-resolved recalls dropped dramati-
assessed with sufficient quality by each tool and anno- cally when deviating from the baseline scenario. Although
tated as PASS in the FILTER field of the VCF file. the discovery of insertion sites was not much impacted by
Delage et al. BMC Genomics (2020) 21:762 Page 9 of 17
Table 2 Sequence-resolved recall of several short-read insertion callers according to different simulation scenarios. Cells of the table are
colored according to the variation of the recall value of the given tool with respect to the recall obtained with the baseline simulation
(first line, colored in blue): cells in red show a loss of recall >10%, cells in grey show no difference compared to baseline recall at +/- 10%
Sequence-resolved recall (%)
GRIDSS Manta SvABA MindTheGap
Baseline simulation: 250 bp novel seq. in exons 81 100 96 100
50 bp 56 100 100 100
Scenario 1
Insertion 500 bp 0 0 0 99
size
1,000 bp 0 0 0 98
Dispersed duplication 0 0 16 96
Tandem duplication 0 0 0 0
Scenario 2
Insertion Mobile element 0 0 61 58
type
Tandem repeat (6 bp pattern) 0 0 1 0
Tandem repeat (25 bp pattern) 0 0 0 0
10 bp 99 100 92 0
20 bp 100 100 78 0
Scenario 3
Junctional 50 bp 6 46 10 0
homology
100 bp 0 11 0 0
150 bp 0 0 0 0
Non repeat 80 99 98 96
Simple repeat (<300 bp) 77 98 97 73
Scenario 4 Simple repeat (>300 bp) 77 93 90 58
Genomic
location SINE 77 99 94 53
LINE 76 97 95 89
Clustered insertions (<150 bp) 75 73 2 77
Novel sequences at real locations 64 73 67 37
Scenario 5
Real Real insertions in exonic regions 11 14 14 9
insertions
Real insertions at real locations 6 23 30 6
the insertion size, all tools but MindTheGap were not able False positive amount variations
to recover any of the inserted sequences when it was larger The tools with the largest recalls were also the tools pro-
than 500 bp (Table 2). On the contrary, MindTheGap ducing the largest amounts of false positive discoveries (in
assembled correctly nearly all simulated novel sequences, the order of several hundreds for GRIDSS and SvABA, see
even those of 1 Kb. Concerning the other insertion types, Supplementary Table S4). More surprisingly, the amount
tools were not able to provide sequence resolved calls, of false positives was not constant for most tools between
except for MindTheGap and SvABA for some dispersed the different simulation scenarios. It increased when sim-
duplications and mobile element insertions (Table 2). In ulated insertions presented a duplicative pattern (mobile
the case of tandem repeats, GRIDSS which detected all element, dispersed duplication and junctional homologies
insertion sites, reported inserted sequences of at most above 50 bp). For those, some SV callers predicted vari-
150 bp (instead of 250), corresponding to the simu- ants not only at the insertion site but also at the locations
lated read size. The increase of junctional homology of homologous copies of the inserted sequences. Remov-
size reduced the sequence resolution of GRIDSS and ing the quality filter led to a large increase of the amount
SvABA. Insertions located in repeated regions were less of false positive discoveries for GRIDSS and SvABA (5 to
resolved than in the baseline simulation for every tools. 17 times more respectively).
Finally, the sequence resolution of real insertions simu-
lated at their real locations decreased compared to the Unions and intersections of SV callers
insertion site recall, GRIDSS suffering the greatest loss A classical strategy to report SVs on real data is to recon-
(-33%). cile several SV callsets keeping only variants that are sim-
Delage et al. BMC Genomics (2020) 21:762 Page 10 of 17
ilarly called by different SV callers. This strategy ensures a 3 insertions at their real locations. Concerning sequence
balance between true and false disovery rate. On the last resolution, although Sniffles calls contained systemati-
simulation scenario, only 12% of the insertion sites were cally a full inserted sequence, the latter was imprecise and
validated by the three tools, GRIDSS, Manta and SvABA, contained sequencing errors leading to sequence-resolved
and 39% by at least two tools. However, the union of all recalls around only 20% when requiring at least 90% of
three methods comprised 65% of the real insertion sites, sequence identity. When relaxing the identity threshold
which represented an increase of 20% of the best recall to 80% or using the dedicated benchmark tool SVanalyzer
obtained by a single method (Fig. 5). from GiaB which relies on a less stringent validation, the
sequence-resolved recall was similar to the insertion site
Evaluation of insertion recall with long read simulated data recall for most insertion scenarios (Supplementary Tables
For each of these short-read simulated datasets, we also S5 and S6). These results reveal that long read technolo-
simulated a corresponding PacBio long read dataset, with gies enable the discovery of every types of insertion but
40 X coverage and 16% error rate. We then applied a the calls remain imprecise.
state-of-the-art long-read SV caller, Sniffles [17], on each
of them to assess whether the previously identified dif- Discussion
ficulty factors for short read data have also an impact The discovery of genomic variants is an important
on the recall with long read data (see Supplementary step towards the understanding of genetic diseases and
Table S5). For most insertion scenarios, Sniffles reported species evolution [21, 22]. The detection of insertions too
accurately 100% of the insertions sites, except for the tan- small (<1kb) to be detected using comparative genomic
dem duplication type and for the insertions with large hybridization array (CGH array) but larger than indel
junctional homologies (recall below 20%). In these cases, size (>50 bp) to be detected by the gold standard small
insertions were in fact reported but at more than 10 bp variant discovery pipeline (GATK), remained a challenge
from the simulated insertion site. This is probably due with short read technology [4]. Thus these variations were
to imprecise sequence resolution preventing the correct poorly characterised in databases as compared to other
left normalization of breakpoint positions. Another diffi- SVs such as deletions. Numerous variant callers have been
culty factor was the close proximity of insertion locations, developed to overcome this issue but without resolving
for which Sniffles reported one complex event instead of it [7]. Long read technologies or the crossing of vari-
several close insertions. This mainly explained the low ous sequencing technologies overcome these limitations
recall of 58% for the dataset with the real chromosome but are not affordable for many applications such as rou-
tine diagnosis of genetic diseases [18]. Thus, to improve
current and future SR based SV callers, a better under-
standing of the actual insertion variants present in human
populations is required.
We have presented here one of the most detailed and
comprehensive analyses of actual insertion variants in the
human genome looking for factors impacting their detec-
tion with short read re-sequencing data. This could be
possible thanks to the publication of two exceptional SV
callsets by Chaisson et al. [18] and Zook et al. (GiaB)
[19]. These catalogs of insertions are considered as the
most exhaustive for a given human individual and are
qualified as gold standards thanks to their extensive val-
idation by extensive and cross technology sequencing
datasets. Unlike in the Chaisson et al study, the GiaB
callset contained two categories of variants : 7640 inser-
tions that were reported with a higher confidence (PASS
in the FILTER field) and 6210 other insertions. As men-
tioned by the authors, the first category is likely to
Fig. 5 Intersections of true positive insertion callsets between
different SV callers. Intersections of true positive insertion callsets
be biased towards easier to discover variants. Because
between GRIDSS, SvABA and Manta on the scenario 5 simulation (real we did not want to introduce this potential bias, and
insertions at real locations). In this scenario, the 889 insertions located after checking that these two categories showed similar
on the chromosome 3 from the NA19240 callset were simulated as insertion feature distributions (see Supplementary Figure
described in the vcf file. Insertion calls were validated and compared S1), we decided to conduct our analyses on the whole
based solely on the insertion site prediction
callset.
Delage et al. BMC Genomics (2020) 21:762 Page 11 of 17
Not only, these catalogs of insertion variants are consid- results show that real insertion variants harbour substan-
ered as the most exhaustive for a given human individual, tially larger junctional homologies than insertions that
but they are also the first sets with sequence-resolved would be drawn randomly. Our measures allowed us to
events for any size and type of insertions. The fine resolu- compare such feature between insertion types and all
tion of the inserted sequences, present in these datasets, insertion types have been found to have a substantial
enabled us to propose a refined classification of inser- proportion of variants with large junctional homologies
tion variants. In the two datasets, insertion types were not (greater than 20 bp). Results showed also that large inser-
formally defined and the classifications differed between tions tended to carry larger junctional homologies. As
the datasets. Our classification allowed to normalize these expected by their tandem nature, tandem repeats and tan-
heterogeneous annotations and was a direct application dem duplications had larger homology sizes than other
of variant definitions from the dbVar database which is insertion types.
based on the sequence ontology (SO) [12]. We based our All the features of insertions characterized in our study
insertion type annotation on a minimal sequence cov- (ie. nature and size of the inserted sequence, insertion
erage threshold, that was set to a relatively high value, site genomic context and junctional homologies) showed
80%, in order to ensure a good specificity of our annota- to impact the ability of SR-based SV callers to discover
tion. Increasing this value led to many more unassigned these variants, as defined by method annotations in the SV
insertions, as the annotations were based on sequence callsets. However, an important difference was observed
alignments that were affected by potential remaining between the two studies, with the GiaB study being able
sequencing errors in the inserted sequences, polymor- to detect with short reads almost twice as many insertions
phism with the reference genome and the usage of align- in proportion than in the Chaisson et al study. The dif-
ment heuristics. If the amount of unassigned insertions ference in SR-based recalls between the two studies can
decreased with the coverage threshold value, proportions certainly be explained by the difference in the read depths
of the different insertion types remained quite stable of sequencing datasets (77X vs 300X for Chaisson et al and
(Supplementary Table S2). Among the 12% of unassigned GiaB studies respectively), by the different SR-based tool
insertions, some could correspond to a mixture of several sets used and by the different callset filtering and merg-
insertion types, which particular case was not considered ing methodologies. The two studies used roughly the same
in this study. number of SV-callers (13 and 15), but with a poor inter-
As previously reported in the Chaisson et al and GiaB section: only one SV-caller (Manta) was common to both
studies, we observed a highly heterogeneous distribution studies. Additionally, the method annotation of each vari-
of insertion types and locations along the genome. The ant is highly dependant on the study methodology to filter
vast majority of insertions consisted in tandem repeats and merge the numerous callsets obtained for the same
(63%) and most insertion sites were located in simple individual with different sequencing technologies and SV
repeat regions (70%). These regions of low complexity, callers. For instance, it is not clear if the presence of an SR-
although representing a small proportion of the genome based tag for a given variant does necessarily mean in both
(1.2%), are therefore a major source of inter-individual studies that the latter can be sequence-resolved solely
variability. using short reads. However, both studies showed simi-
The sequence-resolution provided in these SV callsets lar weaknesses to detect tandem repeats, large insertions
also enabled us to analyze precisely the breakpoint junc- and insertions located inside simple repeats. These obser-
tions of each insertion variant. Junctional homology has vations are in-line with the already known difficulties of
been shown to be a frequent feature of SVs, that can be mapping short reads in such contexts.
used to infer the rearrangement molecular mechanism These disparities between studies and the fact that most
[14, 15]. Although, it has been previously described for identified factors responsible of low SR-based recall are
human SV callsets (around 2,000 SV breakpoints, includ- intertwined with one another in real insertion variants led
ing less than 400 insertions) [15], this is, to our knowledge, us to pursue these investigations with simulated data. Our
the first exhaustive quantification of junctional homology simulations did not aim at providing an exhaustive bench-
for such a large and almost complete set of insertions in mark of SV callers but at identifying the precise genomic
a human individual. However, our measure of homology factors of insertion variants that prevent their correct dis-
size is highly dependent on the callset precision of the covery with short reads. As a consequence, we selected
insertion site location and of the inserted sequence. As a small but diverse set of SV callers and we deliberately
SVs are often difficult to precisely localize, are subject to ran them with their default parameters. We based our
left-normalization processes, and their inserted sequences selection of SV callers on a recent and comprehensive
were mostly obtained from error-prone long reads, our benchmark study by Kosugui et al. [7]. SV callers selected
measures may likely result in an under-estimation of the in our study were chosen for their good performance in
actual homology sizes. Despite these potential biases, our this benchmark, for their diversity of algorithms and for
Delage et al. BMC Genomics (2020) 21:762 Page 12 of 17
their ease of installation and usage. MindTheGap was not detectable but often not reported with a sufficient qual-
among the best insertion callers identified by Kosugui et al ity due to this lack of resolution. Furthermore, sequence
but was the only one not based on read mapping and using resolution is essential for the comparison and genotyp-
intensively de novo assembly with the whole read dataset. ing of SVs in many individuals. As these tasks are the
Simulations remain a powerful approach to identify the basis for association studies and medical diagnosis, efforts
strengths and weaknesses of SV callers but they were not should be directed towards a better resolution of the
meant to reflect perfectly real situations. In our simula- sequence of these variants [8, 23]. Results obtained with
tions, several features may be far from the real complexity the local assembly tool MindTheGap showed that the
of human genome re-sequencing, such as some sequenc- use of the whole read dataset allowed many insertions
ing technology biases, the use of one chromosome instead and even large ones to be assembled. The restriction
of the whole genome, and the absence of other polymor- to a small subset of reads to perform local assembly
phisms than insertion variants (SNPs, small indels and may therefore be the shortcoming of the other tested
other SVs). As a consequence, the reported recalls are SV callers. Resolving the inserted sequence is possible
likely to be over-estimations of the ones obtained with to some extent, but tandem repeats larger than the read
real data. Although absolute values should be interpreted size will remain difficult to resolve with short reads
with caution, they can readily be compared between SV technology.
callers and between simulation scenarios. As a matter of Interestingly, sequence resolution appeared also to be
fact, we often observed strong differences in recalls allow- an issue with long read sequencing data. In this case,
ing to provide interesting insights in terms of impacting the tested long read SV caller did report full inserted
factors and SV caller behaviors. Our simulation protocol sequences but with a poor sequence precision, due the
enabled to study each difficulty factor independently and higher sequencing error rate. This issue also prevented
highlighted the larger impact of insertion type compared the correct left normalization of insertion sites leading
to insertion location. However, all studied factors taken to erroneous insertion locations. This low accuracy of
independently could not explain the whole loss of recall predicted calls is likely to hamper the genotyping and
when simulating the real insertions at their real locations comparison of SV calls between individuals. Our results
and there is probably an important synergetic effect of therefore showed that there is also a need to improve long
combining in a single insertion event several of the stud- read SV callers as well.
ied factors. For instance, the discovery of novel sequences Overall, the different SV callers did not performed
in repeated regions was not a problem for almost every well in every situation and in every aspects of inser-
tested tools. However, the change of novel sequences to tion calling. Each caller showed its own strengths and
real inserted sequences, most of them corresponding to weaknesses, often different from the other tools. Precisely
tandem repeats, reduced by half the recall of SV callers. identifying these in terms of insertion variant features
Our simulations revealed that junctional homologies as and genomic contexts will enable each tool to be used to
small as 10-20 bp impacted the recall of all tested tools. its best advantage. To do so, benchmark studies should
Such repeated sequences are likely to alter the mapping take into account the wide variability of variant features
signature targeted by SV callers. Although such features of that this present work has highlighted. Two recent SV
SV breakpoints and their relation to the molecular mecha- benchmarks have raised awareness of the variability in
nisms generating SVs have long been described, they seem the performances of SV callers depending on data sets
to be rarely taken into account in the design of SV caller and approaches [7, 9]. They looked at several factors that
algorithms. Our study of the real insertions showed that could be responsible for this variability. Technical factors
such junctional homology sizes are relatively common, (reads size, insert size and sequencing coverage) and bio-
with almost 40% of insertions with junctional homologies logical factors (nearby SNVs or indels, genomic context,
larger than 10 bp. Therefore, SV callers algorithms would and variant size) showed to impact the recall of SV callers.
benefit from taking into account such properties of the However, the latter factors were analyzed for all SV types
breakpoints, that are likely to generate very specific signals combined and none of these studies took into account
in terms of read mapping. the different types of insertion variants. Best practices for
One striking result of our simulations is the absence of benchmarking small variant calling have been suggested
sequence resolution for most of the simulated insertion based on gold standard callsets in high confidence regions,
features and most of the tested SV callers. In addition to leaving structural variation in the fog [24]. However, it is
the obvious loss of information about the variant event, precisely this type of variation that requires best practices
this also limits the identification of the insertion type, the for benchmarking and a standardization of annotation as
genotyping and the validation of the predicted call. As a they are harder to identify and report. We hope that the
matter of fact, we observed that most insertions regard- present fine characterization of gold standard human SV
less of their type and insertion genomic context were callsets will help in the development of better practices
Delage et al. BMC Genomics (2020) 21:762 Page 13 of 17
for benchmarking SV callers, for both short and long read association studies to variants other than punctual ones,
sequencing data. allowing for instance the development of personalised
Advises to improve the detection using short read tech- medicine and the resolution of diagnostic bottlenecks for
nology have already been described such as the careful many rare diseases.
combination of complementary SV callers [7]. Meta SV
callers such as Meta-sv, Parliament2 or sv-callers recon- Methods
cile SV calls produced by different SV callers [25–27]. Data origin
However, only the calls that are discovered concordantly SV callsets from the Chaisson et al. study [18] were
between different tools are returned. This strategy allows obtained from dbVar with the accession nstd152. The
the precision to be increased, but at the expense of the HG002 SV callset, Tier 1 version v0.6, from the GiaB
recall. Our simulations showed that the intersection of study [19] was used (see the full ftp links in the Declara-
only three SV callers reduced the recall of 30%, whereas tions section). Only insertions from the core genome, that
taking their union could increase the recall by at least 20%. were larger than 50 bp and sequence resolved (ie. with
Considering unions of callsets would require a careful an inserted sequence entirely defined) and called also in
control of false positive rates. A better control could prob- at least one of the parents were kept. No filtering related
ably be achieved with sequence-resolved variants and by to quality or coverage was applied. In the HG002 callset,
taking into account the observed characteristics of the dif- insertion calls containing the “LongHomRef ” tag in the
ferent insertion types. Another alternative, less described, FILTER field were removed because they were not con-
could be the use of dedicated tools for each type of inser- firmed by long read genotyping methods and they had
tion, instead of using only general-purpose SV callers. thus a higher probability to be false positive discoveries
Among them, Expansion Hunter has been designed to (359 insertions). The human reference genome version for
detect tandem repeats, Pamir and Popins for novel inser- this study was Hg38 (GRCh38). To compare the callsets
tions and TARDIS for large duplications [28–31]. on the same reference genome, the HG002 callset pro-
duced on hs37d5 build was converted into Hg38 build
Conclusion using Picard, the hs37d5 to hg19 and the hg19 to hg38
In this work, we produced a detailed characterization chain files from GATK public chain files. Noteworthy, this
of the insertion variants in a given human individual. process can have some impacts on a few SV calls, since
We identified many factors of human insertion variants some genomic regions can differ between the reference
that explain their low recall with SR-based SV callers, versions. In particular, the conversion (liftover) induced a
including complex insertion types, difficult genomic con- loss of 60 SV calls.
texts, large insertion sizes and junctional homologies at
the breakpoints. The significant variability in the char- Comparison of the callsets
acteristics of the insertion variants, as well as the fact As a rough estimation of the amount of shared inser-
that all difficulties were handled differently by the dif- tion variants between callsets, insertion locations were
ferent tested SV callers, call for a better characterization compared regardless of the insertion type or sequence.
and comparison of SV callers according to the targeted Insertion variants located less than 1,000 bp apart from
variant features. The comparison results presented here one another were considered as the same variant.
already provide some concrete suggestions to improve
insertion variant calling with short reads. First, insertion Insertion type annotation
site detection could be improved by taking into account TandemRepeatFinder (TRF) was used to annotate tan-
the atypical mapping signals generated by large junctional dem repeats within each inserted sequence [32]. Recom-
homologies. Then, sequence-resolution recall could be mended parameters were used, except for the maximum
improved by using the whole read set instead of recruited expected TR length (-l) which was set to 6 millions. In
read subsets for the assembly of the inserted sequence. order to annotate mobile elements in inserted sequences,
Our simulation protocol also allowed us to identify com- we used dfam, one of the annotation tools of Repeat-
plementarities between different SV callers and showed Masker [33]. Each inserted sequence was scanned by dfam
that insertion recall could be significantly improved by with the standard HMM profile database of human mobile
taking the union of calls. Finally, based on these com- elements provided by the tool. For the annotation of dis-
plementarities and with improved sequence-resolution, persed duplications and the occurrence count of their
smarter consensus selections, than simply callset unions, copies in the reference genome, each inserted sequence
taking into account insertion type, size and context, could was locally aligned against the Hg38 genome using Blat
be designed to reach a high recall while controlling with default parameters [34]. Only the alignments with at
the False Discovery Rate. Such improvements are cru- least 90% of sequence identity were kept. For the annota-
cial for the generalization of population genomics and tion of tandem duplications, the two sequences on either
Delage et al. BMC Genomics (2020) 21:762 Page 14 of 17
side of the insertion site and of the same size as the these random insertions were identified using the same
insertion were aligned against the inserted sequence using previously described methodology as for real insertions.
Blat.
We used a minimal sequence coverage threshold, Genomic context characterization
Mincov , to annotate the insertions. To be assigned to a To study the genomic context of insertions, we used the
given sub-type, the inserted sequence had to contain at repeat content annotations of RepeatMasker from the
least one contiguous segment annotated with the cor- UCSC genome browser for the Hg38 genome and the gene
responding type and covering at least Mincov % of the annotations from the Gencode v32 [35–38]. Simple repeat
inserted sequence. Novel sequence insertions were a spe- location were extracted from the dedicated simple repeat
cial case where the contiguity of the annotation was not file from the UCSC genome browser.
required: more than Mincov % of the inserted sequence
should not be covered by any alignment with the reference SR-based recall of the gold standard callsets
genome nor with the mobile element reference sequences, Each callset was partitioned in two parts based on the
nor contain tandem repeats. When several types fulfilled discovery technology. The first part, referred as Short
the minimal coverage requirement, only one type was read technology, contained insertion calls that carried the
assigned according to the decision tree described in Fig. 1. Illumina (short reads) tag or a SR-based caller tag. For
Chaisson et al callsets (NA19240, HG00514 and HG0733),
Junctional homology detection the selection was performed on the vcf INFO field and
Junctional homology, as referred to and defined in [16], is the UNION variable. The UNION variable can take three
a DNA sequence that has two identical or nearly identi- potential values, Pacbio,Bionano or Illumina, that cor-
cal copies at the junctions of the two genomic segments responded to the sequencing technology allowing the
involved in the rearrangement. In the case of an inser- variant to be discovered. For the GiaB callset (HG002),
tion, a junctional homology is a sequence segment at insertions that could be discovered with short reads were
the left (resp. right) side of the insertion site which is identified by the Ill tag contained in the ExactMatchID
nearly identical to the end (resp. beginning) of the inserted located in the INFO field of the vcf file. Insertion calls
sequence. Small junctional homologies (<10 bp on each that were labelled Ill only with refining methodologies
side) were searched in a strict manner by scanning simul- and not any discovery methodologies were not taken
taneously the 10 bp sequence at the left (resp. right) side into account for the Short read technology part. The sec-
of the insertion site and the 10 bp end (resp. beginning) ond part, referred as Other technologies, contained all the
of the inserted sequence, counting the number of iden- remaining insertions. It should be noted that all insertion
tical nucleotides starting from the insertion site until a calls in the first part carried also at least one long read
mismatch is encountered. For larger homologies, both technology tag and were not discovered using only short
the 100% identity and strict adjacency to the insertion read technology.
site constraints were relaxed. We used the local align-
ments between the breakpoint junctions and the inserted Simulations
sequence that were previously obtained with BLAT. Only Twenty two sequencing datasets were simulated to char-
the alignments with at least 90% identity and occurring at acterize the impact of the different insertion features
a maximum of 10 bp before (resp. after) the insertion site on SR-based insertion variant calling. Each dataset was
and at a maximum of 10 bp from the end (resp. begin- obtained by altering the human chromosome 3 with 200
ning) of the inserted sequence were retained. In case of insertions. Sequencing reads were generated using ART
multiple candidates hits at one side of the junction, the with the following parameters : 2x150 bp reads, at 40 X
one located at the closest position from the extremities coverage, with insert size of 300 bp on average and 20 bp
was kept. If homologies (small or large) were found at standard deviation [39].
both sides of the junction, the homology size was obtained
by summing both homology sizes after removing poten- Baseline simulation
tial overlap on the inserted sequence. To compute the The simulation referred as the baseline was meant to
expected distribution of junctional homology sizes that represent the easiest type of insertions to detect, where
could be observed by chance, we generated 2,000 ran- inserted sequences contained very few repeats and are
dom insertions on the human chromosome 3 sequence. novel in the genome, the genomic context of inser-
Inserted sequences were generated by concatenating 250 tion was also simple and repeat-free, and breakpoint
nucleotides sampled uniformly on the A,C,G,T alphabet. junctions did not have any homology. To do so, we
The insertion sites were sampled uniformly along the simulated 250 bp novel sequence insertions located in
chromosome 3 sequence. Junctional homology sizes of exons without any homology at the breakpoint junc-
Delage et al. BMC Genomics (2020) 21:762 Page 15 of 17
tions. Novel sequences were extracted from randomly 3. Finally, the 889 insertions were simulated as described
chosen exonic regions of the Saccharomyces cerevisae in the vcf file.
genome.
Insertion calling and benchmarking
Scenario 1: varying the insertion size Simulated reads were aligned with bwa against the hg38
Insertion locations used in the baseline simulation were reference genome, and read duplicates were marked with
kept and the 200 inserted sequences were alterna- samblaster v.0.1.24 and converted into bam file with sam-
tively replaced by sequences extracted from Saccha- tools v1.6 [40, 41]. Bam index and reference dictionary
romyces cerevisae exons of 3 different sizes: 50, 500, and were obtained by picard tools v2.18.2. GRIDSS v2.8.0,
1000 bp. Manta v1.6.0, MindTheGap v2.2.1 and SvABA v1.1.0 were
all run using recommended, or otherwise default, param-
Scenario 2: varying the insertion type
eters [6, 10, 11, 20]. Only “PASS” insertions, that were
Insertion locations were identical to the baseline simula-
larger than 50 bp, were kept for the recall calculation. Two
tion, but the 250 bp inserted sequences were alternatively
types of recalls were computed depending on the pre-
replaced by dispersed duplications, tandem repeats, tan-
cision and information given for each call: insertion-site
dem duplications and mobile elements. Two types of tan-
only recall and sequence-resolved recall. The insertion-
dem repeats were simulated, with a pattern size of 6 bp
site only recall was assessed solely based on the inser-
or 25 bp, the pattern originating from the left breakpoint
tion site location prediction with a 10 bp margin around
junction. As mobile elements, 200 Alu mobile element
the expected location. As a more stringent evaluation,
sequences with a size ranging between 200 and 300 bp
the sequence-resolved recall took also into account the
were randomly extracted from the human genome based
inserted sequence. When it was reported, the inserted
on the RepeatMasker annotation. Tandem duplications
sequence had to share at least 90% of sequence identity
were generated by duplicating the 250 bp right breakpoint
to the simulated one and had to have a similar size of +/-
sequence. The inserted sequences of simulated dispersed
10%, to be considered as a true positive. In case of absence
duplications were extracted from exons of the chromo-
of alternative sequence in the vcf file but the provided
some 3.
annotation of the event allowed us to extract the insertion
Scenario 3 : varying the junctional homology size sequence from the reference genome (for instance for dis-
The 250 bp insertion sequences produced in the baseline persed duplication with the duplicated copy coordinates),
simulation were altered with junctional homology. To sim- it was evaluated similarly as for alternative sequences.
ulate junctional homologies, we replaced the X first bases Recall was computed as the ratio between the amount
of each insertion with the same size sequence originat- of true positive discoveries and the amount of simulated
ing from the right breakpoint sequence. We simulated five insertions. We compared the absolute amounts of false
junctional homology sizes (X value): 10, 20, 50, 100 and positive discoveries between tools and simulations, rather
150 bp. the precision or FDR metrics, as the latter are dependant
of the amount of true positive discoveries.
Scenario 4 : varying the genomic context of insertion
The 250 bp insertions from the baseline simulation were Long read simulation and benchmark
alternatively inserted in specific genomic contexts : either For each short read simulated dataset, a correspond-
inside different types of mobile elements, namely SINEs ing PacBio long read simulated dataset was produced,
and LINEs, in small (<300 bp) and large (>300 bp) simple using Simlord at 40 X coverage with probabilities of dele-
repeats or in other regions not annotated by Repeat- tion, insertion and substitution equal to 11%, 4% and
Masker (non repeated). A dataset with closely located 1% respectively [42]. Reads were aligned with Minimap2,
variants was simulated by adding insertions closed to the alignments were sorted with samtools and variants were
insertions simulated in the baseline scenario. The distance called with Sniffles [17, 43]. The evaluation of insertion
between insertions varied uniformly from 5 to 150 bp. site recalls followed the same process than for short read-
based variant callers. For the sequence-resolved recall,
Scenario 5 : Real insertions two sequence identity thresholds, 90 and 80%, were used
The 889 insertions located on the chromosome 3 from to validate the inserted sequences. We also used the eval-
the NA19240 callset were used to simulate three addi- uation tool from GiaB, SVBenchmark module from SVan-
tional datasets. Novel sequences were first simulated at alyzer tools suite, with parameters similarly set as our
the real chromosome 3 locations, then the real insertions benchmark method: -minsize set to 50 bp and -maxdist to
were simulated inside exonic regions of the chromosome 10 bp.
Delage et al. BMC Genomics (2020) 21:762 Page 16 of 17
Supplementary Information 5. Pirooznia M, Kramer M, Parla J, Goes FS, Potash JB, McCombie WR,
The online version contains supplementary material available at Zandi PP. Validation and assessment of variant calling pipelines for
https://doi.org/10.1186/s12864-020-07125-5. next-generation sequencing. Hum Genomics. 2014;8(1):14.
6. Wala JA, Bandopadhayay P, Greenwald N, Rourke RO, Sharpe T, Stewart
C, Schumacher S, Li Y, Weischenfeldt J, Yao X, Nusbaum C, Campbell P,
Additional file 1: Supplementary Figures and Tables. Getz G, Meyerson M, Zhang C-Z, Imielinski M, Beroukhim R. SvABA:
genome-wide detection of structural variants and indels by local
assembly. Genome Res. 2018;28(4):581-91. https://doi.org/10.1101/gr.
Abbreviations
221028.117.
SV: Structural variation; TE: Transposable element; ME: Mobile element;
SR-based: Based on short reads 7. Kosugi S, Momozawa Y, Liu X, Terao C, Kubo M, Kamatani Y.
Comprehensive evaluation of structural variation detection algorithms for
Acknowledgments whole genome sequencing. Genome Biol. 2019;20:117. https://doi.org/
We are thankful to the Genouest bioinformatics platform, computations have 10.1186/s13059-019-1720-5.
been made possible thanks to their computing resources. We are grateful to 8. Mahmoud M, Gobet N, Cruz-Dávalos DI, Mounier N, Dessimoz C,
Justin Zook for his helpful advices to filter the GiaB callset. Sedlazeck FJ. Structural variant calling: the long and the short of it.
Genome Biol. 2019;20:246. https://doi.org/10.1186/s13059-019-1828-7.
Authors’ contributions 9. Cameron DL, Stefano LD, Papenfuss AT. Comprehensive evaluation and
WD, JT and CL conceived the study. WD developed the annotation and characterisation of short read general-purpose structural variant calling
simulation scripts and carried out the analysis of the results. All author(s) software. Nat Commun. 2019;10:324. https://doi.org/10.1038/s41467-
contributed to the writing and read and approved the final manuscript. 019-11146-4.
10. Rizk G, Gouin A, Chikhi R, Lemaitre C. Mindthegap : integrated detection
Funding and assembly of short and long insertions. Bioinformatics. 2014;30(24):
Not applicable. 3451–7. https://doi.org/10.1093/bioinformatics/btu545.
11. Cameron DL, Schröder J, Penington JS, Do H, Molania R, Dobrovic A,
Availability of data and materials Speed TP, Papenfuss AT. Gridss: sensitive and specific genomic
The human reference genome version for this study was Hg38 (GRCh38). SV rearrangement detection using positional de bruijn graph assembly.
callsets analysed in this study are publicly available (dbVar accessions: nstd152 Genome Res. 2017;27(12):2050–60. https://doi.org/10.1101/gr.222109.
and nstd175 for Chaisson et al and Genome in a Bottle callsets respectively), 117.
and were downloaded from the following links: 12. Lappalainen I, Lopez J, Skipper L, Hefferon T, Spalding JD, Garner J, Chen
NA19240: http://ftp.ncbi.nlm.nih.gov/pub/dbVar/data/Homo_sapiens/by_ C, Maguire M, Corbett M, Zhou G, et al. Dbvar and dgva: public archives
study/genotype/nstd152/NA19240.BIP-unified.vcf.gz. for genomic structural variation. Nucleic Acids Res. 2013;41(D1):936–41.
HG00514: http://ftp.ncbi.nlm.nih.gov/pub/dbVar/data/Homo_sapiens/by_ 13. Abnizova I, te Boekhorst R, Orlov Y. Computational errors and biases of
study/genotype/nstd152/HG00514.BIP-unified.vcf.gz. short read next generation sequencing. J Proteomics Bioinform.
HG00733: http://ftp.ncbi.nlm.nih.gov/pub/dbVar/data/Homo_sapiens/by_ 2017;10(1):1–17.
study/genotype/nstd152/HG00733.BIP-unified.vcf.gz. 14. Conrad DF, Bird C, Blackburne B, Lindsay S, Mamanova L, Lee C, Turner
HG002: http://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/AshkenazimTrio/ DJ, Hurles ME. Mutation spectrum revealed by breakpoint sequencing of
analysis/NIST_SVs_Integration_v0.6/HG002_SVs_Tier1_v0.6.vcf.gz. human germline cnvs. Nat Genet. 2010;42(5):385.
Reference Genome,GRCh38/hg38: http://hgdownload.soe.ucsc.edu/ 15. Kidd JM, Graves T, Newman TL, Fulton R, Hayden HS, Malig M, Kallicki J,
goldenPath/hg38/bigZips/hg38.fa.gz Kaul R, Wilson RK, Eichler EE. A human genome structural variation
Custom scripts used in this study are freely available at https://github.com/ sequencing resource reveals insights into mutational mechanisms. Cell.
WesDe/DeepAn and at https://github.com/WesDe/InserSim. 2010;143(5):837–47.
16. Ottaviani D, LeCain M, Sheer D. The role of microhomology in genomic
Ethics approval and consent to participate structural variation. Trends Genet. 2014;30(3):85–94.
Not applicable. 17. Sedlazeck FJ, Rescheneder P, Smolka M, Fang H, Nattestad M,
von Haeseler A, Schatz MC. Accurate detection of complex structural
Consent for publication variations using single-molecule sequencing. Nat Methods. 2018;15(6):
Not applicable. 461–8.
18. Chaisson MJP, Sanders AD, ..., Marschall T, Korbel J, Eichler EE, Lee C.
Competing interests Multi-platform discovery of haplotype-resolved structural variation in
The authors declare that they have no competing interests. human genomes. Nat Commun. 2019;10:1784. https://doi.org/10.1038/
s41467-018-08148-z.
19. Zook JM, Hansen NF, ..., Chaisson MJ, Spies N, Sedlazeck FJ, Salit M, the
Author details
1 Univ Rennes, Inria, CNRS, IRISA, F-35000 Rennes, France. 2 Inserm U1209, CNRS Genome in a Bottle Consortium. A robust benchmark for detection of
germline large deletions and insertions. Nat Biotechnol. 2020. https://doi.
UMR 5309, Univ. Grenoble Alpes, Institute for Advanced Biosciences, Grenoble,
org/10.1038/s41587-020-0538-8.
France & Genetics, Genomics and Reproduction Service, Centre
20. Chen X, Schulz-Trieglaff O, Shaw R, Barnes B, Schlesinger F, Källberg M,
Hospitalo-Universitaire Grenoble-Alpes, Grenoble, France.
Cox AJ, Kruglyak S, Saunders CT. Manta: rapid detection of structural
variants and indels for germline and cancer sequencing applications.
Received: 7 July 2020 Accepted: 6 October 2020
Bioinforma (Oxford, England). 2016;32:1220–2. https://doi.org/10.1093/
bioinformatics/btv710.
21. Boycott KM, Vanstone MR, Bulman DE, MacKenzie AE. Rare-disease
References genetics in the era of next-generation sequencing: discovery to
1. Zook JM, Chapman B, Wang J, Mittelman D, Hofmann O, Hide W, Salit translation. Nat Rev Genet. 2013;14(10):681–91.
M. Integrating human sequence data sets provides a resource of 22. Wellenreuther M, Mérot C, Berdan E, Bernatchez L. Going beyond snps:
benchmark SNP and indel genotype calls. Nat Biotechnol. 2014;32(3):246. the role of structural genomic variants in adaptive evolution and species
2. Mullaney JM, Mills RE, Pittard WS, Devine SE. Small insertions and diversification. Mol Ecol. 2019;28(6):1203–9.
deletions (indels) in human genomes. Hum Mol Genet. 2010;19(R2):131–6. 23. Chander V, Gibbs RA, Sedlazeck FJ. Evaluation of computational
3. Baker M. Structural variation: the genome’s hidden architecture. Nat genotyping of structural variation for clinical diagnoses. GigaScience.
Methods. 2012;9(2):133–7. 2019;8(9):giz110. https://doi.org/10.1093/gigascience/giz110.
4. Feuk L, Carson AR, Scherer SW. Structural variation in the human 24. Krusche P, Trigg L, Boutros PC, Mason CE, Francisco M, Moore BL,
genome. Nat Rev Genet. 2006;7(2):85–97. Gonzalez-Porta M, Eberle MA, Tezak Z, Lababidi S, et al. Best practices for
Delage et al. BMC Genomics (2020) 21:762 Page 17 of 17
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in
published maps and institutional affiliations.
A.3 Publication 3
SVJedi : Genotyping structural variations with long reads
Lolita Lecompte, Pierre Peterlongo, Dominique Lavenier, Claire Lemaitre.
Bioinformatics 2020 36(17) :4568–4575
92
Bioinformatics, 36(17), 2020, 4568–4575
doi: 10.1093/bioinformatics/btaa527
Advance Access Publication Date: 21 May 2020
Original Paper
Sequence analysis
SVJedi: genotyping structural variations with long reads
Abstract
Motivation: Studies on structural variants (SVs) are expanding rapidly. As a result, and thanks to third generation
sequencing technologies, the number of discovered SVs is increasing, especially in the human genome. At the
same time, for several applications such as clinical diagnoses, it is important to genotype newly sequenced individu-
als on well-defined and characterized SVs. Whereas several SV genotypers have been developed for short read
data, there is a lack of such dedicated tool to assess whether known SVs are present or not in a new long read
sequenced sample, such as the one produced by Pacific Biosciences or Oxford Nanopore Technologies.
Results: We present a novel method to genotype known SVs from long read sequencing data. The method is based
on the generation of a set of representative allele sequences that represent the two alleles of each structural variant.
Long reads are aligned to these allele sequences. Alignments are then analyzed and filtered out to keep only inform-
ative ones, to quantify and estimate the presence of each SV allele and the allele frequencies. We provide an imple-
mentation of the method, SVJedi, to genotype SVs with long reads. The tool has been applied to both simulated and
real human datasets and achieves high genotyping accuracy. We show that SVJedi obtains better performances
than other existing long read genotyping tools and we also demonstrate that SV genotyping is considerably
improved with SVJedi compared to other approaches, namely SV discovery and short read SV genotyping
approaches.
Availability and implementation: https://github.com/llecompte/SVJedi.git
Contact: lolita.lecompte@inria.fr
Supplementary information: Supplementary data are available at Bioinformatics online.
C The Author(s) 2020. Published by Oxford University Press. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com
V 4568
SVJedi 4569
SV2 (Antaki et al., 2018), Nebula (https://www.biorxiv.org/content/ sequences. It outputs a variant file complemented with the individ-
10.1101/566620v1). Though short reads are often used to discover ual genotype information for each input variant in VCF format.
and genotype SVs, this is well known that their short size makes The method consists of four different steps that are illustrated in
them ill-adapted for predicting large SVs or SVs located in repeated Figure 1. The fundamentals of the method lie in its first step, which
regions. SVs are often located alongside repeated sequences such as generates representative allele sequences that represent the two
mobile elements (Kidd et al., 2010), resulting in mappability issues alleles of each SV. Long reads are then aligned on the whole set of
that make the genotyping problem harder when using short read representative allele sequences. An important step consists in select-
data. ing and counting only informative alignments to finally estimate the
Third generation sequencing technology, such as Pacific genotype for each input variant.
Biosciences (PacBio) and Oxford Nanopore Technologies (ONT),
can produce much longer reads compared to NGS technologies.
2.1 Representative allele sequence generation
Fig. 1. SVJedi steps for deletion genotyping. Steps for insertion genotyping are symmetrical and are not shown on the figure for clarity purposes. 1. Two corresponding repre-
sentative allele sequences are generated for each selected SV, one corresponds to the original sequence and the other to the sequence with the deletion. 2. Long read sequenced
data are aligned on these allele sequences using Minimap2. 3. Informative alignments are selected. 4. Genotypes are estimated
4570 L.Lecompte et al.
right sides of the breakpoint, they must satisfy the following condi- Otherwise it is composed of two sequences of size 2 Ladj each cen-
tion in Equation (1) for the alignment to be kept: tered on each breakpoint. We therefore apply the following formula
to compute the normalized read count for the reference allele, c0 , as
ðx > dover Þ & ðy > dover Þ (1) a function of the observed read count for the reference allele, c0, and
Concerning the filtering of spurious false-positive alignments, the deletion size, delsize:
Minimap2 alignments are first filtered based on the mapping quality 8
>
> 2 Ladj
(MAPQ) score. To focus uniquely on mapped reads, the MAPQ < c0 if delsize < 2 Ladj
c0 ¼ ð2 Ladj þ delsize Þ
score of the alignments must be >10. This is not sufficient to filter >
> 1
out alignments due to repetitive sequences since mapping is per- : c0 otherwise
formed on a small subset of the reference genome and these align- 2
ments may appear as uniquely mapped on this subset. (3)
ða1 < dend jj r1 < dend Þ & ða2 < dend jj r2 < dend Þ (2)
where err is the probability that a read maps to a given allele errone-
ously, assuming it is constant and independent between all observa-
The left member of Equation (2) imposes that the unaligned part tions. err was fixed to 5:105 , after empirical experiments on a
at the left of the alignment is small in at least one of the two aligned simulated dataset (see Supplementary Fig. S1).
sequences; the right member imposes the same condition at the right Finally, the genotype with the largest likelihood is assigned and
side of the alignment. all three likelihoods are also output (-log10 transformed) as add-
itional information in the VCF file.
2.4 Genotype estimation
For each variant, the genotype is estimated based on the ratio of 2.5 Implementation and availability
amounts of reads informatively aligned to each allele sequence. In We provide an implementation of this method named SVJedi, freely
the case of deletions and insertions, the allele sequences of a given available at https://github.com/llecompte/SVJedi, under the GNU
variant are of different sizes and contain a different number of Affero GPL license. SVJedi can also be installed from Bioconda.
breakpoints (for a deletion for instance, the reference allele contains SVJedi is written in Python 3, it requires as input a set of SVs (VCF
two breakpoints, whereas the alternative allele contains only one), format), a reference genome (fasta format) and a sequencing read
so even if both alleles are covered with the same read depth, there file (fastq or fasta format). Notably, the main steps are implemented
would be fewer reads that align on the shortest allele sequence and in a modular way, allowing the user to start or re-run the program
that overlap at least one breakpoint. To prevent a bias toward the from previous intermediate results. As an example, the first step is
longest allele, reported read counts for the longest allele are normal- not to be repeated if there are several long read datasets to be geno-
ized according to the allele sequence length ratio, assuming that read typed on the same SV set. Results shown here were obtained with re-
count is proportional to the sequence length. More precisely, in the lease version 1.1.0.
case of a deletion, the reference allele is the longest allele. Its allele
sequence size is the deletion size plus 2 Ladj (cumulative size of the
adjacent sequences) if the deletion is smaller than 2 Ladj . 3 Materials
3.1 Long read simulated dataset
SVJedi was assessed on simulated datasets on the human chromo-
some 1 (assembly GRCh37) based on real characterized deletions
for the human genome. From the dbVar database (Phan et al.,
2016), we selected 1000 existing deletions on chromosome 1
(defined as <DEL> SV type), which are separated by at least
10 000 bp. The sizes of the deletions vary from 50 bp to 10 kb (with
median and average sizes of 950 and 2044 bp, respectively). In this
experiment, deletions were distributed into the three different geno-
Fig. 2. Definition of the different distances used to select informative alignments be- types: 333 deletions are simulated with 0/0 genotype, 334 deletions
tween an allele sequence and a read. The aligned parts of the sequences are illus- with 0/1 genotype and the 333 remaining deletions with 1/1 geno-
trated by vertical bars. The allele sequence is composed of two adjacent sequences
type. Two different sequences were simulated containing each over-
of size Ladj on either side of the breakpoint, which is represented by a red vertical
lapping sets of deletions, representing the two haplotypes of the
thick bar. x and y are the distances of the breakpoint to, respectively, the start and
end coordinates of the alignment on the allele sequence. a1 and a2 (resp. r1 and r2) simulated individual. 1/1 genotype deletions were simulated on both
are the distances of the alignment left and right extremities to the, respectively, left haplotype sequences, whereas deletions of 0/1 genotype were simu-
and right extremities of the allele sequence (resp. read sequence). It follows here, lated each on one randomly chosen of the two haplotype sequences.
that a1 þ x ¼ Ladj and a2 þ y ¼ Ladj Then PacBio data were simulated on both haplotypes, using
SVJedi 4571
SimLoRD (Stöcker et al., 2016) (version v1.0.2) with varying Table 1. Contingency table of SVJedi results on PacBio simulated
sequencing error rates (6%, 10%, 16% and 20%), and at varying data (30x) of human chromosome 1 with 1000 deletions from
total sequencing depths (6, 10, 16, 20, 30, 40, 50 and dbVar
60). Most results presented in the main text are for 16% error rate
and 30 sequencing depth. Ten such datasets were simulated to as- SVJedi predictions
sess the reproducibility of results. 0/0 0/1 1/1 ./.
0/0 331 1 0 1
Truth
3.2 Real data 0/1 3 330 0 1
SVJedi was applied on a real human dataset, from the individual 1/1 0 18 313 2
HG002, son of the so-called Ashkenazi trio dataset. A PacBio Genotyping accuracy: 97.8 %
Table 2. Comparison of several tools and approaches for genotyping the 12 745 deletions and insertions of the GiaB call set in the HG002
individual
Deletions Insertions
Tool Genotyping accuracy Genotyping rate Genotyping accuracy Genotyping rate Time
Notes: Three approaches are compared: using long read genotyping tools (first three tools), using a short read genotyping tool (SVTyper), and using long read
discovery tools (last two tools). Except for the short read genotyping tool (SVtyper) that uses a 30 Illumina sequencing dataset, all other tools were run with a
30 PacBio long read dataset. Runtimes were measured on a 40-CPU computing node.
a
svviz2 is not parallelized.
Fig. 5. Stratified analysis of genotyping accuracy and rate for several genotyping
Fig. 4. Results of genotyping tools for the 12 745 deletions and insertions from the tools with respect to the genomic context of the SVs. Three categories of genomic
GiaB call set in the HG002 individual according to different SV size classes: 50– context are considered: tandem repeats >100 bp (TR), segmental duplications larger
100 bp, 100–250 bp, 250 bp to 1 kb, 1–10 kb and 10 kb. The two figures on top than 10 kb (SeqDup) and all other regions (Other)
represent the genotyping accuracies and the genotyping rates of SVJedi, Sniffles–Ivcf
and svviz2 on a 30 PacBio dataset, and of SVtyper for deletions only on a 30
Illumina dataset. The bottom figure represents the SV count of each SV size class We then compared the genotyping performances with respect to
the genomic context of the SVs (Fig. 5). As shown previously, SVs
falling in a Tandem Repeat >100 bp are harder to genotype with
is observed for small SVs (<100 bp), but, apart from this size class, SVJedi, with a 9% decrease of accuracy for these SVs, compared to
its accuracy is quite robust with respect to the size of SVs. On the
those outside TRs. A larger decrease of genotyping accuracy (14%)
opposite, svviz2 obtained its best genotyping accuracy for the small-
is observed with Sniffles–Ivcf in these regions. Although less fre-
est SVs (<100 bp) and it rapidly drops for SVs larger than 250 bp,
quent (here, only 48 concerned SVs), large segmental duplications,
falling below 30% for SV sizes between 1 and 10 kb. When compar-
typically larger than 10 kb, are also likely to affect long read map-
ing between SV types, svviz2 genotyping accuracy is significantly
ping accuracy and thus genotyping accuracy. SVJedi accuracy
lower for insertions than deletions, with, in particular, <10% of the
seemed not to be affected by these duplications, contrary to Sniffles–
insertions larger than 1 kb that are correctly genotyped (see
Ivcf (Fig. 5).
Supplementary Fig. S4). This inability for genotyping large SVs can
probably be explained by the way svviz2 identifies informative reads
for a given SV: only the reads mapped initially to the reference gen- 4.3.2 Comparison with a short read based genotyping approach
ome are selected before re-aligning them against both the reference For this same individual (HG002), some short read datasets are also
and alternative alleles. Consequently, most reads coming from large available. We, therefore, can compare SV genotyping performances
insertion alternative alleles could probably not be used for estimat- between two approaches and data types, namely long versus short
ing these genotypes. To a lesser extend, Sniffles–Ivcf genotyping ac- reads. SVJedi predictions were compared to an SV genotyping tool
curacy is also lower for large insertions than large deletions (69.6% for short reads, SVtyper, known as a reference tool in the state of
for insertions versus 85.7% for deletions, 1 kb), whereas SVJedi the art (Chiang et al., 2015; Chander et al., 2019). Since SVtyper
genotyping accuracy is unaffected by SV type for all size classes. does not support insertion variants, we focus here only on deletions,
4574 L.Lecompte et al.
and the 5464 deletions from the GiaB call set were genotyped with variants once they have been discovered, at least with long read data
SVtyper using a 2 250 bp 30x Illumina read dataset of HG002. as was shown here. Indeed, without a priori SV discovery is a much
Table 2 shows that more than half of the deletions are genotyped harder task than genotyping. Because the alternative allele is not
differently by SVtyper than in the high confidence GiaB call set, known in discovery, discovery methods rely on fewer or noisier sig-
resulting in a genotyping accuracy of only 46.5%, while this per- nals to identify the SVs than genotyping methods. Consequently,
centage rises to 91.7% for SVJedi with long reads. Remarkably, both approaches would likely benefit from different optimal param-
many of the discrepancies of SVtyper with GiaB are totally contra- eter settings, with for instance discovery methods requiring a more
dictory with 0/0 genotypes instead of 1/1 ones (see Supplementary stringent set of parameters to limit the false discovery rate.
Table S5). We can clearly see, in Figure 5, that short read based gen- However, in this paper, we have no intention to oppose both
otyping is much more impacted by the presence of TRs at the break- approaches but rather argue that they are complementary and are
point. As expected, mapping reads in these regions is much more intended to different purposes: when the aim is strictly to genotype
Li,H. (2011) A statistical framework for SNP calling, mutation discovery, as- Phan,L. et al. (2016) dbVar structural variant cluster set for data analysis and
sociation mapping and population genetical parameter estimation from variant comparison. F1000Research, 5, 673.
sequencing data. Bioinformatics, 27, 2987–2993. Sedlazeck,F.J. et al. (2018) Accurate detection of complex structural variations
Li,H. (2018) Minimap2: pairwise alignment for nucleotide sequences. using single-molecule sequencing. Nat. Methods, 15, 461–468.
Bioinformatics, 34, 3094–3100. Spies,N. et al. (2015) svviz: a read viewer for validating structural variants.
Lupski,J.R. (2015) Structural variation mutagenesis of the human genome: im- Bioinformatics, 31, 3994–3996.
pact on disease and evolution. Environ. Mol. Mutagen., 56, 419–436. Stancu,M.C. et al. (2017) Mapping and phasing of structural variation in pa-
Merker,J.D. et al. (2018) Long-read genome sequencing identifies causal struc- tient genomes using nanopore sequencing. Nat. Commun., 8, 1326.
tural variation in a Mendelian disease. Genet. Med., 20, 159–163. Stöcker,B.K. et al. (2016) SimLoRD: simulation of long read data.
Nielsen,R. et al. (2011) Genotype and SNP calling from next-generation Bioinformatics, 32, 2704–2706.
sequencing data. Nat. Rev. Genet., 12, 443–451. Zook,J.M. et al. (2016) Extensive sequencing of seven human genomes to
Norris,A.L. et al. (2016) Nanopore sequencing detects structural variants in characterize benchmark reference materials. Sci. Data, 3, 160025.
Etat civil
Claire Lemaitre, née le 7 novembre 1982 à Grenoble (38), vie maritale, 2 enfants, nationalité
française.
Coordonnées professionnelles
Equipe Genscale
Inria Rennes Bretagne Atlantique
Campus de Beaulieu
35042 Rennes cedex
tel : 02 99 84 71 16
claire.lemaitre@inria.fr
Page personnelle : http://people.rennes.inria.fr/Claire.Lemaitre/
Situation actuelle
Études et diplômes
1
Institut National des Sciences Appliquées
1
Expérience professionnelle
Activités d’enseignement
Encadrement
Développement de logiciels
Suite du CV
2
Activités d’enseignement
Actuellement
2 demies UEs d’algorithmique des séquences (pattern matching, structures d’indexation, aligne-
ment de séquences, assemblage de séquences), pour des publics différents :
Enseignements passés
• Algorithmique des séquences, niveau M1, parcours informatique, ENS de Rennes (12 heures
par an), 2017.
• TPs de biostatistiques sous R, niveau L3, licence de biologie de l’université de Rennes 1 (12
heures par an), 2013 à 2017.
3
Encadrements
• Encadrements en cours :
• Encadrements passés :
4
– Anaı̈s Gouin, ingénieure, 3 ans (2012-2015) : gestion et analyse des données de séquençage
d’insectes ravageurs de culture.
Co-encadrement avec Fabrice Legeai, taux réel 70 %.
Production : 2 publications (P18: Heredity, P12: Scientific Reports).
Devenir : ingénieure dans le privé, Aix en Provence.
– Liviu Ciortuz, post-doctorant, 1 an (2012-2013) : développement de nouveaux algo-
rithmes pour la détection de variants complexes (variants de structure) dans des données
NGS non assemblées.
Co-encadrement avec Pierre Peterlongo, taux réel 70 %.
Production: 1 publication (C1: Conférence AlCoB), 1 logiciel (L11: TakeABreak).
Devenir: Maı̂tre de Conférence à l’université de Iasi en Roumanie.
– Erwann Scaon, doctorant, 1 an (2012-2013) : Algorithmes pour l’assemblage de novo
de génomes fortement redondants.
Arrêt de la thèse à la fin de la 1ère année pour des raisons personnelles.
Co-encadrement avec Dominique Lavenier, taux réel 90 %.
Devenir: CDD ingénieur bioinformaticien en France.
– Thomas Derrien, post-doctorant, 1 an (2011-2012) : détection de variants structuraux
dans les génomes re-séquencés du puceron du pois.
Encadrement à 100 %.
Devenir: CR CNRS titulaire à l’Institut de Génétique de Rennes (IGDR), depuis 2012.
• Encadrements de stagiaires :
– Arthur Le Bars, stage de 6 mois niveau M2 (2019) : assemblage de génomes avec des
données linked-reads.
– Wesley Delage, stage de 6 mois niveau M2 (2017) : reconstruction de génomes par
assemblage guidé dans des données métagénomiques.
– Charlotte Mouden, stage de 6 mois niveau M2 (2017) : Adaptation et application de
l’outil discoSnp à des données de RAD-seq (logiciel L3: DiscoSnpRAD).
– Mael Kerbiriou, stage de 6 mois niveau M2 (2017) : Amélioration du logiciel discoSnp.
(taux 30%)
– Mathieu Perrotin, stage de 6 mois niveau M2 (2016) : Développement d’une méthode de
requêtage dans des fichiers de séquençages métagénomiques. (taux 30 %)
– Pierre Marijon, stage de 4 mois niveau M1 (2015) : amélioration d’un logiciel existant
(L1: MindTheGap).
– Gaëtan Benoit, stage de 6 mois niveau M2 (2014) : développement d’une méthode de
compression des fichiers de séquençage génomique.
Production: 1 publication (P16: BMC Bioinformatics, 1 logiciel (L9: Leon), .
– Julien Erabit, stage de 6 mois niveau M2 (2013) : développement d’un environnement de
benchmark d’assemblage de génomes polyploı̈des.
– Gaëtan Benoit, stage de 3 mois niveau M1 (2013) : développement d’une méthode de
correction des reads.
Production: un logiciel (L10: Bloocoo).
– Bastien Hervé, stage de 5 mois niveau L2 (2013) : évaluation d’une méthode de barcoding
à partir de données NGS plein-génome non assemblées.
– Elise Larsonneur, stage de 6 mois niveau M2 (2012) : intégration de données “omics”,
application sur le puceron du pois.
– Quentin Oliveau, stage de 3 mois niveau M1 (2012) : détection d’insertions virales et de
points de jonction dans des données de séquences à haut débit.
5
– Olivier Rué, stage de 6 mois niveau M2 (2011) : détection de variants (SNP et SV) dans
plusieurs génotypes du puceron du pois re-séquencés par NGS.
6
Responsabilités, tâches collectives
• Contrats de recherche :
– Responsable de tâches et de partenaire de projets ANRs :
∗ Divalps (2021-2025), appel générique ”Terre vivante”, crédits alloués de 590 Ke,
porté par Laurence Desprès (CNRS LECA, Grenoble). Financement et encadrement
d’une thèse à venir (2021-2024).
∗ Supergene (2020-2024), appel générique ”Terre vivante”, crédits alloués de 650 Ke,
porté par Mathieu Joron (CNRS CEFE, Montpellier). Financement et encadrement
d’un post-doctorant (Pierre Morisse, 2 ans).
∗ SpecRep (2015-2019), appel générique ”Terre vivante”, crédits alloués de 490 Ke,
porté par Marianne Elias (MNHN, Paris). Financement et encadrement d’un post-
doctorant (Jérémy Gautier, 2 ans).
∗ Ada-Spodo (2012-2015), appel générique ”Biodiversité, évolution, écologie, agronomie”;
crédits alloués de 430 Ke, porté par Emmanuelle d’Alençon (INRA, Montpellier). Fi-
nancement et encadrement d’une ingénieure (Anaı̈s Gouin, 18 mois).
– Participation à d’autres projets ANRs (sans responsabilité de partenariat) :
∗ Hydrogen (2014-2019), porté par Dominique Lavenier, puis Pierre Peterlongo (Inria,
Rennes). Financement et encadrement d’une thèse (Gaëtan Benoit, 3 ans).
∗ Colib’read (2013-2016), porté par Pierre Peterlongo (Inria, Rennes).
∗ SpeciAphid (2011-2014), porté par Jean-Christophe Simon (INRA, Rennes). Finance-
ment et co-encadrement d’une ingénieure (Anaı̈s Gouin, 18 mois).
– Responsable du projet “Mirage” de la Région Bretagne, Stratégie Attractivité Durable,
durée : 1 an (2012-2013). Financement et encadrement d’un post-doctorant (Liviu
Ciortuz, 1 an).
– Responsable du projet PEPS Bio-Math-Info, “barcoding de nouvelle génération”, durée :
2 ans (2012-2014)
• Relectures :
– Relecture d’articles pour des journaux (Bioinformatics, Nature Reviews Genetics, Nu-
cleic Acids Research, BMC Evolutionary Biology, PLoS One, IEEE/ACM Transactions on
Computational Biology and Bioinformatics), et des conférences internationales de bioin-
formatique (RECOMB, WABI, ISMB).
– Membre du Comité de Programme des conférences nationales JOBIM (2012, 2013, 2017,
2019, 2021), SeqBio (2011, 2013, 2014), seqBIM (2019, 2020).
– Relecture d’1 projet ANR, programme blanc SVSE6 en 2012.
• Organisation de conférences :
– Co-présidente du Comité de Programme de la conférence JOBIM 2022. Cette manifesta-
tion est le rendez-vous annuel de la communauté bioinformatique francophone (∼500 par-
ticipants, 6 conférenciers invités, ∼ 250 soumissions (dont 50 proceedings, 200 posters)).
– Organisation (co-responsable des comités d’organisation et de programme) des journées
annuelles du Groupe de Travail seqBIM en 2019 et 2020 (2 jours, 60-100 participants, 2
orateurs invités et 15 exposés sélectionnés, chaque année).
– Membre du comité d’organisation de l’école Jeunes Chercheurs de Bioinformatique Moléculaire
((JC)2 BIM), 2018 et 2020-2021 (5 jours, 30 élèves par session).
– Organisation de journées de formation à la librairie GATB, 2018 et 2019 (1 jour, 10-15
participants par session).
7
– Co-présidente du Comité d’Organisation de JOBIM 2012. Cette manifestation est le
rendez-vous annuel de la communauté bioinformatique francophone (∼400 participants,
budget > 100 000 euros).
– Membre du comité d’organisation du workshop Seq BI 2011 à Rennes.
– Organisation de formations continues de bioinformatique (2009-2010, au sein de la plate-
forme de bioinformatique de Bordeaux).
– Membre de 5 jurys de thèses en tant que examinatrice (Joseph Lucas, 2016; Arnaud
Meng, 2017; Yoann Seeleuthner, 2018, Lyam Baudry, 2019, Guillaume Gautreau, 2020).
– Membre de comités de recrutement d’un poste de Maı̂tre de Conférence (Université de
Lyon, 2014) et de 2 postes d’ingénieurs de recherche (INRA et Université de Bordeaux,
2015), d’1 poste d’ATER (Université de Rennes, 2017).
– Membre de la commission CORDIS de recrutement de doctorants ou post-doctorants du
centre Inria Rennes Bretagne Atlantique en 2013, 2014, 2017.
– Membre de 7 comités de thèses entre 2015 et 2020.
– Membre de jurys de soutenances de stage (niveau master) : école d’ingénieur ESTBB en
2010 à Bordeaux, master de bioinformatique de l’université Rennes 1 en 2012 et 2017.
• Animation de communautés :
– Compilation du rapport annuel d’activité Inria de l’équipe Symbiose, puis Genscale (depuis
2011).
8
Publications & communications
Dans mon domaine de recherche, la bioinformatique (et en biologie également), l’ordre des auteurs
a de l’importance : les premières et dernières places montrent une contribution plus importante. En
particulier, le ou les auteurs en dernières positions sont les personnes qui ont supervisé le travail.
En bioinformatique, les publications se font majoritairement par des soumissions dans des journaux
plutôt que des participations à des conférences, notamment pour avoir plus d’impacts auprès des biolo-
gistes utilisateurs des méthodes. J’ai surligné ci-dessous en gris souligné les revues majeures de bioinformatique,
et en vert, les revues appartenant à des domaines applicatifs (génomique, écologie).
9
divergence of pea aphid host races.
Molecular Ecology, 2018, 27(16):3287-3300.
P12 A. Gouin, ... , C. Lemaitre, F. Legeai, E. d’Alençon, P. Fournier (Plus de 20 auteurs, co-
dernier auteur).
Two genomes of highly polyphagous lepidopteran pests (Spodoptera frugiperda, Noctuidae)
with different host-plant ranges.
Scientific Reports, 2017, 7:11816.
10
P21 N. Maillet, C. Lemaitre, R. Chikhi, D. Lavenier, Pierre Peterlongo.
Compareads: comparing huge metagenomic experiments.
BMC Bioinformatics, 2012, 13 (Suppl 19):S10.
P26 C. Lemaitre, M. Braga Dias Vieira, C. Gautier, M.-F. Sagot, E. Tannier, G.A.B. Marais.
Footprints of inversions at present and past pseudoautosomal boundaries in human sex chro-
mosomes.
Genome Biology and Evolution, 2009, 1(1):56–66.
Chapitres de livre
? C. Guyomar, C. Lemaitre.
Métagénomique et métatranscriptomique.
Dans Du texte aux graphes : méthodes et structures discrètes pour la bioinformatique, ISTE
Science Publishing, à paraı̂tre 2021.
11
Communications à des conférences internationales avec comité de lecture (sans
actes)
C10 C. Lemaitre, M. Braga Dias Vieira, C. Gautier, M.-F. Sagot, E. Tannier, G.A.B. Marais.
Footprints of inversions at present and past pseudoautosomal boundaries in human sex chro-
mosomes.
Integrative Post Genomics (IPG) 2008, Lyon, 19-20 novembre 2008.
12
C12 W. Delage, J. Thevenon, C. Lemaitre.
Towards a better understanding of the low discovery rate of short-read based insertion variant
callers.
JOBIM 2020, Montpellier, juillet 2020 (8 pages).
Séminaires invités
? Comparaison (massive) de (nombreux) metagénomes. Passons par les kmers pour passer à
l’échelle.
Journée scientifique sur ”le Microbiome” organisée par Biogenouest, Rennes, décembre 2016.
? Détection et analyse des points de cassure de réarrangements dans les génomes de mammifères.
Conférence invitée au Centre de BioInformatique de Bordeaux, Bordeaux, avril 2008.
13
Manuscrit de thèse
? C. Lemaitre
Réarrangements chromosomiques dans les génomes de mammifères : caractérisation des points
de cassure.
Université Claude Bernard Lyon 1, 2008.
https://tel.archives-ouvertes.fr/tel-00364265/
Pédagogie
? E. Billoir, C. Lemaitre, S. Charles.
La modélisation des réseaux de gènes : une situation-problème pour l’apprentissage de méthodes
mathématiques avancées.
Poster au 25e Congrès de l’Association Internationale de Pédagogie Universitaire (AIPU), Mont-
pellier, Mai 2008.
Article dans les actes de la conférence.
Médiation scientifique
? S. Alizon, F. Cazals, S. Guindon, C. Lemaitre, T. Mary-Huard, A. Niarakis, M. Salson, C.
Scornavacca, H. Touzet.
SARS-CoV-2 Through the Lens of ComputationalBiology: How bioinformatics is playing a key
role in the study of the virus and its origins.
Mars 2021, déposé sur HAL (https://hal-cnrs.archives-ouvertes.fr/hal-03170023).
Rapport de médiation scientifique produit sous l’égide du GDR BIM, à destination de non
scientifiques ou scientifiques d’autres disciplines que la bioinformatique.
? C. Lemaitre
Le génome humain, la bioinformatique et le métier de bioinformaticien(ne).
Dans le cadre du dispositif “La découverte de la recherche”.
2 Interventions pour des élèves de 1ère du Lycée Bertrand d’Argentré à Vitré, mai 2013.
14
Logiciels
En bioinformatique, les logiciels sont un élément très important pour la visibilité de la recherche
auprès des biologistes.
Tous ces logiciels sont distribués sous licence GNU GPL ou Affero GPL en Open Source, la
majorité d’entre eux sont déposés à l’Agence de Protection des Programmes (APP).
• Logiciels majeurs (C++, > 5000 lignes de code), maintenus sur la durée et qui évoluent
encore :
15
L9 Leon (2014) : logiciel de compression des fichiers de séquençage d’ADN, codé en C++
Evolution : intégré en 2018 dans la librairie GATB-core, en cours de valorisation indus-
trielle par l’entreprise Enancio (actuellement Illumina) (créée par un des auteurs du logiciel
: Guillaume Rizk).
https://github.com/GATB/leon
L10 Bloocoo (2013) : logiciel de correction des lectures de séquençage à très faible emprunte
mémoire, codé en C++.
https://github.com/GATB/bloocoo
L11 TakeABreak (2014) : logiciel de détection d’inversions sans génome de référence, codé
en C++
Preuve de concept, usage essentiellement interne.
https://github.com/GATB/TakeABreak
L12 Compareads (2012) : logiciel de comparaison de séquences, codé en C++.
http://alcovna.genouest.org/compareads
L13 Cassis (2008) : logiciel de génomique comparative (codé en perl et R).
http://pbil.univ-lyon1.fr/software/Cassis
16