3 added 2967 characters in body
source | link

To make this customizable change:

sizes <- c(whatever, numbers, you, like)

output=NA (to a file path prefix to get a .csv file with the metadata)

genome = DNAStringSet() To replace the sample sequence

bpSeqSize = number (to change how many of the flanking region is shown in the metadata)

To make this customizable change:

sizes <- c(whatever, numbers, you, like)

output=NA (to a file path prefix to get a .csv file with the metadata)

genome = DNAStringSet() To replace the sample sequence

bpSeqSize = number (to change how many of the flanking region is shown in the metadata)

2 added 2967 characters in body
source | link

EDIT with RSVSim solution!

Finally found a solution based on RSVSim

genome = DNAStringSet(
  c(chr1 = "AAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTT",
    chr2 = "GGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCCCCC"))

length_seq_chr1 = width(genome[1]) - 10
length_seq_chr2 = width(genome[2]) - 10
sizes <- c(3,3,5,5)
start_chr1 <- c("list", length(sizes))
start_chr2 <- c("list", length(sizes))
end_chr1 <- c("list", length(sizes))
index = 1
for(i in sizes){
  newstartchr1 <- sample(1:length_seq_chr1,1)
  newstartchr2 <- sample(1:length_seq_chr2,1)
  start_chr1[[index]] <- newstartchr1
  start_chr2[[index]] <- newstartchr2
  end_chr1[[index]] <- newstartchr1 + (i-1)
  index <- index + 1}
start_chr1 <- as.integer(start_chr1)
end_chr1 <- as.integer(end_chr1)
start_chr2 <- as.integer(start_chr2)

knownInsertion2 = GRanges(IRanges(c(start_chr1),c(end_chr1)),seqnames="chr1", chrB="chr2", startB=c(start_chr2), copied=TRUE)
knownInsertion2

sim = simulateSV(output=NA, genome=genome, regionsIns=knownInsertion2, bpSeqSize=5, random=FALSE, verbose=FALSE)
sim
metadata(sim)

Gives as an result:

> knownInsertion2
GRanges object with 4 ranges and 3 metadata columns:
      seqnames    ranges strand |        chrB    startB    copied
         <Rle> <IRanges>  <Rle> | <character> <integer> <logical>
  [1]     chr1  [21, 24]      * |        chr2        30      TRUE
  [2]     chr1  [16, 19]      * |        chr2        21      TRUE
  [3]     chr1  [ 6, 11]      * |        chr2        10      TRUE
  [4]     chr1  [ 9, 14]      * |        chr2        27      TRUE
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths
> 
> sim = simulateSV(output=NA, genome=genome, regionsIns=knownInsertion2, bpSeqSize=5, random=FALSE, verbose=FALSE)
> sim
  A DNAStringSet instance of length 2
    width seq                                                                                                                           names               
[1]    40 AAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTT                                                                                      chr1
[2]    60 GGGGGGGGGAAAAAAGGGGGGGGGGGAAAACCCCCCAAAAAACCCTTTTCCCCCCCCCCC                                                                  chr2
> metadata(sim)
$insertions
  Name ChrA StartA EndA ChrB StartB EndB Size Copied BpSeqA BpSeqB_5prime BpSeqB_3prime
1    1 chr1     21   24 chr2     30   33    4   TRUE                 CCTT          TTCC
2    2 chr1     16   19 chr2     21   24    4   TRUE                 GGAA          AACC
3    3 chr1      6   11 chr2     10   15    6   TRUE                 GGAA          AAGG
4    4 chr1      9   14 chr2     27   32    6   TRUE                 CCAA          AACC

EDIT with RSVSim solution!

Finally found a solution based on RSVSim

genome = DNAStringSet(
  c(chr1 = "AAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTT",
    chr2 = "GGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCCCCC"))

