{"id":3357,"date":"2017-05-02T15:06:19","date_gmt":"2017-05-02T19:06:19","guid":{"rendered":"http:\/\/bioinfo.iric.ca\/?p=3357\/"},"modified":"2017-05-03T09:21:41","modified_gmt":"2017-05-03T13:21:41","slug":"flux-de-donnees-et-programmation-reactive","status":"publish","type":"post","link":"https:\/\/bioinfo.iric.ca\/fr\/flux-de-donnees-et-programmation-reactive\/","title":{"rendered":"Flux de donn\u00e9es et programmation r\u00e9active"},"content":{"rendered":"

Qu’est-ce que tout cela ?<\/h2>\n

ReactiveX est la combinaison des meilleures id\u00e9es du mod\u00e8le observateur, du mod\u00e8le it\u00e9rateur et de la programmation fonctionnelle.
\n\u00c0 l’aide des librairies Rx, vous pouvez ais\u00e9ment:
\n– Cr\u00e9er des flux de donn\u00e9es ou d’\u00e9v\u00e8nements \u00e0 partir de sources diverses comme des fichiers ou des services web
\n– Fusionner ou transformer ces flux gr\u00e2ce \u00e0 divers op\u00e9rateurs
\n– Souscrire aux flux et \u00ab\u00a0r\u00e9agir\u00a0\u00bb \u00e0 leurs \u00e9missions pour produire de nouvelles donn\u00e9es<\/p><\/blockquote>\n

L’int\u00e9r\u00eat pour la programmation r\u00e9active a progress\u00e9 ces derni\u00e8res ann\u00e9es. Peut-\u00eatre avez-vous seulement entendu le terme employ\u00e9 ici et l\u00e0, peut-\u00eatre avez-vous adopt\u00e9 ce paradigme de programmation dans votre travail quotidien… Personnellement, je n’y ai \u00e9t\u00e9 expos\u00e9 qu’au cours du dernier mois lors d’un bootcamp Angular2. J’ai imm\u00e9diatement \u00e9t\u00e9 attir\u00e9 par l’id\u00e9e de r\u00e9agir aux \u00e9v\u00e9nements \u00e9mis par des flux de donn\u00e9es.
\nNous autres, bioinformaticiens, aimons nos tuyaux (du type shell<\/em> :)) Nous canalisons constamment la sortie d’un programme vers un autre programme.
\nVous voulez savoir combien d’entr\u00e9es un fichier fasta contient? Trop facile:<\/p>\n

%> grep '>' my_fasta_file.txt | wc -l<\/code><\/pre>\n
Nous sommes donc d\u00e9j\u00e0 \u00e0 l’aise avec la notion de flux de donn\u00e9es (stdout dans ce cas) et leur consommation (avec stdin) pour produire des r\u00e9sultats.
\n\u00c0 ce point-ci, j’avais \u00e9videmment besoin de jouer avec la programmation r\u00e9active pour mieux me familiariser avec ses concepts. En recherchant un petit projet \u00e0 mettre en \u0153uvre, je me suis souvenu de cette publication de Patrick<\/a> et me suis dit: Ooooohhh je vais \u00e9crire quelque chose qui traite un fichier FastQ!
\nJustement, un projet de recherche sur lequel je travaille exige que j’identifie les fragments d’un transcriptome pouvant potentiellement coder pour un ensemble de peptides … Parfait!
\nJ’ai donc d\u00e9but\u00e9 mon p\u00e9riple en me disant que celui-ci serait assez simple. Le web fourmille d’exemples tels que \u00abAbonnez-vous \u00e0 un flux d’indice boursiers en temps r\u00e9el et affichez un avertissement si le cours d’une action a augment\u00e9 de plus de X % au cours des derni\u00e8res Y secondes, le tout avec une seule ligne de code!\u00a0\u00bb. Je devais donc atteindre mon objectif assez facilement, n’est-ce pas? Eh bien, il s’av\u00e8re que ce fut un peu plus compliqu\u00e9 que pr\u00e9vu et j’esp\u00e8re d’autres pourront b\u00e9n\u00e9ficier des solutions que j’ai trouv\u00e9es en cours de route.<\/p>\n
Contexte et sources de documentation<\/h2>\n
Premi\u00e8rement, je me dois de mentionner que la librairie ReactiveX est disponible pour de multiples langages de programmation mais puisque je suis un grand fan du langage Python, j’ai choisi de travailler avec RxPY. (Je vous laisse le soin de consulter la documentation<\/a> afin de mettre en place un environnement de d\u00e9velopmenet fonctionnel, c’est plut\u00f4t trivial).
\nIl existe de nombreuses sources de documentation sur le web mais voici celles qui m’ont le plus aid\u00e9:
\nRxPY sur github<\/a>
\nRx Marbles<\/a>
\nLe site officiel de ReactiveX<\/a>
\nLes page de documentation de Dag Brattli’s<\/a>
\nEt n’oublions pas:
\nL’introduction qu’il vous manquait<\/a> (qui est en fait une traduction de The intro you’ve been missing<\/a>)<\/p>\n
Au travail !<\/h2>\n
Alors voici le plan d’action que je propose:
\nPour chaque entr\u00e9e d’un fichier FastQ, traduire la s\u00e9quence d’ADN de l’entr\u00e9e en fonction des 6 cadres de lecture et v\u00e9rifier si l’une des 6 traductions contient un des peptides stock\u00e9es dans un fichier texte s\u00e9par\u00e9. Si l’un des peptides est trouv\u00e9, basculer l’entr\u00e9e vers un autre fichier FastQ. Que des choses simples, quoi !<\/p>\n
Lire un fichier<\/h3>\n
Le mantra de la programmation r\u00e9active est: \u00ab\u00a0Tout est un flux<\/strong>\u00a0\u00bb et RxPY offre de nombreuse fonctions utilitaires pour g\u00e9n\u00e9rer des objets Observables<\/a> (Flux) depuis diverses sources. Celle qui m’int\u00e9resse dans ce cas-ci est la fonction from_<\/code> qui construit un Observable<\/code> \u00e0 partir de n’importe quel it\u00e9rateur Python.<\/p>\n
import gzip\r\nfrom rx import Observable\r\n\r\nfastq_entries = Observable.from_(gzip.open('my.fastq.gz', 'rt'))\r\n<\/code><\/pre>\n
Comme vous pouvez le constater, les fichiers FastQ sont g\u00e9n\u00e9ralement \u00ab\u00a0gzipp\u00e9s\u00a0\u00bb car ils sont plut\u00f4t lourds.
\nMais tout ceci importe peu, nous sommes sur la bonne voie!<\/p>\n
Ajoutons quelques op\u00e9rateurs \u00e0 notre flux (Observable<\/code>) afin de retirer les retours de chariots en fin de ligne (\u00e0 l’aide de la transformation\u00a0map<\/a>) et filtrons<\/a> les lignes vides:<\/p>\n
import gzip\r\nfrom rx import Observable\r\n\r\nfastq_entries = Observable.from_(gzip.open('my.fastq.gz', 'rt')) \\\r\n    .map(lambda line: line.strip()) \\\r\n    .filter(lambda line: line != '')\r\n<\/code><\/pre>\n
Bien, bien \ud83d\ude42
\nEt soit dit en passant, nous ferons un grand usage des fonctions lambda.<\/p>\n
Les entr\u00e9es d’un fichier FastQ prennent un format particulier: elles sont toutes compos\u00e9es de 4 lignes et ressemblent \u00e0 ceci:<\/p>\n
@HWI-ST942:38:C2N4BACXX:5:1101:1228:2134 1:N:0:ATCACG\r\nTTTTCAGCCTACATCAAGGAGGTGGAGGAACGGCCGGCACCCACCCCGTGGGGCTCCAAGATGCCCTTTGGGGAACTGATGTTCGAATCCAGCAGTAGCT\r\n+\r\n@CC?DDFFHHHHGJJIJJJCHIAFHEFBACGGHGIIGIIJIIIGHHDD?CCDDBDA@DDC@C@CDDD<ACC9;7>A:?ABDB@CCCAB::AC:\r\n<\/pre>\n
Alors assurons-nous de n’\u00e9mettre que des entr\u00e9es compl\u00e8tes (ie: les 4 lignes) \u00e0 l’aide d’une transformation de type buffer<\/code><\/a>:<\/p>\n
import gzip\r\nfrom rx import Observable\r\n\r\nfastq_entries = Observable.from_(gzip.open('my.fastq.gz', 'rt')) \\\r\n    .map(lambda line: line.strip()) \\\r\n    .filter(lambda line: line != '') \\\r\n    .buffer_with_count(4)\r\n<\/code><\/pre>\n
C’est si simple ! Je commence \u00e0 tomber sous le charme des outils de cette librairie !<\/p>\n
Attaquons-nous maintenant \u00e0 la traduction des s\u00e9quences d’ADN.
\nPermettez-moi d’introduire quelques bouts de codes pour g\u00e9rer l’aspect traduction<\/em> en tant que tel..<\/p>\n
codonTable = {\r\n    'TTT': 'F', 'TCT': 'S', 'TAT': 'Y', 'TGT': 'C',\r\n    'TTC': 'F', 'TCC': 'S', 'TAC': 'Y', 'TGC': 'C',\r\n    'TTA': 'L', 'TCA': 'S', 'TAA': '*', 'TGA': '*',\r\n    'TTG': 'L', 'TCG': 'S', 'TAG': '*', 'TGG': 'W',\r\n\r\n    'CTT': 'L', 'CTC': 'L', 'CTA': 'L', 'CTG': 'L',\r\n    'CCT': 'P', 'CCC': 'P', 'CCA': 'P', 'CCG': 'P',\r\n    'CAT': 'H', 'CAC': 'H', 'CAA': 'Q', 'CAG': 'Q',\r\n    'CGT': 'R', 'CGC': 'R', 'CGA': 'R', 'CGG': 'R',\r\n\r\n    'ATT': 'I', 'ATC': 'I', 'ATA': 'I', 'ATG': 'M',\r\n    'ACT': 'T', 'ACC': 'T', 'ACA': 'T', 'ACG': 'T',\r\n    'AAT': 'N', 'AAC': 'N', 'AAA': 'K', 'AAG': 'K',\r\n    'AGT': 'S', 'AGC': 'S', 'AGA': 'R', 'AGG': 'R',\r\n\r\n    'GTT': 'V', 'GTC': 'V', 'GTA': 'V', 'GTG': 'V',\r\n    'GCT': 'A', 'GCC': 'A', 'GCA': 'A', 'GCG': 'A',\r\n    'GAT': 'D', 'GAC': 'D', 'GAA': 'E', 'GAG': 'E',\r\n    'GGT': 'G', 'GGC': 'G', 'GGA': 'G', 'GGG': 'G'\r\n}\r\n\r\n\r\ndef reverse_complement(sequence):\r\n    \"\"\"Returns the reverse complement of a DNA sequence\"\"\"\r\n    tb = str.maketrans('ACGT', 'TGCA')\r\n    return sequence.translate(tb)[::-1]\r\n\r\n\r\ndef translate_dna(fastq_entry, frame='f1'):\r\n    \"\"\"Translates DNA to amino acids according specified reading frame\"\"\"\r\n    starts = {'f1': 0, 'f2': 1, 'f3': 2,\r\n              'r1': 0, 'r2': 1, 'r3': 2}\r\n\r\n    dna = fastq_entry[1].strip('N').replace('N', 'A')\r\n\r\n    if frame.startswith('r'):\r\n        dna = reverse_complement(dna)\r\n\r\n    sequence = dna[starts[frame]:]\r\n\r\n    protein = \"\"\r\n\r\n    for i in range(0, len(sequence),  3):\r\n        if 3 == len(sequence[i:i+3]):\r\n            protein += codonTable[sequence[i:i+3]]\r\n\r\n    return protein\r\n\r\n\r\ndef translate_dna_6frames(fastq_entry):\r\n    \"\"\"Translates DNA sequence according all 6 reading frames\"\"\"\r\n    frames = ['f1', 'f2', 'f3', 'r1', 'r2', 'r3']\r\n    \r\n    translations = []\r\n    for frame in frames:\r\n        translations.append(translate_dna(fastq_entry, frame))\r\n\r\n    return translations\r\n<\/code><\/pre>\n
Je n’entrerai pas dans les d\u00e9tails de ce que font ces bouts de code, j’en ai seulement besoin pour faire un exemple focntionnel.
\nLe seul point \u00e0 noter est que je retire les ‘N’s en d\u00e9but et fin de s\u00e9quence et que je remplace les ‘N’ pr\u00e9sents dans la s\u00e9quence par des ‘A’.<\/p>\n
Tr\u00e8s bien, nous pouvons maintenant faire appel \u00e0 map<\/code> une fois de plus afin de traduire chaque entr\u00e9e FastQ:<\/p>\n
fastq_entries.map(lambda entry: translate_dna_6frames(entry))<\/code><\/pre>\n
Comment \u00e7a se passe, \u00e0 date ?<\/h3>\n
Il est probablement temps pour moi d’introduire la m\u00e9thode .subscribe()<\/code> qui permet d’observer l’ex\u00e9cution de notre programme.
\n.subscribe()<\/code> est une m\u00e9thode permettant de souscrire<\/em> un Observateur<\/code> \u00e0 notre Observable<\/code> afin de consommer ses \u00e9missions et d’y r\u00e9agir<\/strong>.
\nPour le bien de cet exemple, nous ne ferons qu’imprimer la sortie \u00e0 la ligne de commande<\/p>\n
fastq_entries.map(lambda entry: translate_dna_6frames(entry)) \r\n    .subscribe(lambda translations: print(translations)\r\n<\/code><\/pre>\n
Output<\/h4>\n
['VHYDRSGRSLGTADVHFERKADALKAMKQYNGV', 'CTMIALVAA*EQQTCTLSGRQMP*RP*SSTTA', \r\n 'AL*SLWSQLRNSRRAL*AEGRCPEGHEAVQRR', 'DAVVLLHGLQGICLPLKVHVCCS*AATRAIIVH', \r\n 'TPLYCFMAFRASAFRSKCTSAVPKLRPERS*C', 'RRCTASWPSGHLPSAQSARLLFLSCDQSDHSA']...\r\n<\/pre>\n
La premi\u00e8re chose qui m’a frapp\u00e9 lorsque j’ai commenc\u00e9 \u00e0 bidouiller avec les flux c’est qu’ils peuvent \u00eatre relativement difficiles \u00e0 observer\/d\u00e9bugger en cours d’utilisation.
\nPar exemple:\u00a0si vous ne faites pas usage de la fonction .subscribe()<\/code> pour souscrire \u00e0 vos flux \u00e0 un moment ou \u00e0 un autre dans votre code, RIEN ne se produira \u00e0 l’ex\u00e9cution du script et celui-ci fera preuve d’un mutisme parfait. Pas tr\u00e8s convivial !
\nDans un autre ordre d’id\u00e9e, il serait aussi pertinent de vous familiariser avec la notion de Schedulers<\/a> (surtout si vous comptez faire usage d’Observables \u00ab\u00a0temporaux\u00a0\u00bb tels que Observable.interval()<\/code> ou d\u00e9sirez rendre votre code asynchrone).<\/p>\n
Trouvons les peptides<\/h3>\n
Afin de trouver les peptides, nous aurons besoin de cette fonction simple:<\/p>\n
def find_peptide(peptide, translations):\r\n    for translation in translations:\r\n        if peptide in translation:\r\n            return True\r\n\r\n    return False\r\n<\/code><\/pre>\n
C’est ici que les chose se corsent :).
\nConstruire un flux \u00e0 partir de notre fichier de peptides est simple, nous savons d\u00e9j\u00e0 comment nous y prendre gr\u00e2ce \u00e0 nos tribulations avec le fichier FastQ.<\/p>\n
peptides = Observable.from_(open('peptides.txt', 'r')) \r\n    .map(lambda line: line.strip()).filter(lambda line: line != '')\r\n<\/code><\/pre>\n
Mais nous devons maintenant reproduire l’\u00e9quivalent d’une boucle imbriqu\u00e9e dans le contexte de la programmation r\u00e9active.
\nRappelez-vous, nous souhaitons identifier les fragments qui codent pour n’importe lequel de nos peptides. Hmmm..
\nVoici la solution que j’ai trouv\u00e9e:<\/p>\n
def lookup(peptides, translations):\r\n    \"\"\"Serves as the nested loop\"\"\"\r\n    found = peptides.map(lambda p: find_peptide(p, translations)) \r\n        .reduce(lambda x, y: x or y, False)\r\n\r\n    return found\r\n<\/code><\/pre>\n
lookups = fastq_entries.map(lambda entry: translate_dna_6frames(entry)) \r\n    .map(lambda translations: lookup(peptides, translations))\r\n<\/code><\/pre>\n
Alors, qu’est-ce qui se passe ici ?
\nPremi\u00e8rement, nous \u00ab\u00a0mappons\u00a0\u00bb .lookup()<\/code> sur nos chacune de nos traductions qui, \u00e0 leur tour, \u00ab\u00a0mappent\u00a0\u00bb .find_peptides()<\/code> sur notre flux de peptides. Nous reduisons<\/a> ensuite le r\u00e9sultat de ce flux \u00e0 une \u00e9mission unique prenant la forme d’un bool\u00e9en. Super, nous y sommes presque !<\/p>\n
Mais, attendez une minute.. Si nous \u00ab\u00a0souscrivons\u00a0\u00bb (avec .subscribe()<\/em>) au flux lookups<\/strong> et affichons ses \u00e9missions, il devient vite apparent que nous avons des soucis:
\n1- Notre boucle imbriqu\u00e9e ne s’ex\u00e9cute qu’une seule fois et non pour pour chaque \u00e9mission du flux translations<\/strong>
\n2- .lookup()<\/code> retourne une liste d’Observables<\/code><\/strong> plut\u00f4t que de retourner le bool\u00e9en auquel nous nous attendions (et ce, m\u00eame si .reduce()<\/code> n’\u00e9met que la valeur finale puisque cette appel est toujours encapsul\u00e9 dans un Observable<\/code>)<\/p>\n
Attaquons-nous \u00e0 r\u00e9soudre ces probl\u00e8mes un \u00e0 la fois:<\/p>\n
Probl\u00e8me 1<\/h4>\n
L’explication de ce comportement se trouve dans le fait que lorsque l’on cr\u00e9e un Observable<\/code> \u00e0 partir d’un it\u00e9rateur de fichier \u00e0 l’aide de l’utilitaire from_<\/code>, nous ne pouvons y souscrire qu’une seule fois. Toute souscription subs\u00e9quente n’\u00e9mettra aucune valeur puisque l’it\u00e9rateur a atteint la fin du fichier lors de la premi\u00e8re consommation du flux. Il existe diff\u00e9rentes fa\u00e7ons de contourner ce probl\u00e8me et l’une d’entre elle consiste \u00e0 encapsuler la construction de l’Observable dans un appel \u00e0 .defer()<\/code>.<\/p>\n
peptides = Observable.defer(\r\n    lambda: \r\n    Observable.from_(open('peptides.txt', 'r'))\r\n        .map(lambda line: line.strip())\r\n        .filter(lambda line: line != '')\r\n)\r\n<\/code><\/pre>\n
Assez simple \u00e0 r\u00e9gler me direz-vous, mais encore faut-il comprendre ce qui se passe! Et comme cet aspect n’est pas tr\u00e8s bien document\u00e9, j’ai tourn\u00e9 en rond un bon moment avant de trouver la solution.<\/p>\n
Probl\u00e8me 2<\/h4>\n
Contrairement \u00e0 .translate_dna_6frames()<\/code>, la fonction .lookup()<\/code> retourne des flux (Observables<\/em>) plut\u00f4t qu’une valeur alors nous devons trouver une fa\u00e7on de combiner cette liste de flux en un seul<\/strong> flux auquel nous pourrons alors souscrire pour r\u00e9cup\u00e9rer nos bool\u00e9ens. Il nous faut \u00ab\u00a0croiser les effluves\u00a0\u00bb ! ReactiveX expose bon nombre d’op\u00e9rateurs de combinaison<\/a> mais celui que nous utiliserons ici est switch<\/a> (ou .switch_latest(), l’impl\u00e9mentation dans RxPY)<\/p>\n
lookups = fastq_entries.map(lambda entry: translate_dna_6frames(entry)) \r\n    .map(lambda translations: lookup(peptides, translations)) \r\n    .switch_latest()\r\n<\/code><\/pre>\n
Combinaison et Multicasting<\/em><\/h3>\n
Le seule chose qui nous reste \u00e0 faire est de combiner fastq_entries<\/strong> and lookups<\/strong> en un seul flux que nous pourrons filtrer (avec .filter()<\/code>) and souscrire (avec .subscribe()<\/code>) afin de n’exporter que les entr\u00e9es FastQ qui pr\u00e9sentent le potentiel de coder pour l’un de nos peptides.
\nPour l’\u00e9tape de combinaison, j’ai choisi .zip()<\/a> (ou .zip_array()<\/a> dans RxPY).<\/p>\n
Observable.zip_array(found, fastq_entries)\r\n    .filter(lambda arr: arr[0]) \r\n    .subscribe(print)\r\n<\/code><\/pre>\n
Malheureusement un nouveau probl\u00e8me se pr\u00e9sente: lookups<\/strong> et le r\u00e9sultat de zip_array<\/code> souscrivent tous les deux au flux fastq_entries<\/strong> et re\u00e7oivent donc des entr\u00e9es FastQ distinctes, en alternance.. Pas bon \u00e7a. Nous voudrions plut\u00f4t que  zip_array()<\/code> \u00ab\u00a0zip<\/em>\u00a0\u00bb le m\u00eame<\/strong> fastq_entry qui a \u00e9t\u00e9 consomm\u00e9 et transform\u00e9 par la cha\u00eene d’action lookups<\/em>. Cette m\u00e9canique existe et porte le nom po\u00e9tique de multicasting<\/em> et nous pouvons obtenir le r\u00e9sultat escompt\u00e9 en ajoutant un appel \u00e0 .publish()<\/code> lors de la cr\u00e9ation de notre flux fastq_entries:<\/p>\n
fastq_entries = Observable.from_(gzip.open('small.fastq.gz', 'rt')) \r\n    .map(lambda line: line.strip()) \r\n    .buffer_with_count(4) \r\n    .filter(lambda line: line != '') \r\n    .publish()\r\n<\/code><\/pre>\n
Cet ajout tout simple transforme fastq_entries en un flux connectable<\/strong> qui attendra un appel \u00e0 .connect()<\/code> pour d\u00e9buter ses \u00e9missions \u00e0 tous ses souscripteurs enregistr\u00e9s.<\/p>\n
La b\u00eate dans toute sa splendeur<\/h3>\n
Voici finalement la version finale de notre script:<\/p>\n
import sys\r\nimport gzip\r\n\r\nfrom rx import Observable\r\n\r\ncodonTable = {\r\n    'TTT': 'F', 'TCT': 'S', 'TAT': 'Y', 'TGT': 'C',\r\n    'TTC': 'F', 'TCC': 'S', 'TAC': 'Y', 'TGC': 'C',\r\n    'TTA': 'L', 'TCA': 'S', 'TAA': '*', 'TGA': '*',\r\n    'TTG': 'L', 'TCG': 'S', 'TAG': '*', 'TGG': 'W',\r\n\r\n    'CTT': 'L', 'CTC': 'L', 'CTA': 'L', 'CTG': 'L',\r\n    'CCT': 'P', 'CCC': 'P', 'CCA': 'P', 'CCG': 'P',\r\n    'CAT': 'H', 'CAC': 'H', 'CAA': 'Q', 'CAG': 'Q',\r\n    'CGT': 'R', 'CGC': 'R', 'CGA': 'R', 'CGG': 'R',\r\n\r\n    'ATT': 'I', 'ATC': 'I', 'ATA': 'I', 'ATG': 'M',\r\n    'ACT': 'T', 'ACC': 'T', 'ACA': 'T', 'ACG': 'T',\r\n    'AAT': 'N', 'AAC': 'N', 'AAA': 'K', 'AAG': 'K',\r\n    'AGT': 'S', 'AGC': 'S', 'AGA': 'R', 'AGG': 'R',\r\n\r\n    'GTT': 'V', 'GTC': 'V', 'GTA': 'V', 'GTG': 'V',\r\n    'GCT': 'A', 'GCC': 'A', 'GCA': 'A', 'GCG': 'A',\r\n    'GAT': 'D', 'GAC': 'D', 'GAA': 'E', 'GAG': 'E',\r\n    'GGT': 'G', 'GGC': 'G', 'GGA': 'G', 'GGG': 'G'\r\n}\r\n\r\n\r\ndef reverse_complement(sequence):\r\n    tb = str.maketrans('ACGT', 'TGCA')\r\n    return sequence.translate(tb)[::-1]\r\n\r\n\r\ndef translate_dna(fastq_entry, frame='f1'):\r\n    \"\"\"Translates DNA to amino acids according specified reading frame\"\"\"\r\n    starts = {'f1': 0, 'f2': 1, 'f3': 2,\r\n              'r1': 0, 'r2': 1, 'r3': 2}\r\n\r\n    dna = fastq_entry[1].strip('N').replace('N', 'A')\r\n\r\n    if frame.startswith('r'):\r\n        dna = reverse_complement(dna)\r\n\r\n    sequence = dna[starts[frame]:]\r\n\r\n    protein = \"\"\r\n\r\n    for i in range(0, len(sequence),  3):\r\n        if 3 == len(sequence[i:i+3]):\r\n            protein += codonTable[sequence[i:i+3]]\r\n\r\n    return protein\r\n\r\n\r\ndef translate_dna_6frames(fastq_entry):\r\n    frames = ['f1', 'f2', 'f3', 'r1', 'r2', 'r3']\r\n\r\n    translations = []\r\n    for frame in frames:\r\n        translations.append(translate_dna(fastq_entry, frame))\r\n\r\n    return translations\r\n\r\n\r\ndef find_peptide(peptide, translations):\r\n    for translation in translations:\r\n        if peptide in translation:\r\n            return True\r\n\r\n    return False\r\n\r\n\r\ndef lookup(peptides, translations):\r\n    found = peptides.map(lambda p: find_peptide(p, translations)) \r\n        .reduce(lambda x, y: x or y, False)\r\n\r\n    return found\r\n\r\n\r\ndef main():\r\n    out = gzip.open('new.fastq.gz', 'wt')\r\n\r\n    peptides = Observable.defer(\r\n        lambda:\r\n        Observable.from_(open('peptides.txt', 'r'))\r\n            .map(lambda line: line.strip())\r\n            .filter(lambda line: line != '')\r\n    )\r\n\r\n    fastq_entries = Observable.from_(gzip.open('small.fastq.gz', 'rt')) \r\n        .map(lambda line: line.strip()) \r\n        .filter(lambda line: line != '') \r\n        .buffer_with_count(4) \r\n        .publish()\r\n\r\n    lookups = fastq_entries.map(lambda entry: translate_dna_6frames(entry)) \r\n        .map(lambda translations: lookup(peptides, translations)) \r\n        .switch_latest()\r\n    \r\n    Observable.zip_array(lookups, fastq_entries) \r\n        .filter(lambda arr: arr[0]) \r\n        .subscribe(on_next=lambda arr: print('n'.join(arr[1]), sep='', end='n', file=out, flush=True),\r\n                   on_completed=lambda: out.close())\r\n\r\n    fastq_entries.connect()\r\n\r\nif __name__ == '__main__':\r\n    sys.exit(main())\r\n<\/code><\/pre>\n
Le mot de la fin<\/h2>\n
La Programmation R\u00e9active nous promet de rendre notre code<\/p>\n