Simulating a prokaryotic genome

[aprice67]

I want to simulate a genome, then simulate RNA-seq reads from that genome that will map with varying degrees of multi-mapping. So for instance, I want 10% of reads to map more than once, 20% of reads to map more than once, 30%, etc. This is a prokaryotic genome, so no need to worry about indels or the like. Plan to do the mapping with bowtie2. So far I haven't been able to find any papers that do this or tools that can do this. Can anyone offer some advice or point me in a good direction? Thank you in advance!

[GenoMax]

Why simulate one when you have thousands of prokaryotic genomes to choose from in Genbank?

[aprice67]

I want to show if different levels of reads mapping to multiple locations has a specific effect in a very controlled way. I could go and find a bunch of different genomes and simulate rna seq reads until I found reads with 10, 20, 30% of positions mapped to multiple positions, but I don't have any way of knowing even how many multi-mapping positions there would be in any given genome.

If there is a way I can tell in advance roughly how the proportion of multi-mapped reads there would be for a given genome though, using existing genomes might be a solid workaround. Is there any way to know or predict or control for that?

[GenoMax]

Isn't this going to be partially a characteristic of a particular library? If a library happens to have areas of genome over-represented (e.g. because of too many PCR cycles) then the multimapper problem will be acute but with a different prep the problem may be absent.

Not sure what you exactly mean by "different levels of reads mapping to multiple locations has a specific effect". What sort of effect?

[aprice67]

So I have some real datasets and I'm looking at false discovery of differentially expressed genes. I'm particularly looking at how using closely related reference genomes (but not exactly the one the reads were generated from) influences the false discovery rate. So I've found that some genes are incorrectly identified as DE when there are SNPs, which I would expect, but some genes that have little to no difference are incorrectly identified and I suspect its because of reads mapping to multiple locations and how aligners handle those reads.

So for the data sets I have, one batch had around 10-20% multi-mapped genes, and the other was much higher, like 60-70%. The dataset with fewer multi-mapped reads had significantly fewer false positives. So I want to see if the numbers of false positives I can identify increases as I increase the number of reads that map to multiple locations. I feel like having that sort of analysis in addition to what I've done with my actual data will strengthen the paper I'm working on.

It's the only way I can think of at the moment to see if the mapping is the cause of the false positives I can't predict based on SNPs alone.

[GenoMax]

Which aligner are you using and how it (or you) are handling your multi-mappers? 60-70% seems pretty high. Are reads of equal length in all cases?

I will refer to the BBMap's handling of multi-mappers below

best (use the first best site)
toss (consider unmapped)
random (select one top-scoring site randomly)
all (retain all top-scoring sites)

[aprice67]

Bowtie2. I also think it's pretty high. It's a publicly available dataset I found through a paper. As far as I know reads are all the same length.