Pitanje:
Kako pretvoriti FASTA u krevet
SmallChess
2017-05-18 05:14:58 UTC
view on stackexchange narkive permalink

Imam FASTA datoteku:

  > Sequence_1GCAATGCAAGGAAGTGATGGCGGAAATAGCGTTAGATGTATGTGTAGCGGTCCC ... > Sequence_2GCAATGCAAGGAAGTGATGGCGGAAATAGCGTTAGATGTATGTGTAGCGGTCCC ........  

Želim generirati BED datoteka za svaku sekvencu poput:

  Sequence_1 0 1500Sequence_2 0 1700  

BED regije jednostavno će biti veličine sekvenci.

P: To sam i prije radio s naredbom u jednom retku. Ne sjećam se što je to, bilo je na Biostarsu. Ne mogu pronaći post sada. Koji je najjednostavniji način za pretvorbu?

Mislite li na ovaj post o biostarsu? https://www.biostars.org/p/15476/
@Greg Post mi ne daje naredbu u jednom retku da radim što želim. Sigurna sam 100% da sam je negdje vidjela negdje na biostarsu, ali ne sjećam se naziva programa koji to može.
Što je s ovom? https://www.biostars.org/p/15476/#189089
@Greg Bio je to Python program. Ne sjećam se imena, a time i pojma kako pretraživati.
Trebamo li stvoriti novu oznaku, možda nešto poput "pretvorba formata", "pretvorba datoteke"? Kao što očekujem, bit će mnogo više pitanja s tipa "X u Y".
Sedam odgovori:
#1
+13
bli
2017-05-18 16:33:24 UTC
view on stackexchange narkive permalink

To možete lako učiniti s bioawk, koja je inačica awk s dodanim značajkama koje olakšavaju bioinformatiku:

  bioawk -c fastx '{print $ name " \ t0 \ t "length ($ seq)} 'test.fa  

-c fastx govori programu da podatke treba raščlaniti kao fasta ili fastq format . To čini varijable $ name i $ seq dostupnim u naredbama awk.

#2
+10
Sam Nicholls
2017-05-18 13:40:08 UTC
view on stackexchange narkive permalink

Dobra je praksa indeksirati svoj FASTA , tako da možete iskoristiti .fai koji ćete vjerojatno već imati. Ako ne, indeks možete jednostavno generirati pomoću samtools i pomoću neke awk napraviti svoj BED:

  samtools faidx $ fastaawk 'POČINI {FS = "\ t"}; {print $ 1 FS "0" FS $ 2} '$ fasta.fai > $ fasta.bed  

Ovo će zadržati razdvajanje kartica, ali možete ispustiti izjavu BEGIN koristiti prostore. BED specifikacija zahtijeva samo "razmak" za jednostavni BED format.

#3
+6
neilfws
2017-05-18 06:01:14 UTC
view on stackexchange narkive permalink

Možete prilagoditi ovaj awk one-lineer. Imajte na umu da pretpostavlja da ID-ovi sekvenci nisu duži od 100 znakova i da nema opisa koji slijedi ID-je sekvence u zaglavlju.

  cat myseqs.fasta | awk '$ 0 ~ ">" {ispis c; c = 0; printf substr (0,2 100 USD) "\ t0 \ t"; } $ 0! ~ ">" {c + = dužina ($ 0);} KRAJ {ispis c; } ' 

Inače, bilo koja knjižnica Bio * (Perl, Python, Ruby) nudi parsere formata FASTA koji će izvući ID-jeve i duljine sekvenci.

Istaknuo bih iako ovo podsjeća na BED, strogo gledano nije, jer BED preslikava koordinate na kromosom ili neki duži objekt.

#4
+4
Scott Gigante
2017-05-18 06:08:32 UTC
view on stackexchange narkive permalink

Inspirirani ovim odgovorom na povezano pitanje o distribuciji dužine čitanja, to možete učiniti pomoću Biopythona:

  iz Bio.SeqIO import parsewith open ("region .bed "," w ") kao krevet: za zapis u parseu (" region.fasta "," fasta "): ispis (record.id, 0, len (record.seq), sep =" \ t ", datoteka = krevet)  
