Vraag:
Hoe NGS-uitlezingen te simuleren, waarbij de reeksdekking wordt gecontroleerd?
SmallChess
2017-05-17 21:26:17 UTC
view on stackexchange narkive permalink

Ik heb een FASTA-bestand met 100+ reeksen zoals deze:

  >Sequence1GTGCCTATTGCTACTAAAA ... >Sequence2GCAATGCAAGGAAGTGATGGCGGAAATAGCGTTA ......  

Ik heb ook een tekstbestand zoals dit:

  Sequence1 40Sequence2 30 ......  

Ik zou de volgende generatie gepaarde-end reads willen simuleren voor alle de sequenties in mijn FASTA-bestand. Voor Sequence1 zou ik willen simuleren met 40x dekking. Voor Sequence2 zou ik willen simuleren met 30x dekking. Met andere woorden, ik wil mijn reeksdekking voor elke reeks in mijn simulatie regelen.

V: Wat is de eenvoudigste manier om dat te doen? Welke software moet ik gebruiken? Bioconductor?

[Deze recensie] (http://www.nature.com/nrg/journal/v17/n8/abs/nrg.2016.57.html) zou een goed startpunt kunnen zijn. Het vergelijkt 23 NGS-simulatietools en biedt een beslissingsboom voor het selecteren van een geschikt hulpmiddel op basis van uw behoeften.
Wat is de leeslengte die u gebruikt? Hoe lang zijn de reeksen? Moet u het dekkingsdoel precies of met enige waarschijnlijkheid bereiken?
Klinkt als iets dat je gewoon in tien of twintig regels python zou coderen met de Biopython-module.
Ik zou nog een paar vragen aan Greg's willen toevoegen. Begrijp ik goed dat u de volgorde van de sjabloonreeks vanuit het bestand wilt simuleren? Dus vertegenwoordigen de sequenties genomen? Zou u versterkingsbias willen overwegen? Opeenvolgende fouten? Welk sequencing-platform zou je willen simuleren?
Negen antwoorden:
#1
+6
H. Gourlé
2017-05-30 12:06:03 UTC
view on stackexchange narkive permalink

Ik werk aan een Illumina-sequencing-simulator voor metagenomics: InSilicoSeq

Het is nog steeds in alpha-release en erg experimenteel, maar gezien een multi-fasta en een overvloed-bestand, het genereert reads van uw ingevoerde genomen met verschillende dekkingen.

Uit de documentatie:

  iss genereren --genomes genomes.fasta --abundance abundance_file.txt \ - model_file HiSeq2500 --output HiSeq_reads  

Waar:

  # multi-fasta bestand >genome_AATGC ... >genome_BCCGT ...... # overvloed bestand (totale overvloed moet 1 zijn!) genome_A 0.2genome_B 0.4 ...  

Ik heb het niet ontworpen om te werken met dekking maar eerder een overvloed van het genoom in een metagenoom, dus je moet misschien een een klein beetje wiskunde;)

#2
+4
Ian Sudbery
2017-05-18 13:57:06 UTC
view on stackexchange narkive permalink

Het bioconductorpakket van polyester kan dit doen. Er staat dat het RNA-seq-reads simuleert, maar ik weet niet of dat echt anders is dan andere NGS-reads.

Het kan een reeks fout- en vertekeningsmodellen gebruiken, of ze leren van een dataset.

Kunt u dit een beetje toelichten? Hoe zou het OP het gebruiken voor hun bestanden? Welk bereik van fouten en vooringenomenheid? Hoe zouden we ze gebruiken? Wat zou een redelijke standaard zijn om te gebruiken? Misschien zou een minimaal werkend voorbeeld nuttig zijn.
#3
+4
Kamil S Jaron
2017-05-18 22:32:11 UTC
view on stackexchange narkive permalink

Dit python-script neemt een fasta-bestand en een tsv-bestand met tellingen en drukt de reeksen in de fasta-bestanden zo vaak af als het is gespecificeerd in het tsv-bestand (uitgaande van het formaat in de vraag). Dus als bar.tsv en foo.fasta uw bestanden zijn:

  van Bio import SeqIOrepeat = {} voor regel in open ( "bar.tsv"): seq_id, coverage = line.split () herhaal [seq_id] = int (dekking) voor seq_record in SeqIO.parse (foo.fasta, "fasta"): voor i binnen bereik (repeat.get ( seq_record.name, 0)): print (">", seq_record.name, "_", i, sep = '') print (seq_record.seq)  
