Vraag:
Lees lengteverdeling uit FASTA-bestand
Scott Gigante
2017-05-17 09:38:52 UTC
view on stackexchange narkive permalink
Lees lengteverdeling uit FASTA-bestand
zie https://www.biostars.org/p/72433/; https://www.biostars.org/p/45951/
@Pierre deze oplossingen zijn niet geschikt voor langdurig lezen. Als ik een enkele (tekst) regel per 1/10 basen uitvoer, zou dit leiden tot een tekstuitvoer van meer dan 10.000 regels (en mogelijk meer dan 10x dat) voor een lange lees-fasta.
@ScottGigante zijn deze fasta in een bestand van 10 GB of een reeks bestanden met een totaal van 10 GB?
@amblina Ik heb een enkel fasta-bestand als uitvoer van poretools of iets dergelijks, met> 1M leest van gemiddelde lengte ~ 8Kb
Zeven antwoorden:
#1
+17
Sam Nicholls
2017-05-17 15:43:50 UTC
view on stackexchange narkive permalink

Als je iets snel en vies wilt, kun je de FASTA snel indexeren met samtools faidx en vervolgens de lengtekolom tot en met R (andere talen zijn beschikbaar) op de opdrachtregel plaatsen.

  samtools faidx $ fastacut -f2 $ fasta.fai | Rscript -e 'data <- as.numeric (readLines ("stdin")); samenvatting (gegevens); hist (data) ' 

Dit geeft een statistische samenvatting weer en creëert een PDF in de huidige map genaamd Rplots.pdf, met daarin een histogram.

