Pitanje:
Dobivanje jedinstveno mapiranih čitanja iz BWA mem poravnanja
gringer
2017-06-07 03:23:31 UTC
view on stackexchange narkive permalink

To se temelji na pitanju betsy.s.collins na BioStarsu. Izvorni post možete pronaći ovdje.

Ima li netko prijedloge za druge oznake ili korake filtriranja BAM generiranih BAM datoteka koje se mogu koristiti tako da čita samo mapu na jedno mjesto? Primjer primjene bio bi pronalaženje sjemena za TULIP asembler / skelu, koji najbolje djeluje za čitanje te mape na jedinstvene genomske lokacije.

Koncept "jedinstveno mapirani očitavanja" učitani pojam, a većina izvora sugerira filtriranje prema MAPQ-u. Međutim, čini se da ovaj pristup ne funkcionira kada se BWA koristi kao mapper za čitanje.

Postajem pomalo drzak s ovim postom, uglavnom kako bih procijenio smatra li zajednica dobru ideju ili ne. To je pitanje koje je na BioStarsu upravo postavio potpuno novi korisnik. Pitanje je vrlo specifično, iza kojeg stoji puno dobre priče, pa mislim da se dobro uklapa u model Stack Exchange.
Jedinstveno mapirana očitanja (tj. Očitanja koja se preslikavaju na jedno mjesto u genomu) ponekad se preferiraju prilikom izvođenja analiza koje ovise o kvantifikaciji očitavanja, a ne samo o pokrivenosti (npr. RNASeq). Dvostruka očitavanja zahtijevaju dodatno filtriranje ili normalizaciju na strani analize, a većina downstream programa neće uzimati u obzir koncepte poput "petina pročitanih karata negdje * ovdje * na kromosomu 5 i četiri petine pročitanih karata negdje * ovdje * na kromosom 14. "
Ispitivač misli na genomski jedinstvena poravnanja, koja se razlikuju od više mogućih poravnanja na jednom mjestu. Ako je jedino višestruko poravnanje na istom genomskom mjestu, ta su poravnanja i dalje jedinstvena u genomskom smislu. Nemogućnost pronalaženja točnog, ispravnog poravnanja dobro je poznat problem, a postoji nekoliko nizvodnih metoda (npr. Lijeva normalizacija) kako bi se osiguralo da višestruka lokalna poravnanja i / ili pogreške u slijedu smanje utjecaj.
Hvala. Ako želite doslovno reproducirati tekst s drugog mjesta, uvijek ga stavite u navodnik (upotrijebite uređivač ili dodajte `>` na početku retka) kako biste jasno naznačili da je citiran. Ako ne, riskirate da vaš post bude izbrisan zbog plagijarizma, čak i ako citirate izvor (kao što ste to učinili).
Izvorni poster sada je odgovorio na biostars spominjući da će ići s MAPQ filtriranjem (zajedno s uklanjanjem dodatnih poravnanja).
Također je spomenula da je njezin * stvarni * problem bio u tome što Geneious i BWA mapiranja nisu isti (što je bilo razumljivo ako netko pročita u "zidu teksta" koji je objavljen); ovo čak ni izuzetak višestruko preslikanih čitanja nije riješio.
Dva odgovori:
gringer
2017-06-07 09:21:29 UTC
view on stackexchange narkive permalink

Da biste izuzeli sva moguća višestruko mapirana čitanja iz BWA preslikane BAM datoteke, čini se da trebate koristiti grep na nekomprimiranim SAM poljima:

  prikaz samtools -h mapiran.bam | grep -v -e 'XA: Z:' -e 'SA: Z:' | prikaz samtools -b > unique_mapped.bam  

Objašnjenje slijedi ...

Pretpostavit ću situaciju u kojoj se bioinformatičar prikazuje mapiranom BAM datotekom u produkciji BWA i nikako ne može doći do izvornih čitanja. Jedno rješenje s velikim naporima bilo bi izdvajanje mapiranih očitavanja iz BAM datoteke i ponovno mapiranje s drugim mapperom koji koristi MAPQ rezultat za ukazivanje na više mapiranja.

