Pitanje:
Kako simulirati čitanje NGS-a, kontrolirajući pokrivenost sekvence?
SmallChess
2017-05-17 21:26:17 UTC
view on stackexchange narkive permalink

Imam FASTA datoteku sa više od 100 sekvenci poput ove:

  >Sequence1GTGCCTATTGCTACTAAAA ... >Sequence2GCAATGCAAGGAAGTGATGGCGGAAATAGCGTpre ...... > > > tekstualna datoteka poput ove:  
  Sequence1 40Sequence2 30 ......  

Želio bih simulirati čitanje uparenih krajeva sljedeće generacije za sve sekvence u mojoj FASTA datoteci. Za Sequence1 , želio bih simulirati pokrivenost od 40x. Za Sequence2 , želio bih simulirati pokrivenost 30x. Drugim riječima, želim kontrolirati pokrivenost sekvence za svaku sekvencu u svojoj simulaciji.

P: Koji je najjednostavniji način za to? Bilo koji softver koji bih trebao koristiti? Bioprovodnik?

[Ova recenzija] (http://www.nature.com/nrg/journal/v17/n8/abs/nrg.2016.57.html) može biti dobro polazište. Uspoređuje 23 NGS alata za simulaciju i nudi stablo odluke za odabir odgovarajućeg alata na temelju vaših potreba.
Koju dužinu čitanja koristite? Koliko su duge sekvence? Morate li točno ili s izvjesnom vjerojatnošću pogoditi cilj pokrivenosti?
Zvuči kao nešto što biste samo kodirali u deset ili dvadeset redaka pythona pomoću modula Biopython.
Gregovim bih dodao još nekoliko pitanja. Razumijem li dobro da želite simulirati redoslijed niza predloška iz datoteke? Dakle, predstavljaju li sekvence genome? Želite li uzeti u obzir pristranost pojačanja? Pogreške u sekvenciranju? Koju platformu za sekvenciranje želite simulirati?
Devet odgovori:
#1
+6
H. Gourlé
2017-05-30 12:06:03 UTC
view on stackexchange narkive permalink

Radim na simulatoru sekvenciranja Illumina za metagenomiku: InSilicoSeq

Još uvijek je u alfa izdanju i vrlo je eksperimentalan, ali s obzirom na datoteku s više datoteka i obilje generirat će očitanja iz vaših ulaznih genoma s različitim pokrivanjima.

Iz dokumentacije:

  iss generirati --genomes genomes.fasta --obilje obilja_datoteka.txt \ - model_file HiSeq2500 --izlaz HiSeq_reads  

Gdje:

  # multi-fasta file>genome_AATGC ... >genome_BCCGT ...... # datoteka obilja (ukupno obilje mora biti 1!) genome_A 0.2genome_B 0.4 ...  

Nisam ga dizajnirao da radi s pokrivanjem, već s obiljem genoma u metagenomu, pa ćete možda morati napraviti malo matematike;)

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

Paket bioloških vodiča poliester može to učiniti. Kaže da simulira čitanje RNA-seq, ali ne znam da li se to doista razlikuje od ostalih NGS čitanja.

Može koristiti niz pogrešaka i modela pristranosti ili ih naučiti iz skupa podataka.

Možete li malo proširiti ovo? Kako bi ga OP koristio na svojim datotekama? Koji opseg pogrešaka i pristranosti? Kako bismo ih koristili? Što bi bilo razumno zadati za upotrebu? Možda bi vam bio od pomoći minimalni radni primjer.
#3
+4
Kamil S Jaron
2017-05-18 22:32:11 UTC
view on stackexchange narkive permalink