U kunt nog eenvoudiger gaan met zoiets als: `faidx --transform chromsizes | knippen -f2 | Rscript -e 'data <- as.numeric (readLines ("stdin")); samenvatting (gegevens); hist (gegevens) ''. Dit vereist pyfaidx: `pip install pyfaidx`.
#2
+12
gringer
2017-05-17 22:33:20 UTC
view on stackexchange narkive permalink

Statistieken voor het lezen van nanoporiën zijn lastig vanwege het enorme bereik aan leeslengtes dat in een enkele run aanwezig kan zijn. Ik heb gemerkt dat de beste manier om lengtes weer te geven is door een logschaal te gebruiken op zowel de x-as (lengte) als de y-as (opeenvolgende basen of tellingen, afhankelijk van de voorkeur).

Ik heb geschreven mijn eigen scripts om dit te doen: een voor het genereren van de gelezen lengtes en een andere voor het op verschillende manieren plotten van de lengteverdeling. Het script dat leeslengtes genereert, spuugt ook de basisstatistieken voor de samenvatting van de lengte uit naar de standaardfout:

  $ ~ / scripts / fastx-length.pl > lengths_mtDNA_called.txt Totale reeksen: 2110 Totale lengte: 5,106649 Mb Langste reeks: 107.414 kb Kortste sequentie: 219 bMean Lengte: 2.42 kb Mediaan Lengte: 1.504 kb N50: 336 sequenties; L50: 3,644 kb N90: 1359 sequenties; L90: 1.103 kb $ ~ / scripts / length_plot.r lengths_mtDNA_called.txtlengths_mtDNA_called.txt ... done Aantal sequenties: 2110 Lengte kwantielen: 0% 10% 20% 30% 40% 50% 60% 70% 219,0 506,9 724,4 953,0 1196,2 1503,0 1859.2 2347.3 80% 90% 100% 3128.2 4804.7 107414.0 

Hier zijn een paar van de geproduceerde grafieken:

Digital electrophoresis plot

Read length distribution plot

De scripts om deze te genereren zijn hier te vinden:

#3
+7
xApple
2017-05-17 12:32:21 UTC
view on stackexchange narkive permalink

Het gebruik van Biopython en matplotlib lijkt inderdaad de beste manier, het komt eigenlijk neer op drie regels code om die grafiek te krijgen:

  import Bio, pandaslengths = map (len , Bio.SeqIO.parse ('/ path / to / the / seqs.fasta', 'fasta')) pandas.Series (lengths) .hist (color = 'grey', bins = 1000)  

Natuurlijk wil je misschien een langer script maken dat kan worden opgeroepen vanaf de opdrachtregel, met een paar opties. U bent van harte welkom om de mijne te gebruiken:

  #! / Usr / bin / env python2 "" "Een op maat gemaakt script om de verdeling van lengtes in een fasta-bestand uit te zetten. Geschreven door Lucas Sinclair.Kopimi. Je kunt dit script als volgt uit de shell gebruiken: $ ./fastq_length_hist --input seqs.fasta --out seqs.pdf "" "################### ################################################## ########## Modules #import argparse, sys, time, getpass, localefrom argparse import RawTextHelpFormatter van Bio import SeqIOimport panda's # Matplotlib #import matplotlibmatplotlib.use ('Agg', warn = False) van matplotlib import pyplot # ################################################## ############################ desc = "fasta_length_hist v1.0" parser = argparse.ArgumentParser (description = desc, formatter_class = RawTextHelpFormatter ) # Alle vereiste argumenten # parser.add_argument ("- input", help = "Het fasta-bestand dat moet worden verwerkt", type = str) parser.add_argument ("- out", type = str) # Alle optionele argumenten # parser.add_argument ("- x_log", default = True, typ e = bool) parser.add_argument ("- y_log", default = True, type = bool) # Parse it #args = parser.parse_args () input_path = args.inputoutput_path = args.outx_log = bool (args.x_log) y_log = bool (args.y_log) ########################################### ##################################### Lezen #lengths = map (len, SeqIO.parse ( input_path, 'fasta')) # Report # sys.stderr.write ("Lees alle lengtes (% i reeksen) \ n"% len (lengtes)) sys.stderr.write ("Langste reeks:% i bp \ n" % max (lengths)) sys.stderr.write ("Kortste reeks:% i bp \ n"% min (lengtes))
sys.stderr.write ("Making graph ... \ n") # Data #values ​​= pandas.Series (lengtes) # Plot #fig = pyplot.figure () axes = values.hist (kleur = 'grijs', bakken = 1000) fig = pyplot.gcf () title = 'Verdeling van sequentielengten' assen.set_title (titel) axes.set_xlabel ('Aantal nucleotiden in volgorde') axes.set_ylabel ('Aantal sequenties met deze lengte') assen .xaxis.grid (False) # Log # als x_log: axes.set_yscale ('symlog') if y_log: axes.set_xscale ('symlog') # Aanpassen # breedte = 18.0; hoogte = 10,0; onderkant = 0,1; top = 0,93; links = 0,07; right = 0.98fig.set_figwidth (breedte) fig.set_figheight (hoogte) fig.subplots_adjust (hspace = 0.0, bottom = bottom, top = top, left = left, right = right) # Data en bron # fig.text (0.99, 0.98, time.asctime (), horizontalalignment = 'right') fig.text (0.01, 0.98, 'user:' + getpass.getuser (), horizontalalignment = 'left') # Mooie cijfergroepering #sep = ('x' , 'y') als 'x' in sep: locale.setlocale (locale.LC_ALL, '') seperate = lambda x, pos: locale.format ("% d", x, grouping = True) axes.xaxis.set_major_formatter (matplotlib.ticker.FuncFormatter (apart)) if 'y' in sep: locale.setlocale (locale.LC_ALL, '') seperate = lambda x, pos: locale.format ("% d", x, grouping = True) axes.yaxis.set_major_formatter (matplotlib.ticker.FuncFormatter (seperate)) # Save it # fig.savefig (output_path, format = 'pdf')  

EDIT - een voorbeelduitvoer:

Sequence length distribution

#4
+5
neilfws
2017-05-17 10:09:24 UTC
view on stackexchange narkive permalink

Er zijn verschillende mogelijke benaderingen. Bijvoorbeeld:

Welke van deze zijn " snel en efficiënt "met een bestand van 10 GB ... dat is van tevoren moeilijk te zeggen. Het kan zijn dat u er een paar moet proberen te benchmarken.

#5
+4
bli
2017-05-18 14:23:18 UTC
view on stackexchange narkive permalink

bioawk zou redelijk efficiënt kunnen zijn voor dit soort taken.

  $ bioawk -c fastx '{histo [length ($ seq)] ++} END {voor (l in histo) print l, histo [l]} '\ | sort -n0 332701 15422 11323 33974 87765 118846 124747 143418 131659 1546710 2108911 3046912 4520413 6231114 8874415 11576716 14077017 19181018 31308819 51811120 109786721 472971522 657555723 273406224 101547625 49332326 32382727 16441928 10712029 7248730 4012031 2453832 2229533 1312134 938235 484736 385837 316138 285239 238840 163941 96142 68643 37744 19945 11446 7847 5948 5049 5250 4851 4252 3953 2854 4955 5956 5557 5158 5559 4360 5261 5662 4863 6764 9565488  

De -c fastx vertelt het programma om de gegevens als snel te parseren of fasta. Dit geeft toegang tot de verschillende delen van de records als $ name , $ seq (en $ qual in het geval van fastq-formaat) in de awk-code (bioawk is gebaseerd op awk, dus u kunt alle taalfuncties gebruiken die u wilt van awk).

Tussen de enkele aanhalingstekens staat een reeks <condition> {<action> } blokken.

De eerste heeft geen <condition> -deel, wat betekent dat het voor elk record wordt uitgevoerd. Hier worden de lengtetellingen bijgewerkt in een tabel die ik "histo" noemde. lengte is een vooraf gedefinieerde functie in awk.

In het tweede blok betekent de voorwaarde END dat we willen dat het wordt uitgevoerd nadat alle invoer is verwerkt. Het actiegedeelte bestaat uit het doorlopen van de geregistreerde lengtewaarden en ze samen met de bijbehorende telling af te drukken.

De uitvoer wordt doorgesluisd naar sort -n om de resultaten numeriek te sorteren.

Op mijn werkstation duurde het 20 seconden om de bovenstaande code uit te voeren voor een 1.2G fasta-bestand.

Ik realiseer me dat de uitvoer niet handig zal zijn als het gaat om schaarse lengtewaarden, omdat er geen binning is (of, equivalent, de bakken hebben breedte 1).
#6
+3
Mark Ebbert
2017-06-24 22:29:15 UTC
view on stackexchange narkive permalink

U hebt specifiek naar FASTA-bestanden gevraagd, maar het is belangrijk om altijd de lengte en kwaliteit van het lezen samen in overweging te nemen bij het beoordelen van gegevens met veel fouten en lang gelezen gegevens. FASTA-bestanden bieden niet de kwaliteit. Aan de hand van deze informatie kun je bepalen hoe succesvol de run was, hoeveel reads van 'hoge kwaliteit' waren, enz.

Ik heb hier oorspronkelijk een volledig antwoord gepost, met de suggestie pauvre, maar ik besloot dat het een beetje buiten het onderwerp lag, aangezien je blijkbaar alleen de FASTA-bestanden hebt.

Ik raad aan om FASTQ-bestanden te genereren, maar ik weet niet zeker of je de originele basis hebt -geroepen fast5 bestanden. Als dit het geval is, genereer dan FASTQ's met poretools als volgt ( poretools doc voor het genereren van FASTQ-bestanden):

  poretools fastq fast5 /  

Vervolgens raad ik aan om een heatmap en een histogrammargeplot te genereren met pauvre met zowel leeslengte als kwaliteit, zoals hieronder wordt weergegeven.

My description

[Opmerking: ik ben niet de oorspronkelijke auteur van pauvre, maar ik draag nu bij aan dit project]

Ik zou feedback op prijs stellen over waarom mijn antwoord misschien onaangenaam is.
Beste Mark, ik heb problemen met pauvre. Kun je me helpen?
#7
+1
H. Gourlé
2017-05-17 11:38:50 UTC
view on stackexchange narkive permalink

Het is niet precies wat je vroeg, maar je kunt met poretools

een histogram van de leeslengteverdeling van je nanoporiegegevens rechtstreeks vanuit de HDF5-bestanden genereren.
Vanwege schijfruimtebeperkingen creëert het standaardgedrag voor nanopore-basecallers niet langer de FAST5 (HDF5) -bestanden als uitvoer. Hoewel ik ze eigenlijk heb, zullen de meeste gebruikers dat niet doen, en bovendien zou deze methode niet generaliseren naar andere sequencers zoals PacBio.


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