Блог

Мне нужно найти длину перекрывающейся области на одних и тех же хромосомах между 2 группами(gp1 и gp2). (аналогичный вопрос в stackoverflow отличался от моей цели, потому что я хочу найти перекрывающуюся область, а не ИСТИННЫЙ/ЛОЖНЫЙ ответ).

Например:

 gp1: 
chr   start   end   id1 
chr1  580     600    1
chr1  900     970    2
chr3  400     600    3
chr2  100     700    4
 

 gp2:
chr   start   end   id2
chr1  590     864   1
chr3  550     670   2
chr2  897     1987  3
 

Я ищу способ сравнить эти 2 группы и получить такие результаты:

 id1   id2    chr   overlapped_length
 1     1     chr1     10
 3     2     chr3     50   
 

Должно указать вам правильное направление:

Загрузка библиотек

 # install.packages("BiocManager")
# BiocManager::install("GenomicRanges")
library(GenomicRanges)
library(IRanges)
 

Генерировать данные

 gp1 <- read.table(text = 
"
chr   start   end   id1 
chr1  580     600    1
chr1  900     970    2
chr3  400     600    3
chr2  100     700    4
", header = TRUE)

gp2 <- read.table(text = 
                      "
chr   start   end   id2
chr1  590     864   1
chr3  550     670   2
chr2  897     1987  3
", header = TRUE) 
 

Вычисление диапазонов

 gr1 <- GenomicRanges::makeGRangesFromDataFrame(
    gp1, 
    seqnames.field = "chr", 
    start.field = "start", 
    end.field = "end"
)
gr2 <- GenomicRanges::makeGRangesFromDataFrame(
    gp2, 
    seqnames.field = "chr", 
    start.field = "start", 
    end.field = "end"
)
 

Вычисление перекрытий

 hits <- findOverlaps(gr1, gr2)
p <- Pairs(gr1, gr2, hits = hits)
i <- pintersect(p)
 

Результат

 > as.data.frame(i)
  seqnames start end width strand  hit
1     chr1   590 600    11      * TRUE
2     chr3   550 600    51      * TRUE