Une histoire sur la façon dont les multi-cartes sont gérées par bowtie2.
Dans le comportement par défaut de bowtie2, même si la séquence est frappée plusieurs fois, un seul d'entre eux est signalé. Ensuite, la façon de choisir l'un d'entre eux est généralement choisie au hasard parmi celles avec le score d'alignement le plus élevé. Le document le dit également [^ 1]. … Mais dans certains cas, il semble que ce ne sera pas aléatoire et se solidifiera [^ 2], alors j'aimerais clarifier les conditions.
La plupart des conclusions sont écrites ici [^ 2], mais ce problème est dû à la génération aléatoire de bowtie2. La section option --non-deterministic du manuel d'instructions de Bowtie2 [^ 3] contient l'explication suivante.
Normally, Bowtie 2 re-initializes its pseudo-random generator for each read. It seeds the generator with a number derived from (a) the read name, (b) the nucleotide sequence, (c) the quality sequence, (d) the value of the
--seedoption. This means that if two reads are identical (same name, same nucleotides, same qualities) Bowtie 2 will find and report the same alignment(s) for both, even if there was ambiguity. When--non-deterministicis specified, Bowtie 2 re-initializes its pseudo-random generator for each read using the current time. This means that Bowtie 2 will not necessarily report the same alignment for two identical reads. This is counter-intuitive for some users, but might be more appropriate in situations where the input consists of many identical reads.
En résumé, bowtie2 initialise un générateur de nombres aléatoires pour chaque lecture, et sa valeur de départ semble dépendre des quatre suivantes.
--Nom du chef --Séquence de base
--seedEn bref, les prospects qui ont les mêmes informations, y compris le __name, sont conçus pour être mappés au même emplacement __.
Dans ce qui suit, confirmons que le résultat du mappage sera différent si le read name, quality array et --seed option sont différents.
Préparez le FASTQ suivant basé sur la séquence de base multi-mappée. (Les arrangements, etc. sont créés sur la base de l'exemple de vérification ici [^ 2])
--FASTQ avec le même nom de lecture, la même séquence de base et la même séquence de qualité --FASTQ avec la même séquence de base mais une séquence de nom / qualité différente {l'un ou les deux} pour chaque lecture
Comme il est difficile de se préparer à chaque fois, le générateur suivant sera utilisé à la place.
gen_testseq.py
#!/usr/bin/env python
import argparse
import random
parser = argparse.ArgumentParser()
parser.add_argument("-r", "--reads", type=int, default=10000)
parser.add_argument("-n", "--name", action="store_true")
parser.add_argument("-q", "--quality", action="store_true")
args = parser.parse_args()
for i in range(args.reads):
print "@test_sequence{}".format(random.random() if args.name else '')
print "CAAAGTGTTTACTTTTGGATTTTTGCCAGTCTAACAGGTGAAGCCCTGGAGATTCTTATTAGTGATTTGGGCTGGGGCCTGGCCATGTGTATTTTTTTAAA"
print '+'
print ''.join([chr(random.randrange(33, 127)) for _ in range(101)]) if args.quality else 'I' * 101
Lorsque le réseau sur la deuxième ligne d'impression est mappé sur hg19, 9 emplacements sont multi-hit, et les 3 emplacements suivants ont le score d'alignement le plus élevé et sont candidats pour la carte.
| chrom | position | strand |
|---|---|---|
| chr1 | 11635 | + |
| chr2 | 114359380 | - |
| chr15 | 102519535 | - |
bowtie2 a utilisé la version 2.2.4.
Ainsi, seule la position est retirée du résultat de la carte et le nombre est compté.
$ bowtie2 --no-hd -x hg19 -U <(./gen_testseq.py) | cut -f 3-4 | sort | uniq -c
résultat
10000 reads; of these:
10000 (100.00%) were unpaired; of these:
0 (0.00%) aligned 0 times
0 (0.00%) aligned exactly 1 time
10000 (100.00%) aligned >1 times
100.00% overall alignment rate
10000 chr15 102519435
Tous ont été cartographiés en un seul endroit.
--seed$ bowtie2 --no-hd --seed 17320508 -x hg19 -U <(./gen_testseq.py) | cut -f 3-4 | sort | uniq -c
10000 reads; of these:
10000 (100.00%) were unpaired; of these:
0 (0.00%) aligned 0 times
0 (0.00%) aligned exactly 1 time
10000 (100.00%) aligned >1 times
100.00% overall alignment rate
10000 chr1 11635
L'emplacement change en fonction de la valeur, mais un seul emplacement est mappé.
--non-deterministic$ bowtie2 --no-hd --non-deterministic -x hg19 -U <(./gen_testseq.py) | cut -f 3-4 | sort | uniq -c
10000 reads; of these:
10000 (100.00%) were unpaired; of these:
0 (0.00%) aligned 0 times
0 (0.00%) aligned exactly 1 time
10000 (100.00%) aligned >1 times
100.00% overall alignment rate
3302 chr1 11635
3379 chr15 102519435
3319 chr2 114359280
Cette fois, il a été mappé au hasard. Bien sûr, le résultat sera différent à chaque fois que vous l'exécuterez.
$ bowtie2 --no-hd -x hg19 -U <(./gen_testseq.py -n) | cut -f 3-4 | sort | uniq -c
10000 reads; of these:
10000 (100.00%) were unpaired; of these:
0 (0.00%) aligned 0 times
0 (0.00%) aligned exactly 1 time
10000 (100.00%) aligned >1 times
100.00% overall alignment rate
3379 chr1 11635
3294 chr15 102519435
3327 chr2 114359280
Cette fois, il a été mappé au hasard. ici,
$ ./gen_testseq.py -n > test.fastq
$ bowtie2 --no-hd -x hg19 -U test.fastq | cut -f 3-4 | sort | uniq -c
Si vous le déposez dans un fichier comme décrit ci-dessus, vous obtiendrez le même résultat à chaque fois que vous l'exécuterez. Bien entendu, cela ne s'applique pas si la valeur de l'option --seed est différente ou si l'option --non-deterministic est spécifiée.
Mappé aléatoirement même si le tableau de qualité est différent
$ bowtie2 --no-hd -x hg19 -U <(./gen_testseq.py -q) | cut -f 3-4 | sort | uniq -c
10000 reads; of these:
10000 (100.00%) were unpaired; of these:
0 (0.00%) aligned 0 times
0 (0.00%) aligned exactly 1 time
10000 (100.00%) aligned >1 times
100.00% overall alignment rate
3351 chr1 11635
3312 chr15 102519435
3337 chr2 114359280
--non-deterministic.
--Cependant, la reproductibilité est perdue.
--Utilisez l'option --seed si vous avez besoin de résultats avec différentes destinations multimap et que vous voulez assurer la reproductibilité.Il semble qu'il y ait eu un cas où ça ne variait pas bien [^ 2] mais était-ce un bug?
[^ 1]: Voir la section des options -k [^ 3]
Recommended Posts