Kunt u uitleggen wat dit doet en hoe het werkt? Het lijkt erop dat u alleen de oorspronkelijke dekkingstijden van de reeks $ afdrukt, en dat kan niet kloppen. Het OP moet lezingen simuleren, niet alleen dezelfde reeks N keer herhalen.
Au, toen heb ik de vraag verkeerd begrepen. Ik dacht dat de sequenties de reads zijn die moeten worden gesimuleerd en 40x is lateraal, ik zou deze sequentie 40 keer willen hebben. Ik heb het antwoord bewerkt om in ieder geval verwarring te voorkomen.
Ik weet niet of je het verkeerd hebt begrepen, of ik deed het. Maar een van ons deed het :) Bedankt voor de bewerking, nu kan het OP tenminste zien wat dit doet en zelf beslissen.
#4
+4
Daniel Standage
2017-05-23 00:13:50 UTC
view on stackexchange narkive permalink

Het wgsim-pakket van Heng Li (bekend van BWA en samtools) is mijn go-to-tool voor het simuleren van Illumina-reads. Het biedt geen gemakkelijke manier om differentiële dekking over verschillende reeksen te simuleren, maar het zou niet zo moeilijk moeten zijn om wgsim meerdere keren uit te voeren, waardoor het gewenste dekkingsniveau voor elke van belang zijnde reeks wordt gegenereerd.

