Pitanje:
Kako filtrirati unakrsne poravnave iz BED datoteke?
SmallChess
2017-05-19 10:49:47 UTC
view on stackexchange narkive permalink

Imam BAM datoteku:

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

Postoji čitanje poravnato s chr1 , a povezano je s chrx .

Imam BED datoteku:

  chr1 0 100000 TestOnly  

Želio bih filtrirati sve što spada izvan mog BED-a regija, koja uključuje unakrsno poravnanje. U mom primjeru, iako je moje čitanje usklađeno s chr1 , ali to je mate, nije. Ne želim ovo čitati.

Kada to učinim:

samtools view -L test.bed test.bam

naredba mi daje čitanje jer ne provjerava unakrsna poravnanja.

Moje rješenje:

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

ali ovo je vrlo sporo i nespretno. U svom proizvodnom cjevovodu morao bih učiniti nešto poput:

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

P: Postoji li bolje rješenje?

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

Prema specifikaciji SAM, 3. polje SAM linije ( RNAME ) je:

RNAME: Referentni niz IME poravnanja. Ako su prisutne linije zaglavlja @SQ, RNAME (ako nije ‘*’) mora biti prisutan u jednoj od SQ-SN oznaka. Nekomarirani segment bez koordinate ima '*' na ovom polju. Međutim, neomarirani segment može imati i uobičajenu koordinatu takvu da se nakon sortiranja može postaviti na željeni položaj. Ako je RNAME '*', ne mogu se pretpostaviti o POS-u i CIGARI.

A sedmo polje je (naglasak je moj, nedostaje "njihovom"):

RENEXT: Referentno ime sekvence primarnog poravnanja NEXT-a čitanog u predlošku. Za posljednje čitanje, sljedeće čitanje je prvo čitanje u predlošku. Ako su prisutne linije zaglavlja @SQ, RNEXT (ako nije ‘*’ ili ‘=’) mora biti prisutan u jednoj od SQ-SN oznaka. Ovo je polje postavljeno kao "*" kada su podaci nedostupni, i postavljeno kao "=" ako je RNEXT identično RNAME . Ako nije ‘=’, a sljedeće čitanje u predlošku ima jedno primarno mapiranje (vidi također bit 0x100 u FLAG-u), ovo je polje identično RNAME u primarnom retku sljedećeg čitanja. Ako je RNEXT '*', ne mogu se pretpostaviti PNEXT i bit 0x20

Dakle, želite ukloniti one retke čije sedmo polje nije = i za svaki slučaj oni redovi čije 7. polje nije = i nije isto što i 3. polje. Stoga možete koristiti nešto poput ovog:

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

I, da biste opet spremili kao bam datoteku:

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

U zasebnoj bilješci, vrlo rijetko je potrebno lancirati više grep naredbi na takav način. Možete ih jednostavno koristiti \ | (ili | s opcijama -E ili -P ) da biste ih razdvojili. Nešto poput:

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

Ili

  samtools view -L test.bed test.bam | grep -Ev 'chrx | chr2 | chr10 | chrN'  
Ako to učinite na ovaj način, datoteci `fixed.bam` nedostaje zaglavlje, što prema mom iskustvu stvara puno problema. Preporučujem da zaglavlje uvijek dodate; bilo specificiranjem `-h` prilikom čitanja izvornog BAM-a ili dodavanjem odvojeno:` (samtools view -H infile.bam; samtools view ...)> samtools view -b> outfile.bam`.


Ova pitanja su automatski prevedena s engleskog jezika.Izvorni sadržaj dostupan je na stackexchange-u, što zahvaljujemo na cc by-sa 3.0 licenci pod kojom se distribuira.
Loading...