7
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ Why do you have 4 fields in sap_ids for chr17 and 3 for chr1? Why does that change? Do your positions increase with the chromosome numbers (are all positions on chr2 greater than those on chr1) or do you count each chromosome separately (so each chromosome starts from position 1)? If the latter, it makes sense that your approach fails. $\endgroup$ – terdon Nov 24 '17 at 10:41
  • $\begingroup$ I do not get the part with 15, 20 and 30 kb windows. You want the windows to be of various sizes dependent of number of SNPs inside or you want to have one window size that will correspond on average to ~200 SNPs? $\endgroup$ – Kamil S Jaron Nov 24 '17 at 16:00
  • $\begingroup$ @terdon, I just edited my post slightly. see my edits. yes, actually in my case I count each chromosome separately which means each chromosome starts from position1. Any idea how to do it in R? @ Kamil S Jaron, I just want to do my analysis with different window sizes! that's it. $\endgroup$ – Anna1364 Nov 24 '17 at 16:32
  • $\begingroup$ @Anna To avoid removing the answers I rolled back to the previous version. If your question is not fully answered please edit it or comment explaining why not. You can always deselect the accepted answer. $\endgroup$ – llrs Nov 27 '17 at 16:05
  • 1
    $\begingroup$ Please stop radically changing your question! If you need more help, ask a new question explaining what else you need. But since this has been answered, changing it renders the answers obsolete. $\endgroup$ – terdon Nov 28 '17 at 11:43
6
$\begingroup$

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
$\endgroup$
  • 1
    $\begingroup$ Relatedly, there's the tile() command from GenomicRanges that will produce a GRangesList object of the input. It won't produce exactly what was desired by OP, but it may be a more convenient alternative. $\endgroup$ – Devon Ryan Nov 25 '17 at 16:07
  • $\begingroup$ Thanks so much @Sebastian Müller. yes. This is actually a better approach than writing code for that. $\endgroup$ – Anna1364 Nov 25 '17 at 21:52
  • $\begingroup$ Glad it is useful! Maybe you could ask even ask your SNP question in a different question and take it out here so the the actual question becomes more in-line with the title (which is almost self-explanatory by itself)? $\endgroup$ – Sebastian Müller Nov 26 '17 at 8:47
  • $\begingroup$ Also for completeness, there is also the GenomicRanges::slidingWindows function which takes a GRanges object designed to results in overlapping windows, but can be made to produce non-overlapping windows as well by setting width = step $\endgroup$ – Sebastian Müller Nov 26 '17 at 9:11

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.