length_seq_chr1 = width(genome[1]) - 10
length_seq_chr2 = width(genome[2]) - 10
sizes <- c(3,3,5,5)
start_chr1 <- c("list", length(sizes))
start_chr2 <- c("list", length(sizes))
end_chr1 <- c("list", length(sizes))
index = 1
for(i in sizes){
  newstartchr1 <- sample(1:length_seq_chr1,1)
  newstartchr2 <- sample(1:length_seq_chr2,1)
  start_chr1[[index]] <- newstartchr1
  start_chr2[[index]] <- newstartchr2
  end_chr1[[index]] <- newstartchr1 + (i-1)
  index <- index + 1}
start_chr1 <- as.integer(start_chr1)
end_chr1 <- as.integer(end_chr1)
start_chr2 <- as.integer(start_chr2)

knownInsertion2 = GRanges(IRanges(c(start_chr1),c(end_chr1)),seqnames="chr1", chrB="chr2", startB=c(start_chr2), copied=TRUE)
knownInsertion2

sim = simulateSV(output=NA, genome=genome, regionsIns=knownInsertion2, bpSeqSize=5, random=FALSE, verbose=FALSE)
sim
metadata(sim)

Gives as an result:

> knownInsertion2
GRanges object with 4 ranges and 3 metadata columns:
      seqnames    ranges strand |        chrB    startB    copied
         <Rle> <IRanges>  <Rle> | <character> <integer> <logical>
  [1]     chr1  [21, 24]      * |        chr2        30      TRUE
  [2]     chr1  [16, 19]      * |        chr2        21      TRUE
  [3]     chr1  [ 6, 11]      * |        chr2        10      TRUE
  [4]     chr1  [ 9, 14]      * |        chr2        27      TRUE
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths
> 
> sim = simulateSV(output=NA, genome=genome, regionsIns=knownInsertion2, bpSeqSize=5, random=FALSE, verbose=FALSE)
> sim
  A DNAStringSet instance of length 2
    width seq                                                                                                                           names               
[1]    40 AAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTT                                                                                      chr1
[2]    60 GGGGGGGGGAAAAAAGGGGGGGGGGGAAAACCCCCCAAAAAACCCTTTTCCCCCCCCCCC                                                                  chr2
> metadata(sim)
$insertions
  Name ChrA StartA EndA ChrB StartB EndB Size Copied BpSeqA BpSeqB_5prime BpSeqB_3prime
1    1 chr1     21   24 chr2     30   33    4   TRUE                 CCTT          TTCC
2    2 chr1     16   19 chr2     21   24    4   TRUE                 GGAA          AACC
3    3 chr1      6   11 chr2     10   15    6   TRUE                 GGAA          AAGG
4    4 chr1      9   14 chr2     27   32    6   TRUE                 CCAA          AACC
1
source | link

Answer eventually found on the BioConductor post: https://support.bioconductor.org/p/114354/#114585

It basically comes down to not using RSVSim, but other R packages.

Final code snippit:

genome = DNAStringSet(
  c(chr1 = "AAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTT",
    chr2 = "GGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCCCCC"))

sampleFromGenome <- function(genome, wd) {
  chrProb <- width(genome) / sum(width(genome))
  chr <- sample(names(genome), n, TRUE, prob = chrProb)
  start <- runif(n, 1, nchar(genome[chr]) - wd)
  GRanges(IRanges(start, width = wd), seqnames = chr)
}

n <- 2
wd <- sample(5:10, n, TRUE)
insertFrom <- sampleFromGenome(genome[1], wd)
insertAt <- sampleFromGenome(genome[2], 0)

#writeTable(insertFrom, stdout()) . # replace stdout with a file name

sequence <- genome[insertFrom]
at <- split(insertAt, seqnames(insertAt))[names(genome[2])]
atSequence <- split(sequence, seqnames(insertAt))[names(genome[2])]

replaceAt(genome[2], unname(ranges(at)), atSequence)

  A DNAStringSet instance of length 1
    width seq                                                                                                                           names               
[1]    54 GGGGGGGGGGGGGGGGGGTTTTTTTTTGGCCCCCCCCCCCCTTTTTCCCCCCCC                                                                        chr2