color / clrsvsim Goto Github PK

Currently, simulated reads are injected in the same coordinates as the original reads. However, in practice, the aligner may change the alignment position (e.g. soft-clipped reads, simulated SNPs/indels). Thus, we should add an option to re-align reads from the simulated BAM to make it more realistic.

The cigar package is used for generating the CIGAR string when simulating split reads. However, if the CIGAR string contains a deletion and the breakpoint is within the deleted region, then the CIGAR string no longer matches the length of the read (see below for example). While this looks like a bug in the cigar package, we're tracking it here as well in case we're not able to merge the change to the upstream repo in time and/or this is something the CIGAR API does not support.

This has the effect of generating an incorrect BAM file, where tools such as samtools fail to process the entire BGZF block where such an error is present. We recommend setting split_read_ratio=0 when running the simulation to avoid this bug.

Example (note how the length of the read changes after mask_left(10) is called):

>>> from cigar import Cigar
>>> cigar = Cigar('5M10D10M')
>>> len(cigar)
15
>>> mod = cigar.mask_left(4)
>>> mod
Cigar('4S1M10D10M')
>>> len(mod)
15
>>> mod_incorrect = cigar.mask_left(10)
>>> mod_incorrect
Cigar('20S5M')
>>> len(mod_incorrect)
25

Hi there,

I am trying to use clrsvsim to simulate some dataset to test our pipeline, after I made a python sandbox with all the prerequisites satisfied, I use python setup.py instal to install clrsvsim, no error or warnings occurred, following is the last few lines appeared on my terminal:

Using /share/home/kai.feng/env/ve-clrsvsim/lib/python2.7/site-packages
Finished processing dependencies for clrsvsim==0.0.1

after this, I try to use the command in README to test clrsvsim, and I cannot execute the command, is there a detailed readme about how to use clrsvsim or there are something wrong during my installation?

Thanks,

Kai

In the 'write_copies' function, the code for downsampling reads is as below:
# Probabilistically write the remaining fractional part of num_copies.
if (num_copies == .5 and hash(read.qname) % 2) or random.random() < num_copies:
outsam.write(read)
copies_written += 1

Why does it have this condition check '(num_copies == .5 and hash(read.qname) % 2)' specifically for ratio_change=0.5? This would result in ~75% of the original reads be written to the simulated BAM when ratio_change=0.5. Is it intentional? What is the biology behind this?

Thanks