edit explaining python tag
I would still rather have a solution based on the RSVSim package in R, but while waiting for someone who might have an answer I wanted to look at other solutions as well.
If a result like what metadata(sim) produces can be achieved in python as well, I can do it without RSVSim. I need especially the numbers given there, so I do NOT need "Copied BpSeqA BpSeqB_5prime BpSeqB_3prime". The numbers need to reflect the original position of insert into chr2, because I need to reflect this back to the original (unedited) chr2.
end edit explaining pyton tag
I tried using the BioConductor forums for asking this question Question on BioConductor. However I do not seem to get an answer there and after some help with other tools on this website, i thought I could also ask the question here, in the hope that someone knows something about it.
I have a reference genome, and I need to put in some structural variations. It is really easy to do in RSVSim, and used it quite a lot already.
However, when processing INS, I found out that it just copies a DNA seq from one part, and puts it somewhere else in the same DNA seq (Bascially a DUP then), or if you specify another an extra DNAseq, to cut from, It randomly cuts from either chr1 or chr2 and puts it in the other. This creates deletions in the DNA seq I want to investigate, which is not what I want.
Does anyone know how to get this kind of result?
width seq names
[1] 40 AAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTT chr1
[2] 40 GGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCCCCC chr2
to this:
width seq names
[1] 32 AAAAAAAAAAAAAAAATTTTTTTTTTTTTTTT chr1
[2] 48 GGGGGGGGGGGGGGAAAATTTTGGGGGGCCCCCCCCCCCCCCCCCCCC chr2
And do it 30 times, randomly in the genome? (so take a seq from chr1 and put an insertion in chr2.
I tried this and more already for days, and I just cannot figure it out.
library(RSVSim)
seq_random <- readDNAStringSet("/path/random.fasta", "fasta")
seq_ref <- readDNAStringSet("/path/reference1MB.fasta", "fasta")
genome2 = DNAStringSet(c(seq_random, seq_ref))
names(genome2) = c("chr1","chr2")
y <- c(30, 30, 30, 30, 30, 50, 50, 50, 50, 50, 100, 100, 100, 100, 100, 500, 500, 500, 500, 500, 1000, 1000, 1000, 1000, 1000, 5000, 5000, 5000, 5000, 5000)
length_seq = width(genome2[2])
knownInsertion = GRanges(IRanges(0,length_seq), seqnames="chr2")
knownInsertion = GRanges(IRanges(0,length_seq), seqnames="chr1", chrB="chr2", chrA="chr1")
knownInsertion = GRanges(seqnames="chr1", chrB="chr2")
sim = simulateSV(output='/folder/of/output/', genome=genome2,
ins = 30, sizeIns=y, regionsIns=knownInsertion)
EDIT for comment section
The simple genome example:
genome = DNAStringSet(c("AAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTT","GGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCCCCC"))
names(genome) = c("chr1","chr2")
I believe I already tried those parameters. For example:
knownInsertion = GRanges(IRanges(16,25), seqnames="chr1", chrB="chr2")
sim = simulateSV(output=NA, genome=genome, ins = 1, regionsIns=knownInsertion, bpSeqSize=6, random=TRUE)
results in:
> sim
A DNAStringSet instance of length 2
width seq names
[1] 40 AAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTT chr1
[2] 40 GGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCCCCC chr2
> metadata(sim)
$insertions
Name ChrA StartA EndA ChrB StartB EndB Size Copied BpSeqA BpSeqB_5prime BpSeqB_3prime
1 insertion_1 chr1 3 12 chr1 6 15 10 FALSE AAAAA AAAAAA AAAAAA
so nothing happened at all. If I keep it simple, you get the result I mentioned before, where it is sometimes chr1->chr2 and sometimes chr2<-chr1
sim = simulateSV(output=NA, genome=genome, ins = 3, sizeIns = 5, bpSeqSize=6, random=TRUE)
> sim
A DNAStringSet instance of length 2
width seq names
[1] 40 AAAAAACCCCCAAAAATTTTTTTTTTTTTTAAAATTTTTT chr1
[2] 40 GGGGGGGGAAAAAGGGGGGGGGGGGCCCCCCCCCCCCCCC chr2
> metadata(sim)
$insertions
Name ChrA StartA EndA ChrB StartB EndB Size Copied BpSeqA BpSeqB_5prime BpSeqB_3prime
1 insertion_1 chr1 17 21 chr1 36 40 5 FALSE AAATTT TTTAAA AATTTT
2 insertion_2 chr1 1 5 chr2 9 13 5 FALSE GGGAAA AAAGGG
3 insertion_3 chr2 33 37 chr1 12 16 5 FALSE CCCCCC AAACCC CCCAAA
In this case I dont even know what happened exactly, because an uneven number of Ins, should make one of the chromosomes shorter or longer?
Finally, when setting copy on (so no cut/paste)
sim = simulateSV(output=NA, genome=genome, ins = 3, sizeIns = 5, bpSeqSize=6, random=TRUE, percCopiedIns=1.00)
sim
metadata(sim)
> sim
A DNAStringSet instance of length 2
width seq names
[1] 50 AAAATTTTTAAAAAAAAAAAAAAAATTTTTTATTTTTTTTTTTTTTTTTT chr1
[2] 45 GGGGGGGGGGGGGGGGGAAAAAGGGCCCCCCCCCCCCCCCCCCCC chr2
> metadata(sim)
$insertions
Name ChrA StartA EndA ChrB StartB EndB Size Copied BpSeqA BpSeqB_5prime BpSeqB_3prime
1 insertion_1 chr1 36 40 chr1 5 9 5 TRUE AAATTT TTTAAA
2 insertion_2 chr1 20 24 chr1 27 31 5 TRUE TTTATT TTTTTT
3 insertion_3 chr1 13 17 chr2 18 22 5 TRUE GGGAAA AAAGGG
Again. all the examples above have random chrA and chrB. Setting these specifically, only seems possible with the GRanges and IRanges, for which I tried tons of versions already.
It should be a very simple thing according to the manual, but actually doing it doesn't seem to get the desired results, while my version of Insertions is just the normal one in my opinion? (why would you want to create extra DELS or DUPS while making insertions?).
On their own support page, it is only possible to ask a question on the forums, but they never answer it (only look at it)
simulateSV? You don't get the desired result or something else? From the manual it seems to be possible. BTW knownInsertion should be chromosome 2 if you want all the insertions to be into the chromosome 2. But perhaps you are missing the optionrandom = TRUE$\endgroup$