Dotaz: výběr řádků textu

[Honza]

Potřeboval bych pomoct se zpravováním textových souborů. Mám dva textové soubory, kdy v jednom (CP002667.1.txt) je seznam kontigů (kusů DNA), které patří konkrétnímu bakteriálnímu druhu. Čísla kontigů jsou oddělena řádky.

Ve druhém souboru (G7_ml200contigs.fa) jsou vlastní kontigy, které ovšem patří více bakteriálním druhům. Potřebuji z něho tedy vyextrahovat jen některé kontigy, jejichž seznam je v prvním souboru.

Struktura souboru s kontigy, tedy kusy DNA patřící různým druhům je následující: první řádek je popis s číslem kontigu, následuje RÚZNÝ počet řádků s vlastní sekvencí DNA:

>contig00001 lenght=194432

ATTTAATGCCTCCACCACTAACGTTGAATCGTGCTGTTTCTGCAATCGCCAACCAACCAC TTGACGCCGCGCCAGGTTCATCACGATTGCTAAATAATGCCATTCCCCCTGTACCTTCAA

.

.

.

>contig00002 lenght=187674

TGTTCGTTATTACAATGAGGAAAGAATACAATTAAAATTAAATGGACTGAGCCCTGTTCA ATACAGAATTCAGTCCTTAAAATAACAGTCTAACTTTTCGGGGTCAGATCACAGTAATAG

.

.

.

Kdyby byl počet řádků sekvence DNA stejný, šlo by to grepem (tady dostanu první řádky sekvence):

grep -F -f CP002667.1.txt -A 2 G7_ml200contigs.fa > CP002667.1first2seqlines.txt

Jak mám ale dostat celou sekvenci kontigu, když má proměnlivý počet řádků (i několik tisíc)? Potřeboval bych nějak říct grepu, aby pokračoval s výběrem řádků za vyhledaným patternem (-A) až do dalšího hlavičkového řádku uvozeného >contig0000X.

Díky, Honza

Odpovědi

[nn]

Vstup foo.txt:

>contig00001 lenght=194432
ATTTAATGCCTCCACCACTAACGTTGAATCGTGCTGTTTCTGCAATCGCCAACCAACCAc
BBGACGCCGCGCCAGGTTCATCACGATTGCTAAATAATGCCATTCCCCCTGTACCTTCAA
CTGACGCCGCGCCAGGTTCATCACGATTGCTAAATAATGCCATTCCCCCTGTACCTTCAA
DTGACGCCGCGCCAGGTTCATCACGATTGCTAAATAATGCCATTCCCCCTGTACCTTCAA
>contig00002 lenght=187674
ETGACGCCGCGCCAGGTTCATCACGATTGCTAAATAATGCCATTCCCCCTGTACCTTCAA
FTGACGCCGCGCCAGGTTCATCACGATTGCTAAATAATGCCATTCCCCCTGTACCTTCAA
GGTTCGTTATTACAATGAGGAAAGAATACAATTAAAATTAAATGGACTGAGCCCTGTTCB
>contig00003 lenght=187674
HTGACGCCGCGCCAGGTTCATCACGATTGCTAAATAATGCCATTCCCCCTGTACCTTCAA
IGTTCGTTATTACAATGAGGAAAGAATACAATTAAAATTAAATGGACTGAGCCCTGTTCB

Regexp:

sed -n -e '0,/>contig00002/d;/>contig00003/q;p' foo.txt

Vystup:

ETGACGCCGCGCCAGGTTCATCACGATTGCTAAATAATGCCATTCCCCCTGTACCTTCAA
FTGACGCCGCGCCAGGTTCATCACGATTGCTAAATAATGCCATTCCCCCTGTACCTTCAA
GGTTCGTTATTACAATGAGGAAAGAATACAATTAAAATTAAATGGACTGAGCCCTGTTCB

[V.]

Jestli jsem to dobře pochopil ...
pokud mam čísla řádků, kde se co nachází (=klíč), tak sed by mi měl to, co je mezi tim, dokázat vyexportovat - aneb sed '2,4!d' somefile.txt. jen musim ověřit aktuální meze.

[Honza]

Příklad

Mám soubor s kontigy:

>contig00001 len=194432

ATTTAATGCCTCCACCACTAACGTTGAATCGTGCTGTTTCTGCAATCGCCAACCAACCAC

TTGACGCCGCGCCAGGTTCATCACGATTGCTAAATAATGCCATTCCCCCTGTACCTTCAA

>contig00002 len=187674

TGTTCGTTATTACAATGAGGAAAGAATACAATTAAAATTAAATGGACTGAGCCCTGTTCA