Ova python skripta uzima fasta datoteku i tsv datoteku s brojevima i ispisuje sekvence u fasta datotekama onoliko puta koliko je navedeno u tsv datoteci (pod pretpostavkom formata u pitanju). Pa ako su bar.tsv i foo.fasta vaše datoteke:

  iz Bio import SeqIOrepeat = {} za otvoreni redak ( "bar.tsv"): seq_id, pokrivenost = line.split () repeat [seq_id] = int (pokrivenost) za seq_record u SeqIO.parse (foo.fasta, "fasta"): za i u dometu (repeat.get ( seq_record.name, 0)): print (">", seq_record.name, "_", i, sep = '') print (seq_record.seq)  
Možete li objasniti što to radi i kako to radi? Čini se da ćete samo ispisati originalnu sekvencu $ vremena pokrivanja, a to ne može biti točno. OP mora simulirati čitanje, a ne samo ponoviti isti slijed N puta.
Ajme, onda sam krivo razumio pitanje. Iako su sekvence čitanja koja treba simulirati, a 40x je bočno, volio bih imati ovaj niz 40 puta. Uredio sam odgovor kako bih barem izbjegao zabunu.
Ne znam jeste li pogrešno shvatili ili ja jesam. Ali jedan od nas jest :) Hvala na uređivanju, barem sada OP može vidjeti što to čini i odlučiti sami.
#4
+4
Daniel Standage
2017-05-23 00:13:50 UTC
view on stackexchange narkive permalink

Paket wgsim autora Henga Lija (od slave BWA i samtools) moj je alat za simulaciju čitanja Ilumine. Ne pruža prikladan način za simulaciju diferencijalne pokrivenosti kroz različite sekvence, ali ne bi trebalo biti teško pokretati wgsim više puta, generirajući željenu razinu pokrivenosti za svaki slijed od interesa.

