Domanda:
Come simulare le letture NGS, controllando la copertura della sequenza?
SmallChess
2017-05-17 21:26:17 UTC
view on stackexchange narkive permalink

Ho un file FASTA con oltre 100 sequenze come questa:

  >Sequence1GTGCCTATTGCTACTAAAA ... >Sequence2GCAATGCAAGGAAGTGATGGCGGAAATAGCGTTA ......  
un file di testo come questo:
  Sequence1 40Sequence2 30 ......  

Vorrei simulare letture paired-end di prossima generazione per tutti le sequenze nel mio file FASTA. Per Sequence1 , vorrei simulare con una copertura 40x. Per Sequence2 , vorrei simulare con una copertura 30x. In altre parole, voglio controllare la copertura della mia sequenza per ogni sequenza nella mia simulazione.

D: Qual è il modo più semplice per farlo? Qualche software dovrei usare? Bioconduttore?

[Questa recensione] (http://www.nature.com/nrg/journal/v17/n8/abs/nrg.2016.57.html) potrebbe essere un buon punto di partenza. Confronta 23 strumenti di simulazione NGS e fornisce un albero decisionale per la selezione di uno strumento appropriato in base alle proprie esigenze.
Qual è la lunghezza di lettura che stai utilizzando? Quanto durano le sequenze? Hai bisogno di raggiungere l'obiettivo di copertura esattamente o con una certa probabilità?
Suona come qualcosa che dovresti semplicemente codificare in dieci o venti righe di Python con il modulo Biopython.
Vorrei aggiungere qualche altra domanda a quella di Greg. Ho capito bene che vorresti simulare la sequenza della sequenza del modello dal file? Quindi le sequenze rappresentano i genomi? Vorresti considerare il bias di amplificazione? Errori di sequenza? Quale piattaforma di sequenziamento vorresti simulare?
Nove risposte:
#1
+6
H. Gourlé
2017-05-30 12:06:03 UTC
view on stackexchange narkive permalink

Sto lavorando a un simulatore di sequenziamento Illumina per la metagenomica: InSilicoSeq

È ancora in versione alpha e molto sperimentale, ma dato un file multi-fasta e abbondanza, genererà letture dai tuoi genomi di input con diverse coperture.

Dalla documentazione:

  iss generate --genomes genomes.fasta --abundance abbondance_file.txt \ - model_file HiSeq2500 --output HiSeq_reads  

Dove:

  # multi-fasta file>genome_AATGC ... >genome_BCCGT ...... # file di abbondanza (abbondanza totale deve essere 1!) genome_A 0.2genome_B 0.4 ...  

Non l'ho progettato per funzionare con la copertura ma piuttosto l'abbondanza del genoma in un metagenoma, quindi potresti dover fare un un po 'di matematica;)

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

Il pacchetto bioconduttore poliestere può farlo. Dice che simula le letture RNA-seq, ma non so se sia davvero diverso dalle altre letture NGS.

Può utilizzare una serie di modelli di errore e bias o apprenderli da un set di dati.

Potresti approfondire questo un po '? Come lo userebbe l'OP sui propri file? Quale gamma di errori e pregiudizi? Come li useremmo? Quale sarebbe un valore predefinito ragionevole da utilizzare? Forse sarebbe utile un minimo esempio di lavoro.
#3
+4
Kamil S Jaron
2017-05-18 22:32:11 UTC
view on stackexchange narkive permalink

Questo script python prende un file fasta e un file tsv con i conteggi e stampa le sequenze nei file fasta quante volte è specificato nel file tsv (assumendo il formato nella domanda). Quindi, se bar.tsv e foo.fasta saranno i tuoi file:

  from Bio import SeqIOrepeat = {} for line in open ( "bar.tsv"): seq_id, coverage = line.split () repeat [seq_id] = int (coverage) for seq_record in SeqIO.parse (foo.fasta, "fasta"): for i in range (repeat.get ( seq_record.name, 0)): print (">", seq_record.name, "_", i, sep = '') print (seq_record.seq)  
Potresti spiegare cosa fa e come funziona? Sembra che stamperai solo la sequenza originale $ coverage times, e questo non può essere corretto. L'OP deve simulare le letture, non solo ripetere la stessa sequenza N volte.
Ahi, allora ho frainteso la domanda. Anche se le sequenze sono le letture che dovrebbero essere simulate e 40x è lateralmente, vorrei avere questa sequenza 40 volte. Ho modificato la risposta almeno per evitare confusione.
Non so se hai capito male, o l'ho fatto. Ma uno di noi l'ha fatto :) Grazie per la modifica, almeno ora l'OP può vedere cosa fa e decidere da solo.
#4
+4
Daniel Standage
2017-05-23 00:13:50 UTC
view on stackexchange narkive permalink

