dividing genome into non-overlapping windows using R

Question

I have nearly 10 million SNPs located on 10 chromosomes. I want to divide the genome into non-overlapping windows of 15, 20 and 30 kb. Here is part of my SNP table:

head (sap_ids)
       snp_id    chr  pos 
    Chr01__15043   1 15043    
    Chr01__15079   1 15079   
    Chr01__15139   1 15139    
    Chr01__15165   1 15165    
    ...
    ...
    ...
    Chr17__214708424  17 214708424        
    Chr17__214708451  17 214708451        
    Chr17__214708484  17 214708484        
    Chr17__214708508  17 214708508        
    Chr17__214708574  17 214708574        

I have been using this code, but it gives me wrong output because each chromosome starts from position 1

win_size<-c(15000,30000,50000,100000)
res<- cbind(snp_ids,
       data.frame(lapply(setNames(win_size, paste("window",win_size, sep = "_")),
       function(x)as.numeric(ceiling(snp_ids$pos/x)))))

For instance, as I you see in the example below in window 4 I get SNPs from chromosomes 1, 3 and 5:

snp_id     chrom poistion   window
Chr01__58332   1 58332            4
Chr01__58335   1 58335            4
Chr01__58341   1 58341            4
Chr01__58450   1 58450            4
Chr01__58471   1 58471            4
Chr01__58530   1 58530            4
Chr01__58542   1 58542            4
Chr01__58641   1 58641            4
Chr03__45457   3 45457            4
Chr03__45604   3 45604            4
Chr04__56873   4 56873            4
Chr04__57387   4 57387            4
Chr04__57399   4 57399            4
Chr04__57528   4 57528            4
Chr04__58419   4 58419            4
Chr04__59670   4 59670            4
Chr04__59704   4 59704            4

I guess the solution would be to run the code above by looping through each chromosome I tried something like this:

for (i in 1:10){ 
       res <- cbind(cats,
             data.frame(
               lapply(setNames(win_size, paste("window", win_size, sep = "_")), function(w)
                 paste(cats$chr[cats$chr == i], ceiling(cats$pos/w), sep = "_"))
               ))
       }

But still not working appropriately (column 4 should be 1_1, for the chr 1)!

Chr01__912    1  912          2_1
Chr01__944    1  944          2_1
Chr01__1107   1  1107         2_1
Chr01__1118   1  1118         2_1

Can anyone help me to figure out how can I modify the code to create these non-overlapping windows properly?

Answer 1

The question is a bit confusing to me, at it's core I understand you want non-overlapping windows across chromosomes.

One way to achieve is to use GenomicRanges::tileGenome function, which needs the chromosome lengths as input, e.g.

library(Biostrings)
mygenome <- readDNAStringSet(list.files(mypath,"mygenome.fa$",full=TRUE))
chrSizes <- width(mygenome)
names(chrSizes) <- names(mygenome)
print(chrSizes)
#     Chr1     Chr2     Chr3     Chr4     Chr5 
# 18585056 19698289 23459830 26975502 30427671 

This can than be fed into tileGenomes:

library(GenomicRanges)
bins   <- tileGenome(chrSizes, tilewidth=5e5, cut.last.tile.in.chrom=T)

This creates a GRange object which gives you a whole repertoire of post processing options, i.e. overlapping with SNPs or NGS reads (I use the latter for genome wide coverage plots):

print(bins)
# GRanges object with 240 ranges and 0 metadata columns:
#         seqnames               ranges strand
#            <Rle>            <IRanges>  <Rle>
#     [1]     Chr1   [      1,  500000]      *
#     [2]     Chr1   [ 500001, 1000000]      *
#     [3]     Chr1   [1000001, 1500000]      *
#     [4]     Chr1   [1500001, 2000000]      *
#     [5]     Chr1   [2000001, 2500000]      *
#     ...      ...                  ...    ...
#   [236]     Chr5 [28000001, 28500000]      *
#   [237]     Chr5 [28500001, 29000000]      *
#   [238]     Chr5 [29000001, 29500000]      *
#   [239]     Chr5 [29500001, 30000000]      *
#   [240]     Chr5 [30000001, 30427671]      *
#   -------
#   seqinfo: 5 sequences from an unspecified genome