I bi implementirao Python skriptu da izlupa vašu testnu datoteku i pozvao wgsim (koristeći modul potproces ) za svaku sekvencu. To će vjerojatno zahtijevati da svaki slijed imate u zasebnoj datoteci. :-(

Možete li to proširiti dodavanjem primjera naredbe koja pokazuje kako uzeti ulaznu fasta datoteku i proizvesti fastq? Ovo više sliči na (korisne) upute do mjesta gdje se može pronaći odgovor nego na odgovor.
#5
+2
Gabriel Renaud
2017-05-17 21:36:37 UTC
view on stackexchange narkive permalink

Nisam upoznat s bilo kojim softverom koji to može izravno učiniti, ali podijelio bih datoteku fasta u jednu sekvencu po datoteci, prebacio ih preko BASH-a i pozvao ART simulator sekvence (ili drugi) na svakoj sekvenci.

Možete li proširiti ovo? Ovo zaista nije vrlo koristan odgovor u ovom stanju. Što je ART? Gdje ga mogu naći? Kako da ga koristim? Kako iot djeluje? Uvodi li pogreške? Možemo li kontrolirati stopu pogrešaka? Također, nisu li potrebne datoteke multifasta? Ako ne, pregledavanje nekoliko tisuća datoteka u bashu bit će vrlo, vrlo sporo. Možda postoje bolji načini.
> Što je ART? Teško je definirati, rekao bih da se radi o raznim aktivnostima koje se odnose na produkciju djela vizualnih ili slušnih proizvoda usmjerenih na eksternalizaciju autorove unutarnje ... o, mislite na simulator sekvence? Molimo pogledajte njihov članak ovdje: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3278762/> Gdje ga mogu pronaći? pogledajte papir> Kako ga koristiti? pročitajte priručnik> Kako iot djeluje? pročitati rad.> Uvodi li pogreške? pročitati rad / priručnik Možemo li kontrolirati stopu pogrešaka? čitati papir.
> Također, zar stvarno ne trebaju multifasta datoteke? Ako ne, pregledavanje nekoliko tisuća datoteka u bashu bit će vrlo, vrlo sporo. Možda postoje bolji načiniGdje je netko izjavio da ne treba multifasta? OP je želio vrlo specifičnu pokrivenost po fasta zapisu.
Ne, nisam to mislio objasniti. Mislio sam, molim vas uredite svoj odgovor tako da sadrži ove podatke. Nastojimo da odgovori budu što opsežniji i napisani kako bi mogli samostalno stajati bez vanjskih referenci na web mjestima Stack Exchange. Komentari su teški za čitanje, lako ih je propustiti i rezultiraju užasnim neredom biostarsa ​​gdje nemate pojma gdje je stvarni odgovor. Uredite svoje pitanje da biste pojasnili.
Mislio sam na tvoj odgovor. Oprosti.
#6
+2
Greg
2017-05-17 22:55:25 UTC
view on stackexchange narkive permalink

Ako se slažete s nekom slučajnošću, možete generirati čitanja iz datoteke sekvence pomoću Poissonove slučajne varijable. Morat ćete malo izračunati kako biste shvatili koju vrijednost lambde trebate koristiti kako bi se očekivano pokrivanje svakog osnovnog para u vašem čitanju podudaralo s onim što ste postavili u svojoj tekstualnoj datoteci.

Na primjer, vi imaju niz S dužine 1.000, dužine čitanja 50 i veličine umetka 100. Za svaku bazu b u S generirajte Poissonovu slučajnu varijablu p. Zatim ćete generirati p čitanja od baze b do b + 50. Zatim generirajte upareno čitanje počevši od b + 50 + 100.

Opet, morali biste se igrati s njom da biste shvatili koju lambdu koristiti, ali to bi vam u osnovi dalo ono što želite, sve dok se slažete s time što nemate točno pokrivenost koju ciljate za svako čitanje.

A također će biti realnije simulirati udaljenost umetka od Poissonove distribucije! Stvar je u tome što je ovo ponovno pronalazak bicikla. Postoji puno alata za simulaciju.
Ne bi li izvlačenje iz normalne distribucije bilo bolje za veličinu umetanja? I da, znam da postoji puno alata, ali ovo bi bio jedan od načina da to sami napravite da vas zanima.
iz mog iskustva, normalno je malo manje fleksibilno za ugradnju - gustoća je iskrivljena, npr. http://bioinformatics.ucdavis.edu/docs/2016-june-workshop/_images/picard-insertion-size-histogram.png. Dakle, vjerojatno će Negativni binom biti dobar izbor, ali nikada ga nisam pokušao prilagoditi stvarnoj gustoći veličine umetka ...
#7
+1
Karel Brinda
2017-06-21 19:46:49 UTC
view on stackexchange narkive permalink

Simulacija čitanja NGS-a dok se kontrolira pokrivenost sekvenci sada je lako s RNFtools (od verzije 0.3.1). Pogledajte vodič, posebno odjeljak Izvlačenje sekvence.

Priprema okoline

Prvo instalirajte BioConda i dodajte potrebne kanale. Zatim instalirajte RNFtools u zadano Conda okruženje

  conda install rnftools  

ili stvorite i aktivirajte zasebni Conda okruženje (poželjno)

  conda create -n rnftools rnftoolssource aktiviranje rnftools  

Simulacija

Pretpostavimo da imate referentnu datoteku ref.fa i datoteku pokrivenosti odvojenu karticama pokrivanje.tsv (npr. one iz vašeg primjera) . Tada će sljedeći RNFtools Snakefile obaviti posao koji želite:

  import rnftoolsimport csvrnftools.mishmash.sample ("simulation_with_coverage_control" , read_in_tuple = 1) fa = "ref.fa" tsv = "охват.tsv" s otvorenim (tsv) kao f: table = csv.reader (f, graničnik = '\ t') za seqname, cov u tablici: rnftools .mishmash.DwgSim (fasta = fa, sekvence = [seqname], pokrivenost = float (cov), read_length_1 = 10, # brzi test s superkratkom čitanjem read_length_2 = 0,) uključuje: pravilo rnftools.include (): input: rnftools. input ()  

Kada spremite ovu datoteku ( Snakefile ) i pokrenete snakemake , RNFtools će simulirati čitanje pomoću DWGsim s definiranim pokrićima svoju tekstualnu datoteku i spremite sva simulirana čitanja u simulation_with_coverage_control.fq.

Možete se igrati sa svim parametrima. Konkretno, možete koristiti drugi simulator (npr. Art-Illumina pomoću rnftools.mishmash.ArtIllumina ). Pogledajte RNFtools dokumentaciju za više informacija.

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

Drugi alat za simulaciju čitanja sljedeće generacije je gemsim. Nisam ga testirao, ali zanimalo bi me je li netko imao iskustva s tim.

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

Možete podijeliti FASTA datoteku redom po slijedu pomoću

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

a zatim upotrijebite bilo koji postojeći simulator čitanja koji podržava pokrivenost (ART, DWGsim, itd.). Ako želite da se sva čitanja izmiješaju (nisu poredana po izvornom slijedu), možete koristiti RNFtools.

Uredi 1:

Kao @terdon istaknuto je da prethodna naredba radi samo na OS X. Analogna linija za Linux (ali s malo drugačijom shemom imenovanja koja koristi brojeve, a ne slova) može biti

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

Da bi ova naredba djelovala i na OS X, treba instalirati coreutils (npr. pomoću brew-a), a zatim koristiti gcsplit umjesto csplit .

Uredi 2:

Jednom kada se FASTA podijeli nizovima, simulacija postaje jednostavna i mogu se koristiti mnogi različiti pristupi. Najdraži mi je koristiti GNU Parallel. Zamislite da imate pokrivenosti u tekstualnoj datoteci koja se naziva covs.txt u zasebnim redovima i istim redoslijedom kao nizovi u your_file.fa , npr.

4030...

Tada možete simulirati čitanje iz izvornih sekvenci pomoću DWGsim

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

i spojite dobivene FASTQ datoteke koristeći:

  mačka sim_seq _ *. bwa.read1.fastq > čita.1.fqcat sim_seq _ *. bwa.read2.fastq > čita.2.fq  

Jedna moguća Opasnost ovog pristupa je u tome što smo pretpostavili da je broj datoteka seq_ * jednak broju redaka u covs.txt , što možda nije istina (pogreškom). To bismo trebali provjeriti prije koraka simulacije, npr.,

  [["$ (ls -1 seq_ * | wc -l)" == "$ (cat covs.txt | wc -l) "]] \ || echo "Neispravan broj redaka u covs.txt"  

Još je jedna zamjerka da simulirana čitanja nisu slučajnim redoslijedom (grupirana su prema izvornim sekvencama).

Bojim se da ne odgovarate na pitanje. Prije svega, koristite nestandardne opcije `sortiranja`. Pretpostavljam da koristite Mac? `-P` je opcija samo za BSD i nije dostupna na ostalim * nix sustavima AFAIK. To će reći, iako će ovo vjerojatno uspjeti podijeliti datoteku na MacOS sustavu, ali njezino je dijeljenje trivijalno i postoji mnogo, mnogo načina za to. Tvrdi bit simulira čitanja i zapravo ne objašnjavate kako OP to može učiniti.
Hvala na uređivanju. Međutim, glavno pitanje ovdje je da vaš odgovor nije odgovor na postavljeno pitanje. Odgovorite "kako mogu podijeliti multifasta datoteku u mnogo datoteka s jednim slijedom", ali pitanje je o simulaciji čitanja. Dijeljenje datoteke možda je * dio * odgovora ili ne, ali sigurno nije odgovor.
Oh, i za prijenosni način bez upotrebe bilo kojeg od mnogih specijaliziranih alata za to možete koristiti `awk`. Nešto poput `awk -vRS ="> "-F '\ n' 'NR> 1 {print"> "$ 0 >> $ 1" .fa "}' file.fa`. Vaš je `csplit` (i vaš` split`) pristup uistinu jednostavniji.


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