Vraag:
Hoe FASTA naar BED te converteren
SmallChess
2017-05-18 05:14:58 UTC
view on stackexchange narkive permalink

Ik heb een FASTA-bestand:

  > Sequence_1GCAATGCAAGGAAGTGATGGCGGAAATAGCGTTAGATGTATGTGTAGCGGTCCC ... > Sequence_2GCAATGCAAGGAAGTGATGCGATGGAAATAGCGATGTAGTC wilCGATGATGATGTAGTC <.... bed-bestand voor elke reeks zoals:> 
  Sequence_1 0 1500Sequence_2 0 1700  

De BED-regio's hebben gewoon de grootte van de sequenties.

V: Ik deed dat eerder met een opdracht van één regel. Ik weet niet meer wat dat is, het stond op Biostars. Ik kan de post nu niet vinden. Wat is de eenvoudigste manier om de conversie uit te voeren?

Bedoel je deze biostarspost? https://www.biostars.org/p/15476/
@Greg De post geeft me geen opdracht van één regel om te doen wat ik wil. Ik ben er 100% zeker van dat ik het een tijdje geleden ergens op biostars heb gezien, maar ik herinner me de naam niet van het programma dat het kan.
Wat dacht je van deze? https://www.biostars.org/p/15476/#189089
@Greg Het was een Python-programma. Ik herinner me de naam niet en dus geen idee hoe ik moet zoeken.
Moeten we een nieuwe tag maken, misschien zoiets als "format-conversion", "file-conversion"? Zoals ik verwacht, zullen er veel meer vragen zijn van het type "X naar Y-formaat".
Zeven antwoorden:
#1
+13
bli
2017-05-18 16:33:24 UTC
view on stackexchange narkive permalink

U kunt dit gemakkelijk doen met bioawk, een versie van awk met toegevoegde functies die bio-informatica vergemakkelijken:

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

-c fastx vertelt het programma dat de gegevens moeten worden geparseerd als fasta of fastq formaat . Dit maakt de variabelen $ name en $ seq beschikbaar in de awk-opdrachten.

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

Het is een goede gewoonte om je FASTA te laten indexeren, zodat je gebruik kunt maken van de .fai die je waarschijnlijk al hebt. Als dit niet het geval is, kunt u de index gewoon genereren met samtools en wat awk gebruiken om uw BED:

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

Hiermee wordt de tabscheiding gehandhaafd, maar je kunt de BEGIN -instructie laten vallen om spaties te gebruiken. De BED-specificatie vereist alleen "witruimte" voor het eenvoudige BED -formaat.

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

Je zou deze lastige oneliner kunnen aanpassen. Merk op dat het veronderstelt dat reeks-ID's niet langer zijn dan 100 karakters en dat er geen beschrijving staat na de reeks-ID op de kopregel.

  cat myseqs.fasta | awk '$ 0 ~ ">" {print c; c = 0; printf substr ($ 0,2,100) "\ t0 \ t"; } $ 0! ~ ">" {c + = lengte ($ 0);} END {print c; } ' 

Anders biedt elke Bio * -bibliotheek (Perl, Python, Ruby) parsers in FASTA-indeling die reeks-ID's en lengtes extraheren.

Ik zou erop wijzen dat hoewel dit op BED lijkt, is het strikt genomen niet omdat BED wordt toegewezen aan coördinaten op een chromosoom of een object met een langere reeks.

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

Geïnspireerd door dit antwoord op een gerelateerde vraag over leeslengteverdelingen, zou je dit kunnen doen met Biopython:

  van Bio.SeqIO import parsewith open ("Regions .bed "," w ") als bed: voor record in parse (" Regions.fasta "," fasta "): print (record.id, 0, len (record.seq), sep =" \ t ", file = bed)  
