Domanda:
Unire i record del letto in base al nome
bli
2017-08-10 17:46:40 UTC
view on stackexchange narkive permalink

Ho generato un file che inizia con le seguenti righe:

  $ head -6 /tmp/bed_with_gene_ids.bedI 3746 3909 "WBGene00023193". -I 3746 3909 "WBGene00023193". -I 4118 4220 "WBGene00022277". -I 4118 4358 "WBGene00022277". -I 4118 10230 "WBGene00022277". -I 4220 4223 "WBGene00022277". -  

Vorrei unirli in base al campo nome (la quarta colonna), prendendo il minimo per l'inizio e il massimo per la fine. Gli altri campi dovrebbero essere gli stessi per tutti i record con lo stesso nome.

Risultato previsto:

  I 3746 3909 "WBGene00023193". -I 4118 10230 "WBGene00022277". -  

Ho trovato una potenziale soluzione basata su bedtools groupby qui: https://www.biostars.org/p/145751/#145775


Dati di esempio:

  cat genes.bedchr14 49894259 49895806 ENSMUST00000053290 0.000000 ... chr14 49894873 49894876 ENSMUST00000053290 0.000000. ..chr14 49894876 49895800 ENSMUST00000053291 0.000000 ... chr14 49895797 49895800 ENSMUST00000053291 0.000000 ... chr14 49901908 49901941 ENSMUST00000053291 0.000000 ...  

/ Esempio di output di blocco:

>
  sort -k4,4 genes.bed \ | groupBy -g 1,4 -c 4,2,3 -o count, min, max \ | awk -v OFS = '\ t' '{print $ 1, $ 4, $ 5, $ 2, $ 3}' chr14 49894259 49895806 ENSMUST00000053290 2chr14 49894876 49901941 ENSMUST00000053291 3  

Tuttavia:

  1. Non capisco il comportamento di groupBy (Perché -g 1,4 e non solo -g 4 ?, Perché -c 4,2,3 in questo ordine e poi riorganizzare le cose usando awk?)

  2. Questo codice non funziona per me.

Ecco cosa succede quando provo la soluzione sopra indicata:

  $ head -3 /tmp/bed_with_gene_ids.bed | bedtools groupby -g 1,4 -c 4,2,3 -o count, min, max | awk -v OFS = '\ t' '{print $ 1, $ 4, $ 5, $ 2, $ 3}' 3 3746 4220  

Ecco i tentativi basati su ciò che pensavo potesse funzionare secondo la documentazione:

  $ head -6 /tmp/bed_with_gene_ids.bed | bedtools groupby -g 4 -c 1,2,3,4,5,6 -o primo, min, max, distinto, primo, primoI 3746 10230 "WBGene00022277", "WBGene00023193". - $ head -6 /tmp/bed_with_gene_ids.bed | bedtools groupby -g 4 -c 1,2,3,4,5,6 -o primo, min, max, ultimo, primo, primoI 3746 10230 "WBGene00022277". - $ head -6 /tmp/bed_with_gene_ids.bed | bedtools groupby -g 4 -c 1,2,3,5,6 -o primo, min, max, primo, primoI 3746 10230. -  

Non capisco perché quando gruppo in base alla quarta colonna, per la quale ho due valori distinti, non riesco a ottenere due righe nell'output risultante. >

Capisco, in base ai commenti sulla pagina della documentazione, che la documentazione non è aggiornata. In particolare, esiste un'opzione -full necessaria se si desidera che tutti i campi vengano emessi. Rileggendo la soluzione sopra menzionata, penso di aver capito ora il motivo delle più colonne per l'opzione -g e per il riarrangiamento di awk . Da qui il seguente tentativo.

  $ head -6 /tmp/bed_with_gene_ids.bed | bedtools groupby -g 1,4,5,6 -c 2,3 -o min, max -full I 3746 3909 "WBGene00023193". - 3746 10230  

Ma questo non mi dà ancora due righe.

Esistono altri strumenti che potrebbero fare ciò che voglio in modo efficiente?


Modifica: soluzione

