Une note sur le comportement de bowtie2 lors de plusieurs coups

Une histoire sur la façon dont les multi-cartes sont gérées par bowtie2.

Contexte

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 --seed option. 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-deterministic is 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

En 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éparation

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.

Vérification

Lorsque toutes les informations de prospect sont les mêmes FASTQ

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.

Essayez de spécifier l'option --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é.

Essayez de spécifier l'option --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.

Si le nom ou la séquence de qualité est différent pour chaque lecture

$ 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

Conclusion

référence

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

Une note sur le comportement de bowtie2 lors de plusieurs coups
À propos du comportement de la file d'attente pendant le traitement parallèle
À propos du comportement de yield_per de SqlAlchemy
A propos du comportement de enable_backprop de Chainer v2
Un mémo expliquant la spécification de l'axe de l'axe
À propos du comportement de copy, deepcopy et numpy.copy
Écrire une note sur la version python de python virtualenv
Un mémorandum sur les avertissements dans les résultats de sortie de pylint
Un mémo pour comprendre visuellement l'axe des pandas.
Une histoire sur le changement du nom principal de BlueZ
Un mémorandum sur la mise en œuvre des recommandations en Python
Remarque sur le comportement par défaut de collate_fn dans PyTorch
J'ai essayé un peu le comportement de la fonction zip
À propos des composants de Luigi
À propos des fonctionnalités de Python
[python] Une note que j'ai commencé à comprendre le comportement de matplotlib.pyplot
À propos de la valeur de retour de pthread_mutex_init ()
À propos de la valeur de retour de l'histogramme.
À propos du type de base de Go
À propos de la limite supérieure de threads-max
À propos de la taille des points dans matplotlib
À propos de la liste de base des bases de Python
L'histoire de l'exportation d'un programme
Je voulais faire attention au comportement des arguments par défaut de Python
Une note sur les fonctions de la bibliothèque Linux standard qui gère le temps
Vérifiez le comportement du destroyer en Python
Classe qui atteint l'API de DMM
Une compréhension approximative de python-fire et un mémo
À propos de l'environnement virtuel de Python version 3.7
Mesurer la force de l'association dans un tableau croisé
Mémorandum sur le QueryDict de Django
A propos des arguments de la fonction setup de PyCaret
[python] [meta] Le type de python est-il un type?
À propos de l'équation normale de la régression linéaire
Obtenez le nom de fichier du répertoire (glob)
L'histoire du traitement A du blackjack (python)
Note de problèmes sur la coexistence du système Python 2/3
Mémo sur Sphinx Partie 1 (Création d'un projet)
Notez l'achèvement d'une commande chronophage
Un mémo que j'ai essayé le tutoriel Pyramid