Re: pyslim VCF output nonconforming when reading treeseq

45 views

Skip to first unread message

Peter Ralph

unread,

Jun 14, 2021, 1:03:33 PM6/14/21

to slim-discuss

Follow-up here (I replied off-list, whoops): The problem is that
- the ALT allele you're seeing is the SLiM mutation ID of the SLiM
mutation. (more generally, it might be a comma-separated list of IDs,
if the mutations are stacking)
- the REF allele is always "" (the empty string), because anyone with
the ancestral state doesn't have any (SLiM) mutations.

To have nucleotides in those columns we need to edit the tree sequence
to swap out nucleotides for those SLiM IDs. This will be in the next
pyslim release, but in the meantime, here's code you can adapt for
your purpose:
https://github.com/tskit-dev/pyslim/issues/168
and in case anyone wants to have a preview of how it'll work in the
next pyslim (which I happened to be just working on), here it is:
https://github.com/tskit-dev/pyslim/pull/174

thanks,
Peter

On Sun, Jun 13, 2021 at 6:20 PM gasto...@gmail.com
<gasto...@gmail.com> wrote:
>
>
> Hi all,
>
> I am simulating a classic sweep by first simulating a single adaptive variant and overlaying mutations afterwards with pyslim. I output a VCF file, the neutral overlaid mutations look normal, but the "adaptive" mutation does not conform to VCF. The ALT allele is always missing, and the REF allele is some random number like 52 or 28. Here is an example of only the adaptive mutation.
>
> <pre style='color:#000000;background:#ffffff;'>##fileformat=VCFv4.2
> ##source=tskit 0.3.6
> ##FILTER=<ID=PASS,Description="All filters passed">
> ##contig=<ID=1,length=9999991>
> ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
> #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT p5:i4857 p4:i3461 p3:i2823 p2:i1176 p1:i968
> 1 4999990 . 51 . PASS . GT 1|1 0|0 1|1 0|0 0|0
>
> </pre>
>
>
>
> This is the code I am using for overlaying mutations:
>
>
> <pre style='color:#000000;background:#ffffff;'>import argparse
> import msprime, pyslim, gzip
> import numpy as np
> parser = argparse.ArgumentParser( description = 'Overlay neutral mutations into treeseq file and output VCF of last generation.')
> parser.add_argument("treeseq", type=str, help="trees file to overlay and write")
> parser.add_argument("output", type=str, help="output gzipped VCF file")
> parser.add_argument("sample_size", type=int, help="VCF sample size per population")
> args = parser.parse_args()
>
>
> ts = pyslim.load(args.treeseq)
> ts = msprime.mutate(ts, rate=1e-7, random_seed=1, keep=True)
> ts.dump("overlaid.trees")
> ts = pyslim.load("overlaid.trees")
> pops = np.unique(ts.individual_populations)
> individuals = np.array(0)
> populations = np.array(0)
> for i in pops:
> inds_sample = np.random.choice(ts.individuals_alive_at(0,population=i),args.sample_size,replace=False)
> pop_ids = np.full(len(inds_sample),i)
> # append
> populations = np.insert(populations,1,pop_ids)
> individuals = np.insert(individuals,1,inds_sample)
>
> individuals = np.delete(individuals,0)
> populations = np.delete(populations,0)
>
> ind_names = []
>
> for ind,pop in zip(individuals,populations):
> name = "p" + str(pop) + ":" + "i" + str(ind)
> ind_names.extend([name])
>
>
> with gzip.open(args.output, "wt") as f:
> ts.write_vcf(f,individuals=individuals,individual_names=ind_names)
> </pre>
>
> Any idea of what I might be missing?
>
> Thanks,
>
> gaston
>
>
>
> --
> SLiM forward genetic simulation: http://messerlab.org/slim/
> ---
> You received this message because you are subscribed to the Google Groups "slim-discuss" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to slim-discuss...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/slim-discuss/4ce59a75-02f8-45c9-9f77-9d9725028051n%40googlegroups.com.

Reply all

Reply to author

Forward

0 new messages