Domanda:
Come filtrare gli allineamenti incrociati da un file BED?
SmallChess
2017-05-19 10:49:47 UTC
view on stackexchange narkive permalink

Ho un file BAM:

  @SQ SN: chr1 LN: 248956422 @ SQ SN: chrx LN: 248956423ST-E00110: 348: HGVKKALXX: 1: 1201: 5822: 48670 323 chr1 9999 0 67H66M16H chrx 1000 0 GATAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAAC JJJJJJJJJJJJJJJJAJJJJJJJJJJJJFJJJJJJFJFJJJJJJFJJJJJJJJJJJA77FJFJJJ NM: i: 0 MD: Z: 66 AS: i: 66 XS: i: 65 SA: Z: chr5,18606834, -, 73S76M, 34,0; RG: Z: g1  

C'è una lettura allineata a chr1 , ed è mate allineata a chrx .

Ho un file BED:

  chr1 0 100000 TestOnly  

Vorrei filtrare tutto ciò che non rientra nel mio BED regione, che include allineamenti incrociati. Nel mio esempio, sebbene la mia lettura fosse allineata a chr1 ma non lo è, mate. Non voglio che venga letto.

Quando lo faccio:

samtools view -L test.bed test.bam

il comando mi dà la lettura perché non controlla gli allineamenti incrociati.

La mia soluzione:

samtools view -L test.bed test.bam | grep -v chrx

ma questo è molto lento e goffo. Nella mia pipeline di produzione dovrei fare qualcosa del tipo:

samtools view -L test.bed test.bam | grep -v chrx | grep -v ... | grep -v ... | grep -v ... | grep -v ...

D: Esiste una soluzione migliore?

Una risposta:
#1
+6
terdon
2017-05-19 22:44:29 UTC
view on stackexchange narkive permalink

Secondo la specifica SAM, il terzo campo di una riga SAM ( RNAME ) è:

RNAME: sequenza di riferimento NAME dell'allineamento. Se sono presenti righe di intestazione @SQ, RNAME (in caso contrario "*") deve essere presente in uno dei tag SQ-SN. Un segmento non mappato senza coordinate ha un "*" in questo campo. Tuttavia, un segmento non mappato può anche avere una coordinata ordinaria tale da poter essere posizionato in una posizione desiderata dopo l'ordinamento. Se RNAME è "*", non si possono fare supposizioni su POS e CIGAR.

E il settimo campo è (enfasi mia, manca "a" loro):

RNEXT: nome della sequenza di riferimento dell'allineamento principale del NEXT letto nel modello. Per l'ultima lettura, la lettura successiva è la prima lettura nel modello. Se sono presenti righe di intestazione @SQ, RNEXT (in caso contrario "*" o "=") deve essere presente in uno dei tag SQ-SN. Questo campo è impostato come "*" quando le informazioni non sono disponibili, e impostato come "=" se RNEXT è identico RNAME . Se non è "=" e la lettura successiva nel modello ha una mappatura primaria (vedere anche bit 0x100 in FLAG), questo campo è identico a RNAME nella riga principale della lettura successiva. Se RNEXT è '*', non è possibile fare supposizioni su PNEXT e bit 0x20

Quindi, si desidera rimuovere quelle righe il cui settimo campo non è = e, per ogni evenienza, quelle righe il cui settimo campo non è = e non è uguale al terzo campo. Puoi quindi usare qualcosa del genere:

  samtools view -L test.bed test.bam | awk '$ 7 == "=" || $ 3 == $ 7  

E, per salvare di nuovo come file bam:

  samtools view -L test.bed test.bam | awk '$ 7 == "=" && $ 3 == $ 7 | visualizzazione samtolls -b > fixed.bam  

In una nota a parte, è molto raro che sia necessario concatenare più comandi grep in questo modo. Puoi semplicemente usare \ | (o | con le opzioni -E o -P ) per separarli. Qualcosa come:

  samtools view -L test.bed test.bam | grep -v 'chrx \ | chr2 \ | chr10 \ | chrN'  

Oppure

  samtools view -L test.bed test.bam | grep -Ev 'chrx | chr2 | chr10 | chrN'  
Se lo fai in questo modo, al file `fixed.bam` manca l'intestazione, il che nella mia esperienza genera molti problemi. Raccomando sempre di aggiungere nuovamente l'intestazione; o specificando `-h` durante la lettura del BAM originale, o aggiungendolo separatamente:` (samtools view -H infile.bam; samtools view…)> samtools view -b> outfile.bam`.


Questa domanda e risposta è stata tradotta automaticamente dalla lingua inglese. Il contenuto originale è disponibile su stackexchange, che ringraziamo per la licenza cc by-sa 3.0 con cui è distribuito.
Loading...