ATACAGAATTCAGTCCTTAAAATAACAGTCTAACTTTTCGGGGTCAGATCACAGTAATAG

>contig00003 len=113760

AATTATTTAATTTGCCATGTTTCTTTATTTCTGCACTTTGGCAAAATGGGCAGTTTTTTA

TTCATATTTCAAAAAGAGGGCTTAAAGCCTTGTACTATAAGGCTTTAAGCCACTTTTTAG

a potřebuji z toho dostat kontig 2, včetně hlavičkového řádku.

>contig00002 len=187674

TGTTCGTTATTACAATGAGGAAAGAATACAATTAAAATTAAATGGACTGAGCCCTGTTCA

ATACAGAATTCAGTCCTTAAAATAACAGTCTAACTTTTCGGGGTCAGATCACAGTAATAG

Problém je, že v reálu je v souboru kontigů několik set a seznam kontigů pro výběr několik desítek. Neznám rozsah řádků pro výběr.

[V.]

Takže máte ale v jednom souboru formát 1) "klíč podle kterého mohu hledat", 2) data (x řádků), 3) další klíč, podle kterého mohu hledat?
Takže "grep -n klíč zdroj" a pak extrahuju sedem od indexu k indexu.
Nebo tohle fakt nejde?

[Honza]

"grep -n klíč zdroj" funguje:

grep -F -f CP002667.1.txt -A 2 G7_ml200contigs.fa > CP002667.1first2seqlines.txt

To vypíše první dva řádky sekvence kontigů podle seznamu v souboru CP002667.1.txt. Jak ale zapojit do toho sed, aby mě vypsal všechy řádky (proměnlivý počet) sekvence za klíčovým slovem (názvem kontigu) až do hlavičkového řádku dalšího kontigu?

[Andrej]

Doporučuji vybodnout se na sed a použít awk.

[Josef Kufner]

Vybodni se na Bash a Grep. Načti si to do polí objektů v Pythonu a udělej, co potřebuješ. Ušetří ti to spoustu času až budeš chtít další trochu komplikovanější chování. Shell na tohle není stavěný. Navíc ti to otevře další možnosti vizualizace a pokročilého zpracování.

[Andrej]

Šel bych na to tak, že bych napsal awk skript:

• contig.awk:

  #!/usr/bin/awk -f
  BEGIN { output = 0
          while (getline key < key_file) ++keys[key] }
  $0 in keys { print $0; output = 1; next }
  /^>contig/ { output = 0; next }
  output { print $0 }
  

Skript by měl být spustitelný:

chmod +x contig.awk

Skript bere soubor se seznamem hlaviček jako argument (proměnnou) key_file, zatímco velký soubor, ze kterého má vyfiltrovat sekvence, bere jako standardní vstup a píše pak na standardní výstup.

Když budu mít následující dva soubory, …

• contigs.txt:

  >contig00001 len=194432
  ATTTAATGCCTCCACCACTAACGTTGAATCGTGCTGTTTCTGCAATCGCCAACCAACCAC
  TTGACGCCGCGCCAGGTTCATCACGATTGCTAAATAATGCCATTCCCCCTGTACCTTCAA
  >contig00002 len=187674
  TGTTCGTTATTACAATGAGGAAAGAATACAATTAAAATTAAATGGACTGAGCCCTGTTCA
  ATACAGAATTCAGTCCTTAAAATAACAGTCTAACTTTTCGGGGTCAGATCACAGTAATAG
  >contig00003 len=113760
  AATTATTTAATTTGCCATGTTTCTTTATTTCTGCACTTTGGCAAAATGGGCAGTTTTTTA
  TTCATATTTCAAAAAGAGGGCTTAAAGCCTTGTACTATAAGGCTTTAAGCCACTTTTTAG
  

• keys.txt:

  >contig00001 len=194432
  >contig00003 len=113760
  

…spustím to celé takto:

./contig.awk -v key_file=keys.txt < contigs.txt

To^^^ mi vypíše (celé) dvě sekvence z contigs.txt podle klíčů v keys.txt, včetně těch klíčů a bez ohledu na délku sekvencí. Pokud jsem zadání správně pochopil, tohle by mělo být ono.

Samozřejmě by sis to musel trochu upravit v případě, že formát souborů umožňuje prázdné řádky a/nebo komentáře.

[Honza]

Super, tohle dělá přesně to co potřebuju. Díky moc za pomoc. I všem ostatním.

[PetebLazar]

A pokud by to snad nekoho z mistnich zajimalo, v prikladu patrne slo mj. o bakterii "Gallibacterium anatis UMN179".

[Honza]

Bingo! Blast rulezzzz (https://blast.ncbi.nlm.nih.gov/Blast.cgi)