Domanda:
Convertire un file BAM da un riferimento a un altro?
morgantaschuk
2017-05-17 01:00:43 UTC
view on stackexchange narkive permalink

Ho una serie di file BAM che sono allineati utilizzando il riferimento del genoma umano NCBI GRCh37 (con i nomi dei cromosomi come NC_000001.10) ma voglio analizzarlo utilizzando un file BED che ha i nomi dei cromosomi hg19 UCSC (ad esempio chr1 ). Voglio usare bedtools per estrarre tutte le letture sul target e fuori target.

  1. NCBI e UCSC sono direttamente confrontabili? O devo riallineare il BAM / lift-over the BED al riferimento UCSC?
  2. Devo convertire il file BED o il file BAM? Tutti qui usano i nomi / posizioni dei cromosomi UCSC, quindi dovrò convertire comunque gli eventuali file in UCSC.
Tre risposte:
#1
+23
Devon Ryan
2017-05-17 02:04:43 UTC
view on stackexchange narkive permalink

Sei la seconda persona che abbia mai visto usare i "nomi cromosomici" NCBI (sono più simili a ID supercontig). Normalmente ti indicherei una risorsa che fornisce mappature tra i nomi dei cromosomi, ma poiché nessuno ha aggiunto nomi NCBI (eppure, forse li aggiungerò ora) al momento sei sfortunato lì.

Comunque, il modo più veloce per fare quello che vuoi è visualizzazione samtools -H foo.bam > header per ottenere l'intestazione BAM e poi cambiare ogni "nome cromosoma" NCBI con il suo nome del cromosoma UCSC corrispondente. NON RIORDINARE LE LINEE! Puoi quindi utilizzare samtools reheader e il gioco è fatto.

Perché, potresti chiedere, funzionerebbe? La risposta è che i nomi dei cromosomi / contig nei file BAM non vengono memorizzati in ogni allineamento. Piuttosto, i nomi sono memorizzati in un elenco nell'intestazione e ogni allineamento contiene solo l'indice intero in quell'elenco (gli ID del gruppo di lettura sono simili, per quello che vale). Questo porta anche all'avvertimento di cui sopra contro il riordino delle voci, dal momento che è un modo MOLTO conveniente per iniziare a scambiare gli allineamenti tra i cromosomi.

Per inciso, sarebbe molto utile passare ai nomi dei cromosomi Gencode o Ensembl, essi sono piuttosto più coerenti del pasticcio something_random presente in hg19 da UCSC.

Aggiornamento : perché sono gentile, qui è la conversione tra NCBI e UCSC. Nota che se hai allineamenti alle patch, semplicemente non esiste un equivalente UCSC. Uno dei tanti motivi per non utilizzare UCSC (evitare anche le loro annotazioni).

Quanto bene funziona? Hai eseguito benchmark? Lo chiedo perché ho provato a convertire vari file bed tra genomi hg e GRC e i tre strumenti che ho usato hanno dato tutti risultati molto diversi. Questo tipo di mappatura dovrebbe essere davvero una cosa semplice, ma sembra non essere così semplice.
Per i casi in cui è solo un cambio di nome (la maggior parte dei casi) non c'è nulla da confrontare. Per i casi in cui hai anche cambiamenti di posizione, avresti bisogno di una risorsa diversa (vale a dire, liftOver o crossmap).
Sì, sono pedaggi come liftOver e crossmap che ho usato e con cui ho riscontrato problemi. Mi aspettavo che questo fosse un problema risolto, ma sfortunatamente ognuno dei tre strumenti che ho usato ha dato risultati diversi. Il che mi rende diffidente nell'usare i risultati.
I risultati sono deterministici, dovrebbero essere gli stessi indipendentemente dallo strumento a condizione che si utilizzino le stesse impostazioni.
Lo penseresti, sì. Ma non lo sono. Sarei felice di chattare ([chat]) su questo (questa è una delle cose che devo risolvere al lavoro) se lo desideri, e sarei molto contento di apprendere che stavo semplicemente sbagliando, ma il mio test preliminari e vari thread che ho letto online (vedi [qui] (https://www.biostars.org/p/14187/#91490) per esempio) suggeriscono che non è così semplice come potresti pensare.
Ho il sospetto che questo dipenda da carenze nei file a catena, ma mi piacerebbe sentirlo.
#2
+3
chrisamiller
2017-05-17 23:56:14 UTC
view on stackexchange narkive permalink

La soluzione "giusta" sarebbe il riallineamento, ma è costoso e la maggior parte di noi non seguirebbe quella strada. La mia soluzione preferita sarebbe convertire il file bed, al contrario del bam. Ecco perché:

1) Reintestare il bam significa che potresti avere letture allineate a contig senza una voce corrispondente in UCSC (vedi l'elenco di Devon per le mappature). Questo è un problema perché:

  • Alcune di quelle letture sarebbero state probabilmente mappate altrove se fosse stato usato un riferimento senza quei contig.
  • Non sono nemmeno sicuro di cosa accadrà a quelle letture dopo la reintestazione: immagino che dovrebbero essere contrassegnate come non mappate? C'è un sacco di potenziale per le cazzate lì.

2) Sembra più pulito convertire il file bed da UCSC-> NCBI, dove hai la garanzia che ogni voce ha una "casa". Quindi, dopo aver estratto le informazioni dal bam, puoi sempre riconvertire i nomi dei cromosomi se necessario.

Gli allineamenti che finiscono per perdere il loro cromosoma finiscono per non essere mappati, sebbene non ottengano il flag SAM non mappato (piuttosto, finiscono per essere rappresentati come se fossero sul cromosoma "*"). Ma sì, "samtools reheader" dovrebbe essere usato con estrema cautela.
Diversi strumenti daranno un tremolio se passi le letture su "*" che non sono contrassegnate come non mappate. "pysam" è uno. Lo so per esperienza amara (correzione di bug).
#3
  0
Ian Sudbery
2017-08-08 16:02:03 UTC
view on stackexchange narkive permalink

Solo per sottolineare che se vuoi seguire la risposta di @Devon Ryan per un diverso organismo / assemblaggio, che non è nella sua risorsa collegata molto utile, puoi scaricare NCBI in UCSC contig per mappature del numero di cromosomi da https : //www.ncbi.nlm.nih.gov/assembly.

Vai al sito e cerca il nome dell'assembly. In fondo alla pagina c'è una casella chiamata "Definizione dell'assembly globale" contenente un collegamento intitolato "Scarica il rapporto della sequenza completa".

Il file scaricato contiene una tabella con:

  • Numeri dei cromosomi in "Sequence-Name" / "Assigned-molecule"
  • Nomi NCBI in Refseq-Accn
  • Nome contig UCSC in "UCSC-style-name"


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...