Domanda

Come posso estrarre le letture da un file bam (prodotto da bwa-mem ) in fastq dato un elenco di sequenze di riferimento da filtrare?

Potenziali difficoltà

• mantenere l'orientamento FR delle letture di fine coppia (in bam tutte le sequenze sono sequenze di riferimento)
• mantenendo le letture R1 e R2
• mantenendo i punteggi di qualità nella stessa codifica dell'originale fastq (punteggi phred di illumina predefiniti nel mio caso)
• bam can essere (ana di solito è) ordinato per coordinate

Quasi ci sono soluzioni

1. soluzione perl di Sujai in blobology fa esattamente l'opposto: ottenere letture da un elenco di riferimenti (quindi potrei semplicemente invertire l'elenco). Lo svantaggio è che lo script restituisce un file fq interleaved; richiede nomi univoci di compagni, altrimenti le informazioni R1 / R2 vengono perse.

2. samtools + grep tutti da file fastq

creare un elenco di nomi di lettura che non si associano a scaffold filtrati. (cut estrarrà solo i nomi letti, uniq comprimerà i nomi di lettura delle coppie se sono uguali). Quindi grep legge i nomi dai file fastq e rimuove il separatore tra gli hit

  samtools view foo.bam | grep -vf list_of_scaffols_filter \ | cut -f 1 | uniq > list_of_reads_to_keepgrep -A 3 -f list_of_reads_to_keep foo_R1.fq | grep -v "^ - $" > foo_R1_filtered_bash.fqgrep -A 3 -f list_of_reads_to_keep foo_R2.fq | grep -v "^ - $" > foo_R1_filtered_bash.fq  

3. filter bam & picard tools

Oppure I potrei fare solo la parte di filtraggio e usare Picard-tools (Picard.SamToFastq), ma come al solito evito java il più possibile. Immagino che

  samtools view foo.bam | grep -vf list_of_scaffols_filter \ | java -jar picard.jar SamToFastq INPUT = / dev / stdin \ FASTQ = foo_R1_filtered_bash.fq SECOND_END_FASTQ = foo_R2_filtered_bash.fq  

La prima soluzione non funziona davvero per me poiché non voglio rinominare tutte le letture nel file bam e voglio mantenere le informazioni R1 / R2 (poiché R1 e R2 hanno profili di errore diversi). Entrambe le soluzioni 2 e 3 le trovo un po 'goffe e non sono sicuro che siano generali, potrei avere un comportamento imprevisto se una delle letture è mappata per secondo non lo è ... Entrambe si basano sullo stesso passaggio di filtraggio.

Mi stavo chiedendo quale fosse una soluzione pysam. Immagino che sarà molto più lento, ma almeno sarà molto più chiaro e forse più generale. Qualcosa di simile a Converti file Bam in file Fasta: c'è una soluzione pysam per fasta (non fastq), quasi lì ...

Storia

Ho genoma di riferimento molto frammentato. Alcuni scaffold sono troppo piccoli per lavorare, e alcuni di loro sono contaminanti (identificati utilizzando blobtools). Voglio separare le letture che vengono mappate a diversi gruppi per separare contaminanti, scaffold brevi e scaffold che verranno utilizzati per l'analisi a valle. Il motivo è che se rimappiamo tutte le letture a riferimento filtrato (0,7 - 0,8 del genoma originale), la maggior parte di esse (0,95 - 0,99) troverà comunque un posto dove mappano, quindi ci sono 0,2 - 0,3 di letture fuori luogo che ovviamente dovrà influire sull'analisi a valle, come l'identificazione di varianti.

Questa idea di filtraggio è basata sulla logica che se la regione genomica duplicata filtrata conterrà alcune piccole differenze, attireranno le loro letture (e se le filtro le faccio migliorerà la definizione delle varianti) e se saranno esattamente gli stessi, verranno assegnate letture casuali, quindi non c'è nulla di male nel farlo.