Il pacchetto wgsim di Heng Li (di BWA e fama di samtools) è il mio strumento di riferimento per simulare le letture Illumina. Non fornisce alcun modo conveniente per simulare la copertura differenziale tra sequenze diverse, ma non dovrebbe essere troppo difficile eseguire wgsim più volte, generando il livello di copertura desiderato per ogni sequenza di interesse.

I implementerebbe uno script Python per consumare il file di test e richiamerebbe wgsim (usando il modulo subprocess ) per ogni sequenza. Questo probabilmente richiederà di avere ogni sequenza in un file separato. :-(

Potresti espandere questo aggiungendo un comando di esempio che mostra come prendere un file fasta di input e produrre il fastq? Sembra più un'indicazione (utile) per trovare una risposta che una risposta.
#5
+2
Gabriel Renaud
2017-05-17 21:36:37 UTC
view on stackexchange narkive permalink

Non sono a conoscenza di alcun software in grado di farlo direttamente, ma dividerei il file fasta in una sequenza per file, eseguirò un loop su di essi in BASH e richiamerei ART il simulatore di sequenza (o un altro) su ciascuna sequenza.

Potresti approfondire questo? Questa non è davvero una risposta molto utile in questo stato. Cos'è l'ARTE? Dove posso trovarlo? Come lo uso? Come funziona lo iot? Presenta errori? Possiamo controllare il tasso di errore? Inoltre, davvero non accetta file multifasta? In caso contrario, il loop di diverse migliaia di file in bash sarà molto, molto lento. Potrebbero esserci modi migliori.
> Cos'è l'ARTE? Difficile da definire, direi che si tratta di attività varie relative alla produzione di lavori di prodotti visivi o uditivi volti ad esternare l'interiorità dell'autore ... oh, intendi il simulatore di sequenze? Si prega di consultare il loro documento qui: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3278762/> Dove posso trovarlo? vedi carta> Come si usa? leggi il manuale> Come funziona lo iot? leggere la carta.> Presenta errori? leggere carta / manuale Possiamo controllare il tasso di errore? leggere la carta.
> Inoltre, davvero non accetta file multifasta? In caso contrario, il loop di diverse migliaia di file in bash sarà molto, molto lento. Potrebbero esserci modi migliori Dove qualcuno ha affermato che non richiede multifasta? OP voleva una copertura molto specifica per il record di fasta.
No, non volevo spiegarmelo. Intendevo modificare la tua risposta in modo che fornisca queste informazioni. Cerchiamo di avere risposte il più complete possibile e scritte in modo che possano stare da sole senza riferimenti esterni sui siti Stack Exchange. I commenti sono difficili da leggere, facili da perdere e si traducono nell'orribile pasticcio di biostars in cui non hai idea di dove sia la risposta effettiva. Modifica la tua domanda per chiarire.
La tua risposta, volevo dire. Scusate.
#6
+2
Greg
2017-05-17 22:55:25 UTC
view on stackexchange narkive permalink

Se sei d'accordo con un po 'di casualità, puoi generare letture dal tuo file di sequenza utilizzando una variabile casuale di Poisson. Dovrai fare un po 'di matematica per capire quale valore di lambda usare affinché la copertura prevista per ogni coppia di basi nella tua lettura corrisponda a ciò che hai impostato nel tuo file di testo.

Per esempio tu hanno una sequenza S di lunghezza 1.000, una lunghezza di lettura di 50 e una dimensione di inserimento di 100. Per ogni base b in S generare una variabile casuale di Poisson p. Genererai quindi letture p dalla base b a b + 50. Quindi, genera la lettura accoppiata a partire da b + 50 + 100.

Ancora una volta, dovresti giocarci per capire quale lambda usare, ma questo ti darebbe fondamentalmente quello che vuoi, a patto che tu stia bene con non avere esattamente la copertura che stai mirando per ogni lettura.

Inoltre, sarà più realistico simulare la distanza di inserimento da una distribuzione simile a quella di Poisson! Il fatto è che questa è una reinvenzione di una bicicletta. Ci sono molti strumenti di simulazione.
Non sarebbe meglio tirare da una distribuzione normale per la dimensione dell'inserimento? E sì, so che ci sono molti strumenti, ma questo sarebbe un modo per farlo da solo se fossi interessato.
dalla mia esperienza, normale è un po 'meno flessibile per l'adattamento - la densità è distorta, ad esempio http://bioinformatics.ucdavis.edu/docs/2016-june-workshop/_images/picard-insertion-size-histogram.png. Quindi probabilmente un binomio negativo sarà una buona scelta, ma non ho mai provato ad adattarlo alle densità delle dimensioni effettive dell'inserto ...
#7
+1
Karel Brinda
2017-06-21 19:46:49 UTC
view on stackexchange narkive permalink

Simulare le letture NGS mentre si controlla la copertura della sequenza ora è facile con RNFtools (dalla versione 0.3.1). Vedi il tutorial, in particolare la sezione Estrazione di sequenze.

Preparazione dell'ambiente

Innanzitutto, installa BioConda e aggiungi i canali richiesti. Quindi installa RNFtools nell'ambiente Conda predefinito

  conda install rnftools  

o crea e attiva una Conda separata ambiente (preferibile)

  conda create -n rnftools rnftoolssource activation rnftools  

Simulazione

Supponi di avere un file di riferimento ref.fa e un file di copertura separato da tabulazioni coverage.tsv (ad esempio, quelli del tuo esempio) . Quindi il seguente RNFtools Snakefile farà il lavoro che desideri:

  import rnftoolsimport csvrnftools.mishmash.sample ("simulation_with_coverage_control" , reads_in_tuple = 1) fa = "ref.fa" tsv = "coverage.tsv" with open (tsv) as f: table = csv.reader (f, delimiter = '\ t') for seqname, cov in table: rnftools .mishmash.DwgSim (fasta = fa, sequences = [seqname], coverage = float (cov), read_length_1 = 10, # quick test with supershort reads read_length_2 = 0,) include: rnftools.include () regola: input: rnftools. input ()  

Quando salvi questo file ( Snakefile ) ed esegui snakemake , RNFtools simulerà le letture usando DWGsim con le coperture definite il tuo file di testo e salva tutte le letture simulate in simulation_with_coverage_control.fq.

Puoi giocare con tutti i parametri. In particolare, è possibile utilizzare un simulatore diverso (ad esempio, Art-Illumina utilizzando rnftools.mishmash.ArtIllumina ). Consulta la documentazione di RNFtools per ulteriori informazioni.

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

Un altro strumento di simulazione della lettura di nuova generazione è gemsim. Non l'ho testato, ma sarei interessato se qualcuno avesse avuto esperienza con esso.

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

Puoi dividere il tuo file FASTA sequenza per sequenza usando

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

e quindi utilizzare qualsiasi simulatore di lettura esistente che supporti la copertura (ART, DWGsim, ecc.). Se vuoi che tutte le letture siano miste (non ordinate dalla sequenza originale), puoi usare RNFtools.

Modifica 1:

Come @terdon sottolineato, il comando precedente funziona solo su OS X. Un one liner analogico per Linux (ma con uno schema di denominazione leggermente diverso che utilizza numeri anziché lettere) può essere

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

Per far funzionare questo comando anche su OS X, è necessario installare coreutils (ad esempio, utilizzando brew) e quindi utilizzare gcsplit invece di csplit .

Modifica 2:

Una volta che FASTA è diviso per sequenze, la simulazione diventa semplice e possono essere utilizzati molti approcci diversi. Il mio preferito è usare GNU Parallel. Immagina di avere le tue copertine in un file di testo chiamato covs.txt su righe separate e nello stesso ordine delle sequenze in your_file.fa , ad esempio

4030...

Quindi puoi simulare letture dalle sequenze originali usando DWGsim

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

e unisci i file FASTQ ottenuti usando:

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

Uno possibile Il pericolo di questo approccio è che abbiamo assunto che il numero di file seq_ * sia uguale al numero di righe in covs.txt , il che potrebbe non essere vero (per errore). Dovremmo verificarlo prima della fase di simulazione, ad esempio:

  [["$ (ls -1 seq_ * | wc -l)" == "$ (cat covs.txt | wc -l) "]] \ || echo "Numero errato di righe in covs.txt"  

Un altro avvertimento è che le letture simulate non sono in ordine casuale (sono raggruppate per sequenze di origine).

Temo che tu non risponda alla domanda. Prima di tutto, stai usando opzioni di "ordinamento" non standard. Presumo che usi un Mac? `-P` è un'opzione solo BSD e non è disponibile su altri sistemi * nix AFAIK. Detto questo, mentre questo probabilmente funzionerà per dividere il file su un sistema MacOS, dividerlo è banale e ci sono molti, molti modi per farlo. La parte difficile è simulare le letture e in realtà non stai spiegando come l'OP può farlo.
Grazie per la modifica. Tuttavia, il problema principale qui è che la tua risposta non risponde alla domanda posta. Stai rispondendo "come posso dividere un file multifasta in molti file a sequenza singola" ma la domanda riguarda la simulazione delle letture. La divisione del file potrebbe essere * parte * di una risposta, oppure no, ma certamente non è una risposta.
Oh, e per un modo portatile senza usare nessuno dei tanti strumenti specializzati per questo, puoi usare `awk`. Qualcosa come `awk -vRS ="> "-F '\ n' 'NR> 1 {print"> "$ 0 >> $ 1" .fa "}' file.fa`. Il tuo approccio `csplit` (e il tuo` split`) è davvero più semplice.


Questa domanda e risposta è stata tradotta automaticamente dalla lingua inglese. Il contenuto originale è disponibile su stackexchange, che ringraziamo per la licenza cc by-sa 3.0 con cui è distribuito.
Loading...