3
$\begingroup$

I have asked this question at BioStars but have not gotten any suggestions so far. So, I am asking the community here, perhaps someone here has an idea how I can solve my problem.

I am dealing with a bunch of inversions and running some analysis on each of them. One of them is located on chromosome one, the starting position is 6199548 and ending position is 9699260. So the length of this variant is 3499712 bp. I want to generate non-overlapping windows based on the length of these focal positions (3499712) and assign the focal position to one of these windows that overlap with the starting and ending position of my focal position. I am achieving this by using the GenomicRanges::tileGenome function in R:

seqlength <-c(chr1=159217232,chr2=184765313,chr3=182121094,chr4=221926752,chr5=187413619)
bin_length<-c(3499712) 
 bins<-tileGenome(seqlength,tilewidth=bin_length,cut.last.tile.in.chrom
        = T)
            #bins_u<-unlist(bins)
            write.table(bins,file=paste("bins_good_",bin_length,sep=""),col.names
        = FALSE, row.names = TRUE) }

This generates non-overlapping 3499712 bps windows across the genome.

window_id chromosome start end length
"1" "chr1" 1 3499712 3499712
"2" "chr1" 3499713 6999424 3499712 
"3" "chr1" 6999425 10499136 3499712 
"4" "chr1" 10499137 13998848 3499712 
"5" "chr1" 13998849 17498560 3499712

However, I have a major problem here that I have not been able to sort out. I want one of these windows to overlaps exactly with the starting and ending position of my focal position (6199548:9699260). Only the focal position not any other position outside this region. This means 2 windows (one before and one after the focal position) will have a different length which my analysis is fine with. I just do not know how to do it?

This is my desired output:

window_id chromosome start end length
"1" "chr1" 1 3499712 3499712 "*"                  
"2" "chr1" 3499713 6199547 2699834 "*"
"3" "chr1" 6199548  9699260 3499712 "*"
"4" "chr1" 9699261 13998848 4299587 "*"
"5" "chr1" 13998849 17498560 3499712 "*"

As you see, the the second window is now smaller, the third window accommodates the focal position and the forth window is a little larger than the bin size. Does anyone have any idea how can I do this task in R?

$\endgroup$
1
  • $\begingroup$ substr(chromosome, start pos, finish pos) ? $\endgroup$ – M__ Mar 7 '19 at 17:19
5
$\begingroup$

I'm not sure how you did your calculation but for me it looks like your region is 3499713 bp in length.

site <- GRanges("chr1:6199548-9699260")

seqlength <- c(chr1=159217232,chr2=184765313,chr3=182121094,chr4=221926752,chr5=187413619)

bin_length <- width(site)

## check length of region
print(bin_length)
##3499713

bins <- tileGenome(seqlength, tilewidth=bin_length, cut.last.tile.in.chrom=TRUE)

## extract overlapping bins
flank.bins <- subsetByOverlaps(bins,site)

## find new coordinates excluding region of interest
new.flank.bins <- setdiff(reduce(flank.bins),site)

## combine new bins with your region of interest
new.bins <- c(new.flank.bins,site)

## remove old bins
bins <- subsetByOverlaps(bins,new.bins,invert=T)

## add new bins
bins <- sort(c(bins,new.bins))

## check that it makes sense
##table(width(bins[seqnames(bins) == 'chr1']))

write.table(bins,file=paste("bins_good_",bin_length,sep=""),col.names= FALSE, row.names = TRUE)
$\endgroup$
3
  • $\begingroup$ I sincerely appreciate your help. I did not know that the GenomicRanges options! I am now working on it to learn about the many other useful functions. $\endgroup$ – Anna1364 Mar 8 '19 at 0:21
  • $\begingroup$ I have been digging into your code, but I guess it does not give the outputs that I want. This is the output that I get from your code:"1" "chr1" 1 3499713 3499713 "" "2" "chr1" 3499714 10499139 6999426 "" "3" "chr1" 6199548 9699260 3499713 "" "4" "chr1" 10499140 13998852 3499713 "". this is what I get , the focal window(window 3) overlaps with window 2! $\endgroup$ – Anna1364 Mar 8 '19 at 21:45
  • $\begingroup$ If you copy and paste my code exactly and look at the first 5 lines of the resulting file, you should get the following: "1" "chr1" 1 3499713 3499713 "" "2" "chr1" 3499714 6199547 2699834 "" "3" "chr1" 6199548 9699260 3499713 "" "4" "chr1" 9699261 10499139 799879 "" "5" "chr1" 10499140 13998852 3499713 "*" If you share your code I can help you further. $\endgroup$ – story Mar 10 '19 at 13:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.