Secondo questa risposta, il problema con bedtools è che c'è un bug nell'ultima versione (2.26.0 di agosto 2017). Per avere un bedtools groupby funzionale, è necessario ottenere la versione di sviluppo da github.

Con la versione github di bedtools, ora posso ottenere il risultato atteso come segue:

  $ head -6 /tmp/bed_with_gene_ids.bed | bedtools groupby -g 1,4,5,6 -c 2,3 -o min, max | awk -v OFS = "\ t" '{print $ 1, $ 5, $ 6, $ 2, $ 3, $ 4}' I 3746 3909 "WBGene00023193". -I 4118 10230 "WBGene00022277". -  

Includo i campi 1, 5 e 6 in -g (oltre al campo 4) per poterli stampare. Nel mio file bed, dovrebbero essere gli stessi per un dato valore del campo 4. La parte awk è necessaria perché apparentemente non si ha il controllo totale sull'ordine di output: il -g i campi codice> vengono prima dei campi -c .

Cosa vuoi fare con i campi punteggio e filamento se sono diversi tra le righe o non accade mai?
In realtà, non mi interessa il campo del punteggio e idealmente lo impostare su "." se non è già così. Non posso garantire che il campo del filo sarà sempre lo stesso, ma poiché queste linee del letto provengono da annotazioni di trascrizione di cui gene_id ho inserito nel campo del nome, suppongo che in genere sarà vero che per uno stesso nome, il filo sarà lo stesso. Dovrei controllare questo, però.
Cinque risposte:
#1
+5
OrdiNeu
2017-08-10 22:40:41 UTC
view on stackexchange narkive permalink

Anche se non lo dici, immagino che tu stia usando bedtools v2.26.0. La versione 2.26.0 di groupBy contiene un bug che hai riscontrato (è stato risolto poco dopo il rilascio, quindi dovrai usare una versione prima che il bug fosse introdotto o compilare il codice sorgente corrente da https://github.com/arq5x/bedtools2)

v2.26.0:

  local10: ~ / Documents / tmp $ cat asdf. letto I 3746 3909 WBGene00023193. -I 3746 3909 WBGene00023193. -I 4118 4220 WBGene00022277. -I 4118 4358 WBGene00022277. -I 4118 10230 WBGene00022277. -I 4220 4223 WBGene00022277. -local10: ~ / Documents / tmp $ groupBy -i asdf.bed -g 4 -c 2,3 -o min, max 3746 10230 

v2.26.0-125-g52db654 (IE compilando il codice sorgente da github):

  local10: ~ / Documents / tmp $ bedtools2 / bin / groupBy -i asdf.bed -g 4 -c 2,3 -o min, maxWBGene00023193 3746 3909WBGene00022277 4118 10230  

Per rispondere alle tue domande:

1) Potresti notare che il mio output sopra fornisce prima le colonne raggruppate; dovrai riordinare l'output tramite awk per rimetterlo in ordine. Per quanto riguarda il motivo per cui hanno scelto di raggruppare su entrambe le colonne 1 e 4: se hai lo stesso nome su più cromosomi, potresti volerli trattare come caratteristiche separate.

2) Differenze di versione, come indicato nella prima parte della mia risposta.


Per unire effettivamente il file:

Assicurati di eseguirlo con una versione diversa dalla v2.26.0 (come scrive Devon Ryan nei commenti, potresti voler aggiungere colonna 6 a -g per renderlo specifico del filamento):

  ./bedtools2/bin/groupBy -i asdf.bed -g 1,4 -c 2 , 3,5,6 -o min, max, primo, primo \ | awk -v OFS = '\ t' '{print $ 1, $ 3, $ 4, $ 2, $ 5, $ 6}' I 3746 3909 WBGene00023193. -I 4118 10230 WBGene00022277. -  
Se includi 6 in `-g 1,4`, trarrai vantaggio dal non unire geni su filamenti diversi. L'UCSC a volte ha quelli e in realtà non sono lo stesso gene e non dovrebbero essere fusi insieme. Non hai bisogno di 1 in `-c`, o 6 se lo aggiungi a` -g`.
#2
+4
Ian Sudbery
2017-08-10 18:45:38 UTC
view on stackexchange narkive permalink