Ik hou van Python. Ik geef de voorkeur aan dit antwoord, simpelweg omdat het in Python is en dat is wat we zouden moeten gebruiken in 2017. Maar waarom lijkt dit zo veel op Perl? Leesbaarheid telt. Ik wilde een schoner fragment voorstellen, maar het lijkt erop dat opmerkingen geen code toestaan. Zou het als een goede toon worden beschouwd om uw antwoord te bewerken?
Punt gemaakt! Ik heb het geprobeerd, maar je bent welkom om het te verbeteren. Momenteel is er niemand op de site met voldoende reputatie om een ​​bewerking te accepteren. Ik ben blij dat je een bewerking voorstelt via opmerkingen, of bewerk het bericht en wacht tot iemand 350 herhalingen bereikt. We kunnen de meta-opmerkingen verwijderen nadat de bewerking is voltooid.
OK! Ik heb het antwoord bewerkt; lijkt tenslotte een StackExchange-y manier om dingen te doen.
Ik ben ook dol op Python, maar ik heb het gevoel dat je geen script hoeft te schrijven om zoiets triviaals als dit te doen - er zijn veel opdrachtregelprogramma's die dit sneller kunnen doen.
@SamStudio8 Op het maaiveld ben ik het ermee eens. Tegelijkertijd, als een persoon een Python REPL constant open heeft (wat veel mensen zijn), heeft dit minimale overhead, en naar mijn bescheiden mening moet snel algoritmisch denken te allen tijde in dit beroep worden gepromoot. Bovendien, als deze conversie deel uitmaakt van een Snakemake-pijplijn, kan het zelfs natuurlijker zijn dan shell () aan te roepen of een shell-blok te schrijven.
@KirillG Ik heb meestal een REPL openstaan, maar ik denk dat ik nog steeds de voorkeur zou geven aan het gebruik van `awk` (of zoiets). Jaren geleden was een ander verhaal, maar het leren gebruiken van de basisprincipes van `awk` is een van de meest productiviteitsverhogende dingen die ik heb gedaan sinds ik ben overgestapt op Linux. Maar zeker, ik ben het er zeker mee eens dat als dit deel uitmaakt van een groter script of pijplijn, ik niet 'shell' zou durven noemen. Hoewel ik waarschijnlijk `pysam` zou gebruiken boven` biopython`, maar dat is slechts persoonlijke voorkeur :)
@KirillG Ik heb een BioPython-oplossing geschreven die aantoonbaar beter leesbaar / pythonisch is en niet het hele invoerbestand in het geheugen plaatst https://bioinformatics.stackexchange.com/a/106/104
@Chris_Rands'-antwoord is veel beter dan het mijne, stem alsjeblieft op!
Bedankt Scott, hoewel je antwoord met @KirillG's bewerken ook meer pythonisch is! Alleen recentere versies van BioPython staan ​​echter een string toe als invoer voor `SeqIO.parse` (in plaats van een bestandsingang) en ik denk niet dat BioPython in dit geval ooit expliciet de invoerbestandsingang sluit?
Hmm. Goed punt over de Biopython-versie - ik heb alleen de meest recente versie op Python 2.7.11 en 3.5.1 getest. Ik was ook een beetje zenuwachtig over het sluiten van bestandshandvatten, maar SeqIO.parse retourneert een generator, geen bestand, dus zou het niet automatisch het handvat sluiten als we StopIteration bereiken? (disclaimer - ik heb hier geen bewijs van!)
#5
+4
Chris_Rands
2017-05-18 13:04:25 UTC
view on stackexchange narkive permalink

Hier is een benadering met BioPython . De with -instructie zorgt ervoor dat zowel de invoer- als uitvoerbestandshandvatten zijn gesloten en dat er een luie -aanpak wordt gevolgd, zodat er slechts één fasta-record wordt vastgehouden in geheugen per keer in plaats van het hele bestand in het geheugen te lezen, wat een slecht idee is voor grote invoerbestanden. De oplossing gaat niet uit van de lengte van de reeks-ID's of het aantal regels waarover de reeksen zijn verdeeld:

  van Bio import SeqIOwith open ('sequences.fasta') als in_f, open (' sequences.bed ',' w ') as out_f: voor record in SeqIO.parse (in_f,' fasta '): out_f.write (' {} \ t0 \ t {} \ n'.format (record.id, len (record)))  
#6
+4
SmallChess
2017-05-19 10:14:21 UTC
view on stackexchange narkive permalink

We hebben veel uitstekende antwoorden! Dit zal een uitstekende referentie zijn voor toekomstige gebruikers.

Ik vond wat ik precies vroeg in mijn vraag:

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

  $ pip install pyfaidx $ faidx --transform bed test.fasta > test.bed  

Dit is het commando van één regel dat ik vroeg. De andere antwoorden werken ook, maar ik wil mijn eigen antwoord accepteren.

Ontwikkelaar van pyfaidx hier. Ik stond op het punt om dit antwoord te schrijven en betreur het dat ik geen beter uitspreekbare of gedenkwaardige pakketnaam heb gekozen.
#7
  0
gringer
2017-05-18 06:04:49 UTC
view on stackexchange narkive permalink

Als de FASTA-reeksen allemaal op één regel staan, dan zou de volgende perl one-liner moeten werken:

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

Uitleg:

  • Als de regel begint met een '>', druk dan alles af tot aan de eerste spatie (maar plaats geen regeleinde aan het einde)
  • Anders, print "0", gevolgd door de lengte van de regel, gevolgd door een regeleinde


Deze Q&A is automatisch vertaald vanuit de Engelse taal.De originele inhoud is beschikbaar op stackexchange, waarvoor we bedanken voor de cc by-sa 3.0-licentie waaronder het wordt gedistribueerd.
Loading...