-
Notifications
You must be signed in to change notification settings - Fork 191
Simulating reads with vg sim
This page gives an overview of vg sim
and its usage. Note that this process is automated in toil-vg and we especially recommend using toil-vg if working with large graphs.
Most of the time we might want to simulate reads from one of the sample that was used to create the graph. To do this we need to index each haplotype using the GBWT. This index will be used to identify the paths in the graph corresponding to the haplotype or sample we want to simulate the reads from.
Assuming that x.fa contains reference sequence for a chromosome x and variants.vcf.gz contains variants for multiple samples including the sample for which we want to simulate reads (later called SAMP).
vg construct -r x.fa -v variants.vcf.gz -a -f -m 32 > graph.vg
vg index -x graph.xg -G graph.gbwt -v variants.vcf.gz graph.vg
Before calling vg sim
we'll create a new graph with only paths from the desired haplotype.
Haplotypes are represented as threads in the graph. Threads are prefixed by _thread_<SAMPLE>_<CHR>_<HAPLOTYPE>
with <HAPLOTYPE>
being 0 or 1 for diploid samples.
To create a graph with only threads from sample SAMP, chromosome x and haplotype 0:
vg mod -D graph.vg > thread.merge.vg
vg paths --gbwt graph.gbwt --extract-vg -x graph.xg -Q _thread_SAMP_x_0 >> thread.merge.vg
vg mod -N thread.merge.vg > thread_SAMP_0.vg
Now we can index this new graph and simulate, for example, a thousand 100 bp reads:
vg index -x thread_SAMP_0.xg thread_SAMP_0.vg
vg sim -x thread_SAMP_0.xg -l 100 -n 1000 -a | vg view -X - > SAMP.fastq
By default vg sim
outputs reads in a text format. Using -a
it outputs reads in the GAM format that we can then convert into FASTQ using vg view -X
.
Repeat the process above for both haplotypes and concatenate the output FASTQ files.
We can specify the fragment size and its standard deviation using the -p
and -v
arguments. For example:
vg sim -x thread_SAMP_0.xg -l 100 -n 1000 -a -p 500 -v 50 | vg view -X - > SAMP.fastq
A FASTQ file can be used to model the error profile using the -F
argument. For example:
FASTQ=ftp://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/NA12878/NIST_NA12878_HG001_HiSeq_300x/131219_D00360_005_BH814YADXX/Project_RM8398/Sample_U5a/U5a_AGTCAA_L002_R1_007.fastq.gz
vg sim -x thread_SAMP_0.xg -n 1000 -a -p 500 -v 50 -F $FASTQ | vg view -X - > SAMP.fastq
Note that in that case the read length argument (-l
) is not necessary.