käytä picard FilterSamReads http://broadinstitute.github.io/picard/command-line-overview.html#FilterSamReads

READ_LIST_FILE (tiedosto) -luetteloa Tiedosto, joka sisältää lukuja, jotka sisällytetään tai suljetaan pois OUTPUT SAM- tai BAM-tiedostosta. Oletusarvo: null.

(1) Se ei ole erittäin nopea, mutta voit antaa grep: lle QNAMES-tiedoston.

  samtools view file.bam | grep -f 'qnames.txt >-osajoukko.sam  

missä qnames.txt-tiedostossa on

  ESIMERKKI: QNAME1EXAMPLE: QNAME2EXAMPLE: QNAME3EXAMPLE: QNAME4EXAMPLE: QNAME5  

(2) Tämä olisi hieman monimutkaisempi, mutta voitko antaa esimerkin siitä, missä grepillä saattaa olla oikea QNAME?

(3) BAM-otsikon säilyttämiseksi Käytän samtools -H file.bam > header.txt -tunnistetta otsakkeen saamiseksi ja sitten kissaan otsikkotiedoston grep'ed-tiedostolla

Alla on pieni osa python-koodia, joka osoittaa yhden naiivin, mutta hallittavan tavan tehdä tämä (HUOM, odotan, että grep on nopeampi, vaikka sen saaminen otsikon tulostamiseen on ärsyttävää):

  import pysamqnames = set (...) # read names go herebam = pysam.AlignmentFile ("some input file.bam") obam = pysam. AlignmentFile ("jonkinlainen ulostulotiedosto.bam", "w", malli = bam) b: lle bam.fetchissä (asti_eof = Tosi): jos b.kyselyn_nimi qnames: obam.write (b) bam.close () obam. sulje ()  

Nyt se toimii, mutta jos sinulla on paljon nimiä läpi, se loppuu melko hitaasti (kuten grep -f. .. ). Jos sinun on valittava suuri luettelo luettavista nimistä, tehokkaampi strategia on:

1. Nimi lajittelee BAM-tiedoston ( samtools sort -n tai picard) , joka voidaan tehdä useilla ketjuilla.
2. Lajittele luettelo luetuista nimistä vastaavaksi (tämä on monimutkaisempaa kuin voisi naiivisti ajatella, koska samtools ja picard nimeävät lajittelun eri tavalla, joten muista laittaa jotkut ajattelivat tätä).
3. Suorita samanlainen iterointi kuin yllä, mutta vertaa vain luettelosi ylimpään luettuun nimeen vaiheesta 2 (kyseisen elementin poistaminen pinosta ottelun jälkeen).

Grep-suorituksen suorittaminen $ n $ -tasausihin $ m $ -nimiin antaisi sinulle $ O (mn) $ . Jos päätät lajitella ensin BAM- ja qnimesi, voit pienentää sen arvoon $ O (\ min \ {m, n \}) $ avaamalla luku- ja qnimet pois niiden pinot. Et myöskään lue koko BAM-tiedostoa kerralla tällä tavalla, mikä rajoittaa muistin kulutusta. Olen melko varma, että myös Picard tekee jotain tällaista.

Tässä on esimerkki Pythonista, jonka olen kirjoittanut juuri tähän tarkoitukseen:

  tuo pysamdef filter_qname ( bamfile, idfile, outfile): '' 'Suodatin lukee SAM / BAM-tiedostossa kyselynimensä perusteella Parametrit ---------- bamfile: str Lajiteltu BAM-tiedostopolku idfile: str Tekstitiedostopolku, joka sisältää qnimet säilytettäväksi outfile: str Tulostustiedosto, johon kirjoitetaan '' '# luetaan nimen mukaan lajitellussa BAM-tiedostossa bam = pysam.AlignmentFile (bamfile, "rb") jos outfile.endswith ('. sam '): output = pysam.AlignmentFile (outputfile, "w", malli = bam) elif outfile.endswith ('. bam'): output = pysam.AlignmentFile (outputfile, "wb", template = bam) else: nosta ValueError ("Tuntematon ulostulotiedostomuoto tiedostolle" outfile ": {} ". format (outfile)) # lue poistettavat tunnukset (list (set (set (...)) vain yksilöllisten tunnusten säilyttämiseksi) ids = list (set ([l.rstrip () l: lle avoimena (idfile) , 'r'). readlines ()])) # lajiteltavat tunnukset vastaamaan BAM tehokkaaseen käsittelyyn (tuhoisa toiminto) ids.sort (avain = str.lower, käänteinen = True) n_ids = len (ids) # muuttuja, jolla varmistetaan, että BAM on lajiteltu kyselyn nimen mukaan last_q = Ei lukea = bam.fetch (asti_eof = Tosi ) read = next (reads) while True: # jos luettu nimi on suurempi kuin nykyinen pinon yläosa, jos read.query_name > ids [-1]: # jos tämä on viimeinen tunnus, jos n_ids == 1: # kirjoita jäljellä olevat lukut uusi tiedosto ja irtoaa while loop output.write (lue)
