Spørsmål:
Simulering av DNA-sekvensutvikling i R
compbiostats
2018-10-04 00:19:20 UTC
view on stackexchange narkive permalink

Jeg håper noen kan låne tankene sine om koden nedenfor for å generere DNA-sekvenser under Kimura-2-Parameter-modellen for DNA-erstatning.

Problemet er at hver gang koden kjøres og haplotypen distribusjon blir undersøkt, det er alltid en veldig skjev fordeling.

Det jeg ønsker er at det genereres en annen distribusjon hver gang koden kjøres. det vil si at det skal være flere eksemplarer som deler den samme haplotypen (i stedet for at de fleste deler den vanligste haplotypen).

Her er koden min: 

  bibliotek (pegas) set.seed (17) sim.seqs <- TRUElength.seqs <- 500num.seqs <- 100 # nummer av DNA sequencessubst.model <- "K80" # nucleotide substitution modeltransi.rate <- 1e-4 # transition ratetransv.rate <- transi.rate / 2 # transversion rateif (sim.seqs == TRUE) {nucl <- as. DNAbin (c ('a', 'c', 'g', 't')) res <- prøve (nucl, size = length.seqs, erstatt = TRUE, prob = rep (0,25, 4)) if (subst .modell == "K80") {transi.set <- list ('a' = as.DNAbin ('g'), 'c' = as.DNAbin ('t'), 'g' = as.DNAbin ( 'a'), 't' = as.DNAbin ('c')) transv.set <- liste ('a' = as.DNAbin (c ('c', 't')), 'c' = as .DNAbin (c ('a', 'g')), 'g' = as.DNAbin (c ('c', 't')), 't' = as.DNAbin (c ('a', ' g '))) transi <- funksjon (res) {unlist (transi. set [as.character (res)])} transv <- function (res) {sapply (transv.set [as.character (res)], sample, 1)} duplicate.seq <- function (res) {num. transi <- rbinom (n = 1, størrelse = lengde.seqs, prob = transi.rate) # totalt antall overganger hvis (num.transi > 0) {idx <- prøve (lengde.seqs, størrelse = num.transi, erstatt = FALSE) res [idx] <- transi (res [idx])}
num.transv <- rbinom (n = 1, størrelse = length.seqs, prob = transv.rate) # totalt antall transversjoner hvis (num.transv > 0) {idx <- sample (length.seqs, size = num. transv, erstatt = FALSE) res [idx] <- transv (res [idx])} res}} res <- matrix (replikere (num.seqs, duplicate.seq (res)), byrow = TRUE, nrow = num. seqs) klasse (res) <- "DNAbin" # write.dna (res, file = "seqs.fas", format = "fasta") h <- sort (haplotype (res), avtagende = SANT, hva = "frekvenser ") navn (h) <- 1: nrow (h)}

Output 

  hHaplotypes extracted from: resAntall haplotyper: 5 Sekvenslengde: 500Haplotype etiketter og frekvenser: 1 2 3 4596 1 1 1 1

Hvis koden kjøres flere ganger (uten frøet), vises det samme mønsteret: h er alltid skjev mot den mest dominerende haplotype.

Jeg vil for eksempel kunne kjøre koden (uten seed) og få noe sånt som 

  hHaplotypes ekstrahert fra: resAntall haplotyper: 5 Sekvenslengde: 500Haplotype-etiketter og frekvenser: 1 2 3 4 5 35 25 20 15 5

jeg tror jeg må spesifisere en fordeling for haplotypene, men antall haplotyper generert av mutasjonsprosessen er ikke kjent på forhånd.

Noen tanker?

I likhet med hva gringer sa: Kanskje jeg forsto feil modellen din, men hvorfor forventer du bare en overgang og omstilling? Bør ikke det være `rbinom (n = 500, ...`? I en sekvens på 500 forventer du å vite hvor mange substitusjoner som er i hele sekvensen, ikke hvis det er en mutasjon på sekvensen.
Hei Jarrett, det tomme dobbeltlinjeavstanden gjør koden din dessverre vanskeligere å lese enn nødvendig. Det ville hjelpe hvis du bare brukte tomme linjer for å logisk skille blokker, ikke mellom alle linjer.
@KonradRudolph Fin fangst! Jeg har fikset linjeavstanden.
@Llopis Antall erstatninger per sekvens vil ikke være kjent på forhånd. Dette bestemmes med rbinom (...).
En svar:
gringer
2018-10-04 01:54:34 UTC
view on stackexchange narkive permalink

Hva forventer du at disse linjene skal gjøre? 

  num.transi <- rbinom (n = 1, size = length.seqs, prob = transi.rate) ... num. transv <- rbinom (n = 1, size = length.seqs, prob = transv.rate)

Ved å bruke variablene som er spesifisert i koden din, er de ekstremt nullskjevte. Dette er ikke så overraskende gitt mutasjonshastigheten og sekvenslengden.

Jeg forventer at økning av sekvenslengden, eller økning av overgangshastigheten, bør føre til at flere haplotyper genereres.

kjører dette gjennom meg selv, ser jeg at å øke mutasjonsfrekvensen øker antall haplotyper, men skjevheten beholdes. Dette er fordi hver mutasjon skaper en ny haplotype, etter noe som en uendelig allelmodell.

Jeg antar at dette ikke er spesielt overraskende. Det er lite sannsynlig at to tilfeldige mutasjoner fra en bestemt haplotype vil resultere i den samme haplotypen, og derfor vil eventuelle mutasjoner (fra hvilken som helst eksisterende haplotype) generelt skape noe nytt.

Disse linjene beregner antall overganger og transversjoner over sekvensjusteringen av en gitt lengde i henhold til en binomial fordeling. Du har rett i at flere haplotyper sannsynligvis vil bli generert når lengden. Kvoter og / eller transi.rate økes. Imidlertid ønsker jeg å sikre at det er flere sekvenser per haplotype (slik at fordelingen ikke er så skjev).
Du klager på en skjevhet mot den primære haplotypen, dvs. at det ikke er nok mutasjoner. Min anbefalte løsning er å enten øke sekvenslengden eller øke overgangshastigheten.
Har lekt med koden ... får fortsatt den samme effekten når mutasjonshastighet og / eller sekvenslengde økes. Det burde være en bedre måte - jeg vil tenke på det og eksperimentere litt mer.
Min spesifikke anvendelse er i DNA-strekkoding, som i seg selv er avhengig av korte (~ 500 - ~ 650 bp) DNA-fragmenter. Å øke lengden. Kvm er egentlig ikke et alternativ - med mindre jeg fungerer, trimmer jeg justeringen jevnt i både 5 'og 3' ender til ønsket lengde. Jeg grubler på dette.


ⓘ
Denne spørsmålet ble automatisk oversatt fra engelsk.Det opprinnelige innholdet er tilgjengelig på stackexchange, som vi takker for cc by-sa 4.0-lisensen den distribueres under.
Loading...