Volim Python. Više volim ovaj odgovor jednostavno zato što je u Pythonu i to je ono što bismo trebali koristiti 2017. Ali zašto ovo toliko sliči Perlu? Čitljivost se računa. Htio sam predložiti čistiji isječak, ali izgleda da komentari ne dopuštaju kôd. Bi li se smatrao dobrim tonom za uređivanje vašeg odgovora?
Zauzet bod! Pokušao sam, ali možete to poboljšati. Trenutno na web mjestu nema nikoga s dovoljno reputacije da prihvati uređivanje. Drago mi je što ste predložili uređivanje putem komentara ili uredili post i pričekali dok netko ne pogodi 350 ponavljanja. Meta-komentare možemo izbrisati nakon što je uređivanje završeno.
U REDU! Uredio sam odgovor; čini se na kraju kao StackExchange-y način rada.
I ja volim Python, ali smatram da ne biste trebali pisati skriptu da biste učinili nešto tako trivijalno - ovo je mnogo alata naredbenog retka koji to mogu brže postići.
@SamStudio8 Na površini se slažem. Istodobno, ako osoba ima Python REPL stalno otvoren (što je puno ljudi), to ima minimalne troškove i, po mom skromnom mišljenju, u ovoj profesiji treba stalno promicati brzo algoritamsko razmišljanje. Štoviše, ako je ova pretvorba dio Snakemake cjevovoda, onda je to možda još prirodnije od pozivanja shell () ili pisanja bloka ljuske.
@KirillG Obično imam otvoren REPL, ali mislim da bih i dalje favorizirao upotrebu `awk` (ili slično). Prije nekoliko godina bila je druga priča, ali učenje kako koristiti osnove `awka` bila je jedna od stvari koja najviše poboljšava produktivnost koju sam učinio od prelaska na upotrebu Linuxa. Ali sigurno, definitivno se slažem ako je ovo dio veće skripte ili cjevovoda, ne bih se usudio nazvati `shell`. Iako bih vjerojatno koristio `pysam` umjesto` biopython`, ali to su samo osobne preferencije :)
@KirillG Napisao sam rješenje BioPython koje je uvjerljivo čitljivije / pythonic i ne stavlja cijelu ulaznu datoteku u memoriju https://bioinformatics.stackexchange.com/a/106/104
Odgovor na @Chris_Rands' je puno bolji od mog, molim vas za!
Hvala Scotte, iako je s @KirillG's edit vaš odgovor i više piton! Međutim, samo novije verzije BioPythona dopuštaju niz kao ulaz za `SeqIO.parse` (umjesto ručice datoteke) i mislim da u ovom slučaju BioPython nikada nije eksplicitno zatvorio ručicu ulazne datoteke?
Hmm. Dobra poanta oko verzije Biopython - najnoviju verziju testirao sam samo na Pythonu 2.7.11 i 3.5.1. I ja sam bio pomalo nervozan zbog zatvaranja hvataljki datoteka, ali SeqIO.parse vraća generator, a ne datoteku, pa ne bi li automatski zatvorio ručicu kad dosegnemo StopIteration? (izjava o odricanju odgovornosti - za to nemam dokaza!)
#5
+4
Chris_Rands
2017-05-18 13:04:25 UTC
view on stackexchange narkive permalink

Evo pristupa s BioPython . Izraz with osigurava da su ručke ulazne i izlazne datoteke zatvorene i da se pristupi lijenom tako da se u njemu drži samo jedan fasta zapis memorije odjednom, umjesto čitanja cijele datoteke u memoriju, što je loša ideja za velike ulazne datoteke. Rješenje ne pretpostavlja pretpostavke o duljinama ID-a sekvence ili broju linija preko kojih su sekvence raspoređene:

  iz Bio import SeqIOs open ('sequences.fasta') kao in_f, open (' sequences.bed ',' w ') kao out_f: za zapis u SeqIO.parse (in_f,' fasta '): out_f.write (' {} \ t0 \ t {} \ n'.format (record.id, len (zapis)))  
#6
+4
SmallChess
2017-05-19 10:14:21 UTC
view on stackexchange narkive permalink

Imamo mnogo izvrsnih odgovora! Ovo će biti izvrsna referenca za buduće korisnike.

Pronašao sam što sam točno pitao u svom pitanju:

https: //www.biostars. org / p / 191052 /

  $ pip install pyfaidx $ faidx --transformiraj test kreveta.fasta > test.bed  

Ovo je jednoredna naredba koju sam tražio. Ostali odgovori također rade, ali želim prihvatiti svoj odgovor.

Ovdje je programer pyfaidxa. Upravo sam htio napisati ovaj odgovor i žaliti što nisam odabrao izgovorljivije ili pamtljivije ime paketa.
#7
  0
gringer
2017-05-18 06:04:49 UTC
view on stackexchange narkive permalink

Ako su sve FASTA sekvence u jednom retku, tada bi trebala raditi sljedeća perl liner Perl:

  mačka myseqs.fasta | perl -ne 'if (/ ^ > ([^] +) /) {print $ 1} else {print "0", length, "\ n"}'  

Objašnjenje:

  • Ako redak započinje s '>', ispišite sve do prvog razmaka (ali nemojte stavljati prijelom retka na kraj)
  • U suprotnom, ispišite "0", nakon čega slijedi duljina retka, nakon čega slijedi prijelom retka


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