Domanda:
Sottoinsieme BAM più piccolo per contenere diverse migliaia di righe da più cromosomi
EB2127
2018-02-17 18:00:18 UTC
view on stackexchange narkive permalink

Ci sono molti casi in cui vorrei creare un sottoinsieme BAM per creare un piccolo file con cui lavorare (es. test algoritmico, debugging, ecc.)

Normalmente faccio quanto segue, che sarà subset il BAM file.bam e mantieni l'intestazione

  samtools view -H file.bam > header.samsamtools view file.bam | testa -n 5000 | cat header.sam - | samtools view -Sb - > file.unique.bam  

In questo caso, vorrei 5000 righe nel cromosoma 1 e 5000 nel cromosoma 2.

Potrei prima prova il grepping per singolo cromosoma e poi combina i due SAM

es ecco BAM completo con grepped chr1 e intestazione (errata ma completa)

  samtools view -H file.bam > header.samsamtools view file.bam | grep "chr1" | cat header.sam - | samtools view -Sb - > file.unique.bam  

ma poi ho due problemi:

(1) Potrei non ingannare gli allineamenti al cromsoma 2- - potrebbero esistere righe BAM che contengono "chr2" ma non sono allineamenti.

(2) Penso che sia necessario modificare manualmente l'intestazione. Probabilmente non c'è modo di aggirare questo.

Esiste un modo semplice, Bioinformatics SO?

Cinque risposte:
#1
+3
Devon Ryan
2018-02-17 18:29:12 UTC
view on stackexchange narkive permalink

Se non sei troppo bloccato su numeri esatti come 5000 letture, puoi farlo con un singolo comando samtools:

  samtools view -bo subset.bam -s 123.4 alignments. bam chr1 chr2  

Questo selezionerà il 40% (la parte .4 ) delle letture ( 123 è un seme, che è conveniente per la riproducibilità). La parte conveniente di questo è che manterrà gli accoppiamenti accoppiati se hai letture accoppiate. Per 5000 letture per cromosoma, cambia semplicemente la parte .4 in un numero sufficientemente piccolo.

In generale non hai davvero bisogno di sottoinsiemi l'intestazione. Alcuni strumenti funzioneranno un po 'meglio se lo fai, ma generalmente otterrai gli stessi risultati a prescindere.

#2
+3
Karel Brinda
2018-02-27 04:30:15 UTC
view on stackexchange narkive permalink

Puoi utilizzare SAMsift:

  samsift \ -i file.bam \ -0 'c = {"chr1": 5000, "chr2": 5000 } '\ -f' c [RNAME] >0 '\ -c' c [RNAME] - = 1 '\ -m nonstop-remove  

Spiegazione:

  • -i file.bam - file di input
  • -0 'c = {"chr1": 5000, "chr2": 5000}' - inizializzazione (crea un dizionario per il conto alla rovescia per i cromosomi di interesse)
  • -f 'c [RNAME] >0' - criterio di filtraggio (è il contatore del cromosoma corrente ancora> 0?)
  • -c 'c [RNAME] - = 1' - codice che decrementa il contatore del cromosoma corrente (5000 -> 4999 -> ... - > 1 -> 0 -> -1 -> ...)
  • -m nonstop-remove - rimuove le righe che causano errori Python e non fermarsi (in questo caso, un errore può essere causato dall'accesso a un contatore inesistente per un altro cromosoma, ad esempio per chr3)

Per ulteriori informazioni, vedere il Leggimi di SAMsift.

Se potessi aggiungere qualche spiegazione su quali sono queste opzioni e perché sono necessarie. Ad esempio, non capisco perché l'argomento di "f" sia "c [RNAME]" e non altre cose
Ho appena aggiunto la spiegazione.
#3
+2
Ian Sudbery
2018-02-17 21:20:09 UTC
view on stackexchange narkive permalink

In genere consiglierei la risposta di Devon Ryan. Tuttavia, se ti interessa avere lo stesso numero per ogni Chr, puoi usare il seguente codice python / pysam (questo produrrà circa 5000 da ogni Chr):

  from pysam import AlignmentFilefrom random import randomnreads = 5000infile = AlignmentFile (" mybam.bam ") outfile = AlignmentFile (" outbam.bam "," wb ", template = infile) per chr in infile.reference_names: reads_in_chr = infile.count (chr) frac = float (nreads) / reads_in_chr count = 0 # - Sostituisci questo ciclo per esattamente 5000 prime letture - # per read in infile.fetch (chr, multiple_iterators = True): if random () < frac: count + = 1 outfile.write (read) print ("ouputted% i reads from% s"% (count, chr)) outfile.close ()  

Questo usa tutto cromosomi. Se vuoi usare solo, dì chr1 e chr2 sostituisci infile.reference_names con ["chr1","chr2"”.

Se vuoi esattamente 5000 letture da ogni chr, ma non importa se sono le prime o meno, allora potresti sostituire il ciclo for interno con:

  da leggere in infile.fetch (chr, multiple_iterators = True): count + = 1 if count > = nreads: break outfile.write (read)  

se volevi i compagni di quelli si legge, puoi aggiungere quanto segue dopo outfile.write (read) :

  mate_read = infile.mate ( read) outfile.write (mate_read)  

Nota che è lento.

Per velocizzare le cose, `infile.count (chr)` può essere sostituito con l'analisi di `idxstats ()` (a meno che non si stiano usando file CRAM).
Grazie, pensavo che ci sarebbe stato un modo semplice per recuperarlo dall'indice.
"Per velocizzare le cose, infile.count (chr) può essere sostituito con l'analisi di idxstats ()" Non lo sto seguendo. Come sarebbe il precedente?
#4
+1
Pierre
2018-02-19 19:49:23 UTC
view on stackexchange narkive permalink

utilizzando samjdk : http://lindenb.github.io/jvarkit/SamJdk.html

  $ java -jar dist /samjdk.jar --body -e \ 'Map<String, Integer> c = new HashMap<> (); public Object apply (SAMRecord r) {int n = c.getOrDefault (r.getContig (), 0); if (n> = 5000) restituisce false; c.put (r.getContig (), n + 1); return r;} '\ input.bam  
Sarebbe bello se le persone che votano per difetto spiegassero perché lo stanno facendo ...
@Devon (ho visto il commento fino ad ora) Ho downvoted perché trovo che una lunga riga di codice senza spiegare cosa fa internamente è di solito più difficile da capire, il codice così com'è migliorerebbe molto (IMHO) suddiviso in più righe e con commenti . È difficile per me che non conosco Java per capire cosa sta succedendo
#5
+1
Konrad Rudolph
2018-02-20 21:08:42 UTC
view on stackexchange narkive permalink

Se desideri solo creare un file BAM troncato con un'intestazione, puoi semplificare notevolmente il codice originale:

  (samtools view -H input.bam; samtools visualizza input.bam | head -5000) | samtools -bo output.bam  

Questo comando evita il file intermedio e poche invocazioni di comandi intermedie gratuite, al costo di una subshell (invocato da (…) ) .

Come sopra, ma con formattazione:

  (samtools view -H input.bam samtools view input.bam | head -5000 # (*)) \ | samtools -bo output.bam  

... questo può ovviamente essere esteso per filtrare da più cromosomi sostituendo la riga contrassegnata con (*) sopra con uno o più righe che subset dal nome del cromosoma ( samtools view input.bam chr x , non c'è bisogno di grep se hai indicizzato il file BAM originale!) .

È "samtools -bo"?


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