... ali što ako to ne bi bilo moguće?

Moje razumijevanje BWA-inog rezultata je da će se, ako se očitanje savršeno preslika na više genomskih lokacija, dati ocjena visoke kvalitete mapiranja (MAPQ) za obje lokacije. Mnogi ljudi očekuju da čitanje koje mapira na najmanje dvije lokacije može imati (u najboljem slučaju) 50% vjerojatnosti mapiranja na jedno od tih mjesta (tj. MAPQ = 3). Budući da BWA to ne radi, otežava filtriranje višestruko mapiranih očitavanja iz BWA rezultata pomoću MAPQ filtra koji radi za druge poravnave; ovo je vjerojatno razlog zašto trenutni odgovor na Biostars [ samtools view -bq 1 ] vjerojatno neće raditi.

Evo primjera retka iz BWA memorijskog poravnanja koji sam upravo sam napravio. To su Illumina čita mapiran parazit genoma koji ima puno ponavljaju redom:

  ERR063640.7 16 tig00019544 79974 21 21M2I56M1I20M * 0 0 TATCACATATCATCCGACTCAGCTCGACGAGTACAATGCTAATTTAACACTTAGAATGCCCGGCAATGAAATTCGTTTTCCGTCAATTCTTGAAAATTTC <AABBEGABFJKKKIM7GHKKJK>JLKLDGMHLIMIHHCGIJKKLJKLNJGLLLKLILKLMFNDLKGHJEKMKKMIJHGLOJLLLKIJLKKJEJLIGG>D NM: ja : 4 MD: Z: 83A13 AS: i: 77 XS: i: 67 XA: Z: tig00019544, -78808,21M2I56M1I20M, 6; tig00019544, -84624,79M1I20M, 6; tig00019544, -79312,33M4I42M1I20M, 8; kod> 

BWA mem je otkrio da se ovo čitanje, ERR063460.7, preslikava na najmanje tri različita mjesta: tig00019544, tig00019544 i tig00019544. Imajte na umu da je MAPQ za ovo očitanje 21, pa iako se očitane mape postavljaju na više mjesta, MAPQ se ne može koristiti za utvrđivanje.

Međutim, alternativna mjesta prikazana su prisutnošću XA u odjeljku prilagođenih polja na SAM izlazu. Možda će samo filtriranje redaka koji sadrže oznaku XA moći izuzeti višestruko preslikana čitanja. Ručna stranica pogled samtools sugerira da će -x filtrirati određenu oznaku:

  $ samtools view -x XA output.bam | grep '^ ERR063640 \ .7 [[: space:]] ERR063640.7 16 tig00019544 79.974 21 * 0 0 21M2I56M1I20M TATCACATATCATCCGACTCAGCTCGACGAGTACAATGCTAATTTAACACTTAGAATGCCCGGCAATGAAATTCGTTTTCCGTCAATTCTTGAAAATTTC <AABBEGABFJKKKIM7GHKKJK>JLKLDGMHLIMIHHCGIJKKLJKLNJGLLLKLILKLMFNDLKGHJEKMKKMIJHGLOJLLLKIJLKKJEJLIGG>D NM: i: 4 MD: Z: 83A13 AS: i: 77 XS: i: 67  

... pa je filtrirao oznaku (tj. oznaka više ne postoji u SAM izlazu), ali ne i read . U polju FLAG nema korisnih bitova koji bi ukazivali na više genomskih preslikavanja (za koja znam da se mogu filtrirati kako bi se isključilo i čitanje), pa moram pribjeći drugim mjerama.