r: lle bam: output.write (read) break # muuten avaa kyseinen tunnus, yritä uudelleen: ids.pop () n_ids - = 1 # ohita lukut, jotka vastaavat pinon yläosaa. 1]: lue = seuraava (lukee) # jos luettu nimi on pienempi kuin pinon yläosa, kirjoita se ja siirry eteenpäin: output.write (read) read = next (read) output.close ()  

Olen ottanut sen käyttöön omassa erilaisessa bioinformatiikkatyökalupaketissani täällä, jos haluat nähdä sen toiminnassa

Voit käyttää SAMsift:

  samsift -i file.bam -0 'q = open ("qnames .txt "). read (). splitlines () '-f' QNAME in q ' 

-i file.bam määrittää syötetyn BAM-tiedoston. -0 'q = open ("qnames.txt"). read (). splitlines ()' lataa QNAME-luettelosi. -f 'QNAME in q' määrittää kohdistuksen suodattimen.

SAMsiftin avulla voidaan käyttää mitä tahansa Python-lauseketta suodattamiseen. Haittapuoli on, että se on hitaampi kuin muut työkalut.

Suurin osa toistaiseksi lähetetyistä ratkaisuista ovat joko hitaita tai voivat tuottaa vääriä tuloksia joissakin tapauksissa. (En ole testannut Picardia, jonka uskon toimivan tarkoitetulla tavalla. Mutta henkilökohtaisesti minulla on taipumus vetäytyä, kun minun on perustettava JVM!)

Lähestymme ongelmaa eri näkökulmasta:

Hyvillä hash-funktioilla on ~ O (1) monimutkaisuus hakua varten. Jos voimme luoda QNAME-tiedostoistamme hash-taulukon, BAM-tiedoston etsimisen tulisi olla O (n) tietueiden määrässä. Python-kehittäjät ovat erittäin ylpeitä haspy-toiminnosta cpythonissa, joten käytämme sitä:

  #! / Usr / bin / env python3import syswith open (sys.argv [1], 'r') ) hakemistotiedostona: ids = set (l.rstrip ('\ r \ n') l: lle hakemistotiedostossa) riville sys.stdin: qname, _ = line.split ('\ t', 1) jos qname tunnuksissa : sys.stdout.write (rivi)  

Testisarjalla, joka koostuu ~ 3 Gt: n BAM-tiedostosta ja qnames.txt -kohdasta, jossa on ~ 65K merkintää, käynnissä samtools-näkymä input.bam | ./idfilter.py qnames.txt vie testipalvelimellani noin 6,5 sekuntia . Vertailun vuoksi SAM: n piippaaminen fgrep -f qnames.txt -palveluun (käyttäen GNU grep -ohjelmaa) kestää noin 8,5 sekuntia (mutta saattaa tuottaa vääriä tuloksia).

(Miksi grep tuottaa väärät tulokset? Jos qnames.txt -kentässäsi on QNAME foo , mutta BAM-tiedostosi sisältää myös fooX ja fooY , ne kaikki vastaavat fgrep . Ratkaisu on ankkuroida kaikki hakumallisi rivin alkuun ja välilehteen, esim. ^ foo \ t , mutta sitten sinun on käytettävä tavallista grep , jonka on rakennettava NFA jokaiselle mallille, eikä fgrep , joka toteuttaa Aho-Corasickia, ja olet suuruusluokkaa hitaampi.)

Olen kirjoittanut pienen Python-komentosarjoni uudelleen Rustiin käyttämällä std :: collections :: HashMap hash-funktiota varten ja käyttämällä rust_htslib-laatikkoa lukemaan ja kirjoittamaan suoraan BAM: sta, ja tämä on enemmän yli kaksi kertaa nopeammin kuin Python, vaikka BAM-tiedostoa luettaisiin ja kirjoitettaisiin. Rustini ei kuitenkaan todellakaan sovi julkiseen kulutukseen ...