Domanda:
Numero di sequenze di riferimento in un file SAM
Omair Iqbal
2018-05-14 19:40:08 UTC
view on stackexchange narkive permalink

Da un singolo record, posso ottenere l'ID della sequenza di riferimento dal campo rID, ma come posso ottenere il numero totale di diverse sequenze di riferimento memorizzate nell'intero file SAM? Può essere semplice come scorrere tutti i record, ma ho bisogno di un'altra soluzione.

Con `rID` stai parlando del campo` RNAME` (nome della sequenza di riferimento)? È nell'intestazione.
rId è l'ID di riferimento corrente del record e rNextId è l'ID di riferimento della coppia di accoppiamento.
Due risposte:
Konrad Rudolph
2018-05-14 20:45:47 UTC
view on stackexchange narkive permalink

Esistono diversi modi. Eccone due:

  1. samtools idxstats ‹bamfile›
  2. samtools view -H ‹bamfile› | grep '^ @ SQ'

Entrambi i comandi forniscono un elenco dei riferimenti nel file BAM. Per ottenere il loro numero, aggiungi semplicemente | wc -l . Ma tieni presente che idxstats restituisce anche una voce finale, * , che non corrisponde a nessun riferimento ma a tutte le letture non mappate. Quindi sottrai questo dal numero.

Tuttavia, fai attenzione che mentre questo elenca tutti i riferimenti associati a un file BAM, non tutti i riferimenti hanno necessariamente dei record ad esso associati. Quindi potrebbero non esserci letture mappate a un dato riferimento. Questa è la terza colonna nell'output di idxstats . Quindi, per ottenere solo quei riferimenti con almeno una lettura mappata, possiamo filtrare per quelli:

  samtools idxstats ‹bamfile› | awk '$ 4! = 0' | head -n-1 | wc -l  

Oh, e questi (e altri) comandi funzionano sui file BAM , non su SAM. Il comando idxstats necessita inoltre dell'indice BAM, poiché in realtà sta leggendo le informazioni da quel file indice, non dal BAM (questo significa anche che funziona solo su file BAM ordinati). Lavorare su SAM richiederebbe il loop dell'intero file, il che è abbastanza inefficiente. In genere dovresti lavorare con BAM anziché con SAM, quindi si spera che questa non sia una restrizione.

Per fare lo stesso su un file SAM, devi iterare sul file:

  awk '{print $ 3}' ‹samfile› | sort -u | wc -l  

Oppure:

  awk '{refs [$ 3] = 1} END {for (ref in refs) print ref}' ‹samfile ›| wc -l  
winni2k
2018-11-21 21:48:15 UTC
view on stackexchange narkive permalink

Un file SAM memorizza gli allineamenti in una serie di sequenze di riferimento. Tuttavia, quelle sequenze di riferimento non vengono solitamente memorizzate nel file. Invece, saranno in un file di sequenza separato che di solito è in formato FASTA. Puoi ottenere il numero di sequenze di riferimento memorizzate in quel file utilizzando il comando

  grep -c '^ >' reference.fasta  


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