I zou een Python-script implementeren om je testbestand op te slurpen, en wgsim aanroepen (met behulp van de subproces -module) voor elke reeks. Dit vereist waarschijnlijk dat u elke reeks in een apart bestand moet hebben. :-(

Kunt u dit uitbreiden door een voorbeeldopdracht toe te voegen die laat zien hoe u een invoer fasta-bestand neemt en de fastq produceert? Dit voelt meer als een (nuttige) routebeschrijving naar waar een antwoord te vinden is dan als een antwoord.
#5
+2
Gabriel Renaud
2017-05-17 21:36:37 UTC
view on stackexchange narkive permalink

Ik ken geen software die dit rechtstreeks kan doen, maar ik zou het fasta-bestand in één reeks per bestand splitsen, ze in BASH doorlopen en ART de sequentiesimulator (of een andere) op elke reeks oproepen.

Kunt u hier iets over zeggen? Dit is in deze toestand echt geen erg nuttig antwoord. Wat is KUNST? Waar kan ik het vinden? Hoe gebruik ik het? Hoe werkt het? Introduceert het fouten? Kunnen we het foutenpercentage controleren? Zijn er ook echt geen multifasta-bestanden nodig? Als dit niet het geval is, zal het herhalen van enkele duizenden bestanden in bash erg, erg traag zijn. Er zijn misschien betere manieren.
> Wat is KUNST? Het is moeilijk te definiëren, ik zou zeggen dat het verschillende activiteiten zijn die betrekking hebben op het produceren van werken van visuele of auditieve producten die gericht zijn op het externaliseren van het innerlijke van de auteur ... oh, bedoel je de sequentiesimulator? Zie hun paper hier: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3278762/> Waar kan ik het vinden? zie papier> Hoe gebruik ik het? lees handleiding> Hoe werkt het? lees papier.> Introduceert het fouten? lees paper / manual Kunnen we het foutenpercentage controleren? lees papier.
> Zijn er ook echt geen multifasta-bestanden nodig? Als dit niet het geval is, zal het herhalen van enkele duizenden bestanden in bash erg, erg traag zijn. Er kunnen betere manieren zijn Waar heeft iemand gezegd dat het geen multifasta gebruikt? OP wilde een zeer specifieke dekking per fasta-record.
Nee, ik bedoelde het niet aan mij uit te leggen. Ik bedoel, bewerk je antwoord, zodat het deze informatie bevat. We proberen antwoorden zo volledig mogelijk te hebben en zo geschreven dat ze op zichzelf kunnen staan ​​zonder externe referenties op de Stack Exchange-sites. Opmerkingen zijn moeilijk te lezen, gemakkelijk te missen en resulteren in de vreselijke puinhoop van biostars waarvan je geen idee hebt waar het daadwerkelijke antwoord is. Bewerk uw vraag om deze te verduidelijken.
Uw antwoord, bedoelde ik. Sorry.
#6
+2
Greg
2017-05-17 22:55:25 UTC
view on stackexchange narkive permalink

Als je het goed vindt met wat willekeur, kun je reads genereren uit je reeksbestand met behulp van een willekeurige Poisson-variabele. Je zult wat wiskunde moeten doen om erachter te komen welke waarde van lambda je moet gebruiken om ervoor te zorgen dat de verwachte dekking bij elk basenpaar in je read overeenkomt met wat je in je tekstbestand hebt ingesteld.

Bijvoorbeeld jij hebben een reeks S met een lengte van 1.000, een leeslengte van 50 en een invoeggrootte van 100. Genereer voor elke basis b in S een willekeurige Poisson-variabele p. U genereert dan p leest van basis b tot b + 50. Genereer vervolgens de gepaarde uitlezing beginnend bij b + 50 + 100.

Nogmaals, je zou ermee moeten spelen om erachter te komen welke lambda je moet gebruiken, maar dit zou je eigenlijk geven wat je wilt, zolang je het goed vindt dat je niet precies de dekking hebt die je target voor elke gelezen.

En zal ook realistischer zijn om de wisselplaatafstand van Poisson-achtige distributie te simuleren! Het punt is dat dit een heruitvinding is van een fiets. Er zijn veel simulatietools.
Zou trekken vanuit een normale verdeling niet beter zijn voor de insteekmaat? En ja, ik weet dat er veel tools zijn, maar dit zou een manier zijn om het zelf te doen als je geïnteresseerd was.
uit mijn ervaring is normaal iets minder flexibel voor aanpassing - de dichtheid is scheef, bijv. http://bioinformatics.ucdavis.edu/docs/2016-june-workshop/_images/picard-insertion-size-histogram.png. Dus waarschijnlijk is een negatieve binominale keuze een goede keuze, maar ik heb nooit geprobeerd deze aan te passen aan de werkelijke dichtheden van de wisselplaatafmeting ...
#7
+1
Karel Brinda
2017-06-21 19:46:49 UTC
view on stackexchange narkive permalink

Het simuleren van NGS-uitlezingen terwijl de reeksdekking wordt beheerd, is nu eenvoudig met RNFtools (vanaf versie 0.3.1). Zie de tutorial, met name de sectie Sequentie-extractie.

Omgevingsvoorbereiding

Installeer eerst BioConda en voeg de vereiste kanalen toe. Installeer vervolgens RNFtools in de standaard Conda-omgeving

  conda install rnftools  

of maak en activeer een aparte Conda omgeving (bij voorkeur)

  conda create -n rnftools rnftoolssource activeren rnftools  

Simulatie

Stel dat u een referentiebestand ref.fa en een door tabs gescheiden dekkingsbestand coverage.tsv heeft (bijv. die uit uw voorbeeld) . Dan zullen de volgende RNFtools Snakefile het werk doen dat je wilt:

  import rnftoolsimport csvrnftools.mishmash.sample ("simulation_with_coverage_control" , reads_in_tuple = 1) fa = "ref.fa" tsv = "coverage.tsv" met open (tsv) als f: table = csv.reader (f, delimiter = '\ t') voor seqname, cov in tabel: rnftools .mishmash.DwgSim (fasta = fa, sequences = [seqname], coverage = float (cov), read_length_1 = 10, # snelle test met supersnelle leest read_length_2 = 0,) include: rnftools.include () rule: input: rnftools. input ()  

Als u dit bestand opslaat ( Snakefile ) en snakemake uitvoert, zal RNFtools reads simuleren met DWGsim met de gedefinieerde dekkingen je tekstbestand, en sla alle gesimuleerde reads op in simulation_with_coverage_control.fq.

Je kunt met alle parameters spelen. In het bijzonder kunt u een andere simulator gebruiken (bijv. Art-Illumina met rnftools.mishmash.ArtIllumina ). Zie de documentatie van RNFtools voor meer informatie.

#8
  0
winni2k
2017-06-12 01:27:23 UTC
view on stackexchange narkive permalink

Een ander leessimulatietool van de volgende generatie is gemsim. Ik heb het niet getest, maar ik zou geïnteresseerd zijn als iemand er ervaring mee heeft.

#9
-1
Karel Brinda
2017-05-30 00:46:25 UTC
view on stackexchange narkive permalink

U kunt uw FASTA-bestand volgorde-voor-volgorde splitsen met

  split -a 6 -p '^ >' uw_bestand.fa seq_  

en gebruik vervolgens een bestaande leessimulator die dekking ondersteunt (ART, DWGsim, enz.). Als je alle reads gemengd wilt hebben (niet geordend op de originele volgorde), kun je RNFtools gebruiken.

Bewerken 1:

Zoals @terdon erop gewezen, werkt het vorige commando alleen op OS X. Een analoge one-liner voor Linux (maar met een iets ander naamgevingsschema met cijfers in plaats van letters) kan zijn

  csplit -f seq_ -n 6 your_file.fa '/ ^ > /' {* }  

Om dit commando ook op OS X te laten werken, moet je coreutils installeren (bijvoorbeeld door brew te gebruiken) en vervolgens gcsplit gebruiken in plaats van csplit .

Bewerken 2:

Als FASTA eenmaal is opgesplitst in reeksen, wordt de simulatie eenvoudig en kunnen veel verschillende benaderingen worden gebruikt. Mijn favoriete is het gebruik van GNU Parallel. Stel je voor dat je dekkingen in een tekstbestand met de naam covs.txt op aparte regels staan ​​en in dezelfde volgorde als de reeksen in jouw_bestand.fa , bijvoorbeeld

4030...

Vervolgens kun je leesbewerkingen simuleren van de originele sequenties met DWGsim door

  ls -1 seq_ * | plak covs.txt - \ | parallel -v --colsep '\ t' dwgsim -1 100 -2 100 -C {1} {2} sim_ {2}  

en voeg de verkregen FASTQ-bestanden samen met:

  cat sim_seq _ *. bwa.read1.fastq > reads.1.fqcat sim_seq _ *. bwa.read2.fastq > reads.2.fq  

Een mogelijk Het gevaar van deze benadering is dat we hebben aangenomen dat het aantal seq_ * bestanden hetzelfde is als het aantal regels in covs.txt , wat misschien niet waar is (per ongeluk). We moeten dit voorafgaand aan de simulatiestap controleren, bijvoorbeeld door:

  [["$ (ls -1 seq_ * | wc -l)" == "$ (cat covs.txt | wc -l) "]] \ || echo "Incorrent aantal regels in covs.txt"  

Een ander voorbehoud is dat de gesimuleerde reads niet in willekeurige volgorde staan ​​(ze zijn gegroepeerd op bronsequenties).

Ik ben bang dat u de vraag niet beantwoordt. Allereerst gebruikt u niet-standaard 'sorteer'-opties. Ik neem aan dat je een Mac gebruikt? De `-p` is een alleen BSD-optie en niet beschikbaar op andere * nix-systemen AFAIK. Dat gezegd hebbende, hoewel dit waarschijnlijk zal werken om het bestand op een MacOS-systeem te splitsen, is het splitsen ervan triviaal en er zijn veel, vele manieren om dit te doen. Het moeilijke is het simuleren van de leesbewerkingen en je legt niet echt uit hoe het OP dat kan doen.
Bedankt voor de bewerking. Het belangrijkste probleem hier is echter dat uw antwoord geen antwoord is op de gestelde vraag. U antwoordt "hoe kan ik een multifasta-bestand opsplitsen in veel bestanden met een enkele reeks", maar de vraag gaat over het simuleren van leesbewerkingen. Het splitsen van het bestand kan * deel uitmaken van * een antwoord, of misschien niet, maar het is zeker geen antwoord.
Oh, en voor een draagbare manier zonder een van de vele gespecialiseerde tools hiervoor te gebruiken, kun je `awk` gebruiken. Iets als `awk -vRS ="> "-F '\ n' 'NR> 1 {print"> "$ 0 >> $ 1" .fa "}' file.fa`. Uw `csplit` (en uw` split`) benadering is echter inderdaad eenvoudiger.


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