Puoi farlo con il CGAT toolkit:

cgat bed2bed --method = merge --merge-by-name -I bed_with_gene_ids.bed

L'installazione di un pacchetto così massiccio potrebbe essere eccessivo per questo compito.

Succede che cgat sia già installato sul mio computer (anche se ho dimenticato a che scopo). Ho provato il comando che mi suggerisci e mi ritrovo con un duplicato di `I 3746 3909" WBGene00023193 ". -`. Certo, c'erano linee duplicate nel letto originale. Ma è previsto questo comportamento?
Inoltre, se lo eseguo sull'intero file e non solo sulle prime 6 righe, dopo un po 'il programma fallisce su "TypeError:" <"non supportato tra le istanze di" Bed "e" Bed ". Sto aggiornando cgat per vedere se l'errore persiste.
Ho segnalato i problemi qui: https://github.com/CGATOxford/cgat/issues/347
Il tuo primo problema non è, per quanto ne so, il comportamento previsto. E sono abbastanza sicuro che il secondo non sia inteso. Ti suggerisco di presentare una segnalazione di bug.
Abbiamo incrociato i post!
#3
+3
Cotton Seed
2017-08-14 03:27:07 UTC
view on stackexchange narkive permalink

Puoi farlo facilmente con Hail. Hail utilizza principalmente i file BED per annotare i set di dati genetici (vedere l'ultimo esempio di annotate_variants_table), ma puoi manipolare i file BED utilizzando le funzionalità generali di Hail per manipolare i file di testo delimitati. Ad esempio:

  $ cat genes.bedI 3746 3909 "WBGene00023193". -I 3746 3909 "WBGene00023193". -I 4118 4220 "WBGene00022277". -I 4118 4358 "WBGene00022277". -I 4118 10230 "WBGene00022277". -I 4220 4223 "WBGene00022277". -  

Lo script Hail (codice python):

  from hail import * hc = HailContext () (hc .import_table ('genes.bed', impute = True, no_header = True) .aggregate_by_key ('f0 = f0, f3 = f3', 'f1 = f1.min (), f2 = f2.max (), f4 = ".", f5 = "-"' ) .select (['f0', 'f1', 'f2', 'f3', 'f4', 'f5']) .export ('genes_merged.bed', header = False))  

Il risultato:

  $ cat genes_merged.bed I 3746 3909 WBGene00023193. -I 4118 10230 WBGene00022277. -  

Aggrego su chrome e name in modo che questa soluzione non unisca voci su cromosomi diversi. Il select è necessario per riordinare i campi perché aggregate_by_key posiziona prima le chiavi da aggregare.

Divulgazione: lavoro su Hail.

#4
+2
Alex Reynolds
2017-08-10 23:48:21 UTC
view on stackexchange narkive permalink
  $ cut -f4-6 in.bed | sed 's / \ t / _ / g' | ordina | uniq | awk -F'_ '' {system ("grep" $ 1 "in.bed | bedops --merge -"); stampa $ 0; } '| incolla -d "\ t" - - | sed 's / _ / \ t / g' | sort-bed - > answer.bed  

Dato il tuo input di esempio:

  $ more in.bedI 3746 3909 "WBGene00023193". -I 3746 3909 "WBGene00023193". -I 4118 4220 "WBGene00022277". -I 4118 4358 "WBGene00022277". -I 4118 10230 "WBGene00022277". -I 4220 4223 "WBGene00022277". -  

Il file answer.bed :

  $ more answer.bedI 3746 3909 "WBGene00023193". -I 4118 10230 "WBGene00022277". -  

L'ordinamento con sort-bed è utile alla fine, in modo che tu possa collegarlo o lavorarci con altri strumenti BEDOPS o altri strumenti che ora accetta l'input BED ordinato.

Lo streaming è un modo abbastanza efficiente per fare le cose, in generale.


Come funziona

Ecco di nuovo la pipeline:

  $ cut -f4-6 in.bed | sed 's / \ t / _ / g' | ordina | uniq | awk -F'_ '' {system ("grep" $ 1 "in.bed | bedops --merge -"); stampa $ 0; } '| incolla -d "\ t" - - | sed 's / _ / \ t / g' | sort-bed - > answer.bed  

Iniziamo tagliando le colonne da 4 a 6 (id, score e strand), sostituendo le tabulazioni con trattini bassi, ordinando e rimuovendo i duplicati:

  cut -f4-6 in.bed | sed 's / \ t / _ / g' | ordina | uniq  

Ciò che otteniamo da questo è un elenco ordinato di "aghi" - uno per ciascuna combinazione di filo-punteggio ID: un ago-ID - che possiamo usare grep o filtrare il file BED originale.

Questo elenco viene reindirizzato a awk che, per ogni ago ID, esegue grep contro il file BED originale e convoglia il sottoinsieme a bedops --merge - , che unisce gli intervalli sovrapposti.

Nota che l'unione funziona solo per intervalli sovrapposti. L'unione non è necessariamente la stessa che restituire una coppia min-max e questa pipeline si interromperà se ci sono intervalli che non si sovrappongono. Ma potresti modificare l'istruzione awk in elaborare gli intervalli di input e restituire le coordinate dell'intervallo minimo e massimo, se questo è davvero ciò che si desidera, monitorando i valori minimo e massimo su tutti gli intervalli che entrano in awk e stampando un intervallo finale con un END blocco.

Il comando system stampa l'intervallo unito su una riga. La seguente istruzione print $ 0 stampa l'ago sulla riga successiva:

  awk -F'_ '' {system ("grep" $ 1 "in.bed | bedops --merge - "); stampa $ 0; } ' 

Prendiamo ogni coppia di linee alternate e le ri-linearizziamo con paste . Questo risultato ora contiene quattro colonne: le tre colonne di ogni intervallo unito e l'ago ID.

Quindi usiamo sed per sostituire i trattini bassi con tabulazioni, in modo da riportare l'ago ID in tre colonne separate da tabulazioni ID-score-strand:

  incolla -d "\ t" - - | sed 's / _ / \ t / g'  

L'output è ora un file BED a sei colonne, ma è ordinato in base all'ordinamento applicato agli aghi ID più in alto nella pipeline, che non vogliamo. Quello che vogliamo veramente è BED che è ordinato per BEDOPS sort-bed , in modo che possiamo fare più operazioni di set e ottenere un risultato corretto. Quindi lo colleghiamo a sort-bed - per scrivere un file ordinato in answer.bed:

  sort-bed - > answer. letto  
Grazie per la risposta, funziona e credo di aver capito come. Forse potrebbero essere utili alcune spiegazioni sui diversi passaggi.
#5
  0
terdon
2017-08-10 18:59:49 UTC
view on stackexchange narkive permalink

Se sei sicuro al 100% che tutto tranne le posizioni di inizio e fine saranno le stesse per tutte le linee che condividono un nome, potresti farlo da solo. Ad esempio, in Perl:

  $ perl -lane '$ start {$ F [3]} || = $ F [1]; if ($ F [1] < $ start {$ F [3]}) {$ start {$ F [3]} = $ F [1]} if ($ F [2] > $ end {$ F [3 ]}) {$ end {$ F [3]} = $ F [2]} $ chr {$ F [3]} = $ F [0]; $ rest {$ F [3]} = unisciti a "\ t", @F [4, $ # F]; END {foreach $ n (chiavi% chr) {print "$ chr {$ n} \ t $ start {$ n} \ t $ end {$ n} \ t $ n \ t $ rest {$ n}"}} 'file.bed I 3746 3909 "WBGene00023193". -I 4118 10230 "WBGene00022277". -  
Speravo che esistesse già uno strumento efficiente e che mi evitasse di reinventare la ruota in un linguaggio di scripting lento.
@bli assolutamente, questo ha molto più senso. Ho solo pensato che fosse abbastanza semplice, quindi potrei anche fornire una soluzione di scripting. Ma sì, questo sarà lento ed è anche molto ingenuo, quindi si interromperà se i tuoi file sono leggermente diversi.


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