U ovom konkretnom slučaju, Mogu koristiti grep -v na nekomprimiranom SAM izlazu kako bih izuzeo linije za poravnanje koje imaju oznaku XA (i nakon toga ponovno stisnuo u BAM, samo da bih bio uredan):

  $ samtools view -h output.bam | grep -v 'XA: Z:' | prikaz samtools -b > output_filtered.bam $ samtools prikaz output_filtered.bam | grep '^ ERR063640 \ .7 [[: space:]]' <no output>  

Ura! čita filtrirano. Kao malo po strani, ovo grep pretraživanje ima prilično veliko računarsko opterećenje: traži neki niz s tekstom XA: Z: negdje u retku, a ne zapravo zabilježiti svaku situaciju. Neka mazohistička osoba mogla bi se pojaviti kasnije i odlučiti da će nazvati sva svoja čitanja HAXXA: Z: AWESOME!: <readNumber> , u tom slučaju će biti potrebno dotjerivanje ovog grep pretraživanja kako bih bio siguran da postoji razmak (ili preciznije znak taba) prije XA:Z:.

Sada provjeravam bilo duplicirana pročitana imena, samo da budemo sigurni:

  $ samtools view output_filtered.bam | awk '{ispis $ 1}' | razvrstaj | uniq -dERR063640.1194ERR063640.1429ERR063640.1761ERR063640.2336ERR063640.2825ERR063640.3458ERR063640.4421ERR063640.4474ERR063640.4888ERR063640.49ERR063640.4974ERR063640.5070ERR063640.5130ERR063640.5300ERR063640.5868ERR063640.6116ERR063640.6198ERR063640.6468ERR063640.6717ERR063640.6797ERR063640.7322ERR063640.750ERR063640.7570ERR063640. 7900ERR063640.8115ERR063640.8405ERR063640.911ERR063640.9206ERR063640.9765ERR063640.9986  

Oh ... dovraga. Zanima me što su to:

  $ samtools view output_filtered.bam | grep '^ ERR063640.3458 [[: space:]] ERR063640.3458 16 tig00002961 5402 60 * 0 0 58S38M AGGTACCATTCGATAGAGGGAGAAAGGCACTACTAAAGATTTTGCCACATTTGCTATATCCGTATCGCGAAGATCAGGACTTACTCCGCAGAAGAA DD6HFFJBKFH = KDILKLGLJEKLKGFJIH8IKHLLMJEK: L: HBGJIHJKFLLKIHJDHLNKCK; KMKGMFKJILIIIMKI9JLKKHEJFII CC NM: I: 0 MD: Z: 38: I : 38: i XS: 0 SA: Z: tig00002377,202353 -, 14M3I5M1I35M38S, 19,5; ERR063640.3458 2064 tig00002377 202.353 19 0 0 * 14M3I5M1I35M38H AGGTACCATTCGATAGAGGGAGAAAGGCACTACTAAAGATTTTGCCACATTTGCTATA DD6HFFJBKFH = KDILKLGLJEKLKGFJIH8IKHLLMJEK: L: HBGJIHJKFLLKIHJ NM: i: 5 MD: Z : 5G48 AS: i: 35 XS: i: 27 SA: Z: tig00002961,5402, -, 58S38M, 60,0;  

Aha! Dodatna poravnanja koja koriste oznaku službena SA [ostala kanonska poravnanja u himernom poravnanju]. Čini se da su to situacije kada je jedno čitanje podijeljeno i preslikava se na više mjesta. Imajte na umu da u ovom slučaju oba poravnanja i dalje imaju MAPQ ocjene veće od 3 . Zvuči kao da bi se ispitivač također želio riješiti i ovih situacija. Ovaj put postoje i standardna polja s zastavicama za rješavanje ovih situacija (0x800: sekundarno poravnanje). Osim što nije dovoljno samo filtrirati dopunsko poravnanje, jer bi trebalo ukloniti oba pročitana preslikavanja, a ne samo ono (ili one) koje su označene kao sekundarne.

