在GenomicRanges中,一个有趣的问题是基因岛的鉴定。
我正在尝试找到范围的最大子集,其中相邻范围不超过特定距离。为了解决该问题,我尝试根据各个范围之间的差异来分配组。
我在IRanges软件包中搜索了合适的方法,但到目前为止我没有成功。
daf <- data.frame(seqnames="ConA",start=c(10,39,56,1000,5000),end=c(19,44,87,1100,5050),ID=paste0("Range",LETTERS[1:5]))
gr <- makeGRangesFromDataFrame(daf,keep.extra.columns=TRUE)
## Order the data frame by start column
oo <- order(daf$start)
daf <- daf[oo,]
# Calculate the distance
dd <- c(0,diff(daf$start))
daf$diff <- dd
daf$group <- rep(1,nrow(daf))
group <- 1
for(i in 1:nrow(daf)){
if(daf$diff[i] > 500){
group <- group + 1
}
daf$group[i] <- group
}
根据分配的组,可以找到最大的组。您是否知道任何更好的解决方案,避免使用for循环?
我相信这会实现您期望的结果,而无需循环:
## Load library
library(GenomicRanges)
## Create the data
gr <- GRanges(seqnames="ConA",
ranges=IRanges(start=c(10,39,56,1000,5000),
end=c(19,44,87,1100,5050),
names = paste0("Range",LETTERS[1:5])))
## Get the "start sites" = TSS.
## Use the promoters function if you need to take the strand into account:
grTSS <- promoters(gr, upstream = 0, downstream = 1)
## Create islands in which start sites are located less than 500bp apart
## If distance between start sites is <500bp then ranges of size 501bp centered on start sites overlap
islands <- reduce(grTSS + 250, ignore.strand = TRUE)
names(islands) <- paste0("island", 1:length(islands))
## Get the link between start sites and islands:
ovl <- findOverlaps(grTSS, islands, type = "within", ignore.strand = TRUE)
## Add the group info in the original GRanges object:
gr$group <- names(islands)[subjectHits(ovl)]
创建孤岛时可以考虑钢绞线信息:islands <- reduce(grstarts + 250, ignore.strand = FALSE)
不必只关注起始位点/ TSS,我们可以直接使用基因边界(起始和终止):islands <- reduce(gr + 250, ignore.strand = TRUE)
希望这会有所帮助