À ma grande surprise, je n'ai pas retrouvé d'issue parlant de ces barcodes. Existe-t-il des scripts standard pour traiter ces données ? On en ferait un server-pre-process ? Des jeux publics sur lesquels on pourrait tester cela, ou bien un jeu de @mfigeac ?
Idéalement, une fois les reads associés entre eux par barcode moléculaire et la séquence consensus obtenue, l'information du nombre de reads portant le même barcode moléculaire doit être accessible pendant l'analyse. Si on peut aussi garder le nombre d'erreurs par position ca permet d'avoir un score de confiance dans chaque base.
Attention aux outils pour faire le démultiplexage car en fonction du protocole utilisé, la méthode pourra être différente.
Pour retrouver l'information, cela va être difficile (déjà stocker quelque part dans l'header Fasta ce nombre...). Mais on peut déjà faire #2143 pour pouvoir, certes de manière pas très commode, rechercher l'info si on en a besoin.
Oui, ce ne serait qu'une option de pre-process, et il pourrait y avoir plusieurs options selon la méthode.
Mathieu Giraudchanged title from Un pre-process pour traiter les identifiants moléculaires uniques to UMI / Un pre-process pour traiter les identifiants moléculaires uniques
changed title from Un pre-process pour traiter les identifiants moléculaires uniques to UMI / Un pre-process pour traiter les identifiants moléculaires uniques
J'ai rapatrié les fichiers (dossier umi sur mon espace vdb).
En local je réduis la taille des fichiers via le paramètre -U de vidjil. Je ne sais pas si c'est une bonne chose car je ne sais pas a quoi m'attendre sur la comparaison R1/R2 avec cette fonction.
A priori ces données sont faites avec 2x3nt (la combinatoire n'est pas très grande; 4096 possibilités)
J'ai installé calib (via conda) que je vais tester.
J'ai modifié le dernier script pour permettre de faire un export des fichiers fastq nettoyés des duplication de reads.
Avant d'exporter un lot de reads, je vérifie qu'ils sont identiques. Ils le sont, mais il y a parfois des mismatch sur un read (cf exemple).
Je ne sais pas comment faire simplement la comparaison pour choisir le meilleur read. Dans le fichier d'exemple, il y a 4 reads, dont un présente un mismatch. Dans d'autre cas, nous en avons seulement 2. Comment tranché dans ce cas.
Il faudrait développer une heuristique pour ça. Vous savez si on a un snippet qui pourrait faire ça quelques part ?
Au passage. Je suis passé dans le cas présent par un filtre VDJ en amont. En routine, il faudrait probablement passer directement sur les fichiers sources.
Je ne sais pas exactement comment marche calib_cons, mais ce n'est pas très probant. Prenons le cas de l'exemple 1:
nombre reads VDJ
805
calib + scipt
711
calib + calib_cons
105
Si je regarde le résultat de calib_cons, la très grande majorité des séquences sont composées pour moitié de N. Il me semble que ça agrège un peu trop les séquences et que l'on en tirera rien de bon.
En même temps, avec 2x3nt d' UMI, c'est peut-être court. Je n'imagine pas sur les données complètes de séquençage. J'ai lancé calib dessus. je voulais attendre pour y inclure ses résultats, mais ça fait plus d'une heure que ça tourne sans finir...
L'idée n'est pas de redévelopper de nouvelles choses, comme pour Pear/Flash2 on prend un soft qui répond au besoin. Si ça répond pas au besoin, on jette.
J'ai mis les résultats pour l'individu 1 sur vdb, dans umi_nantes sur mon espace. J'ai mis aussi un readme pour expliquer ce que contient chaque fichiers.
On a fait un point avec @mikael-s. Voici quelques observations:
Calib créé des cluster/nodes basés sur les kmer. Par défaut, le minimizer aura une valeur de 7 (pour cette longueur de reads; 150nt), et un seuil de 2. Cela provoque une sortie de calib_cons qui associait des reads très différents surtout avec les VDJ ou il suffit d'avoir un segment similaire.
En augmentant cette valeur, on retrouve bien des reads extrêmement similaire au sein d'un même cluster. D’ailleurs, nous avons maintenant 1 node pour chaque cluster, pas plus.
en regardant le résultat de calib-cons, nous avons jamais plus d'un N apparaissant. En regardant les reads initiaux, il correspondent bien à une différence d'un nt sur une position.
Si dans un cluster on a plus de 2 reads, il tentera de corriger le N en partant de ce qui est le plus représenté.
Dans le fichier de résultat de calib_cons, nous n'avons que 6 séquences avec des N (3 sur les dark base, 1 sur un UMI à la dernière position, et 2 dans la séquence).
La dernière séquence que nous ne retrouvons pas correspond à un cluster avec 4 reads, dont 1 seul présente un mismatch. Celui-ci est correctement corrigé.
Mathieu Giraudchanged title from UMI / Un pre-process pour traiter les identifiants moléculaires uniques to UMI / traitement des identifiants moléculaires uniques
changed title from UMI / Un pre-process pour traiter les identifiants moléculaires uniques to UMI / traitement des identifiants moléculaires uniques
calib serait probablement la bonne solution, merci @flothoni et @mikael-s pour votre évaluation
un server-pre-process n'est peut-être pas le point le plus urgent, extrait dans #4423 : dans un premier temps, il peut y avoir des lancements manuels et/ou intégration ailleurs, en aval du séquenceur
J'ai fait un script pour !1141 (closed), qui est maintenant extrait directement dans le repo contrib. Il s'agit d'un preprocess qui fait le demultiplexage, et qui poursuit sur un flash2.
Il faut trouver où placer ce script preprocess ou en tout cas le monter via docker pour pouvoir l'utiliser.
A voir aussi pour compiler Calib avec les options de compatibilités.