Srećom, čini se da BWA stavlja oznaku SA u sva čitanja koja sadrže dopunska poravnanja (ako to nije slučaj, siguran sam da će me netko ispraviti u vezi s tim). Dakle, u pretragu SA dodajem kao dodatni grep filter:

  $ samtools view -h output.bam | grep -v -e 'XA: Z:' -e 'SA: Z:' | prikaz samtools -b > output_filtered.bam $ samtools prikaz output_filtered.bam | awk '{ispis $ 1}' | razvrstaj | uniq -d<no output>  

Gotovo. Lako graškast! </sarcasm>

... to originalno "visoko naporno" rješenje korištenja drugog poravnavača ne izgleda sada tako loše.

to je lijepo rješenje i dobro argumentirano. Mislim da nije toliki napor ako cijelu stvar sažete u nekoliko redaka skripte ljuske. :)
`SA` označava dopunska poravnanja, a ne sekundarna poravnanja.
Stvar u korištenju MAPQ-a za filtriranje multimapera je u tome što se koristi pretpostavka da je multimapper "podjednako dobro mapiran", što ni za ovdje predstavljene primjere neće biti slučaj. Vaš sam odgovor postavio kao prihvaćeni odgovor na biostarsu, ali tamo sam zatražio da malo pojasni što se točno želi.
Da, nije teško ako znate što možete očekivati, zbog čega sam pokušao proći kroz sve probleme koji bi se mogli dogoditi.
O imenima za čitanje koji možda sadrže "XA: Z:", način postizanja robusnosti bio bi korištenje odgovarajućeg parsera SAM formata, kao što je pysam, ako se programira na pythonu, a ne grep. Ne znam koliko bi ovo bilo sporije.
Ili poput `` samtools view````, koji sam pokušao, ali nisam učinio ono što sam očekivao. Mogu zamisliti da bi samtools u budućnosti mogao implementirati nešto poput filtra oznaka. Možda već ima, a ja samo koristim pogrešan pod-alat.
U redu, dakle, još nije dostupan u samtoolsu, ali ideja o filtriranju oznaka već je pokrenuta u barem [jednom izdanju] (https://github.com/samtools/samtools/issues/665#issuecomment-293926546)
U najgorem slučaju, možete napisati nekoliko redaka pythona ili vjerojatno koristiti Pierreov jvarkit da biste filtrirali stvari.
Samir
2018-06-22 07:09:02 UTC
view on stackexchange narkive permalink

Za mnogo bržu i stabilniju implementaciju $ , sambamba to može riješiti koristeći sljedeću jednostruku liniju:

  sambamba view -t 12 -h -f bam -F "mapping_quality > = 1 a ne (nemapirano ili sekundarno_poravnanje) i ne ([XA]! = null ili [SA]! = null)" test.bwa-mem.bam -o test-uniq.bam  

Pogledajte manpage i sambamba synatx za filtre za detalje o daljnjoj upotrebi.

$: Osim grep temeljem pretraživanja koje je sporo, može vratiti lažne pozitivne negativne pogotke (gornji odgovor @ gringer). Također ne smatram pouzdan pristup zasnovan na awk jer su polja poput XA, SA itd. Neobavezna polja u SAM formatu, a nisu pozicijski fiksna polja.

Na koji su način pogotci bili lažno pozitivni?
Hmm .. Pretpostavljam da bi to trebalo biti lažno negativno jer ste izuzimali one uzorke pretraživanja s `grep -v`, zar ne?
Jedini vjerojatni način na koji bi grep mogao proizvesti pogrešne rezultate je ako je netko izmijenio pročitana imena ili grupu za čitanje kako bi sadržavao jedan od tih nizova, što se čini prilično malo vjerojatnim.
složite se, prilično je malo vjerojatno, ali bez obzira na to, pristup zasnovan na grepu u osnovi je spor.
Spori dio je u konverziji iz BAM u SAM pomoću `samtools view`. Program `grep` vrlo je brz u onome što radi.


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