Сохранение наибольшей строки из перекрытий строк в R

Question

У меня есть объект R data.frame из 5131 строки с четырьмя столбцами, образец которого выглядит следующим образом:

1 X  12830000  12910000  C
2 X  12960000  13510000  C
3 X  13525000  13675000  C
4 X  13670000  13715000  C
5 X  13670000  13770000  E2
6 X  13670000  14050000  E3
7 X  13765000  14050000  E1
8 X  13910000  14050000  E1
9 X  13940000  14050000  C
10 X  15360000  15590000  E3

Это геномный файл .bed, столбец 2 - координата начальной позиции отрезка ДНК, столбец 3 - конечная позиция, а столбец 4 - метаинформация об этом. Так что это можно представить как кусочки струн, которые пересекаются друг с другом в соответствии с их координатами.

Я хочу написать скрипт R, который выполняет следующее:

Если есть перекрытия между рядами более 40000 между собой, выберите строку с самой длинной длиной среди перекрытий. В выборочных данных , последовательно перемещающихся вниз, имеется более 40000 перекрытий между строками 4 и 5, строками 5 и 6, строками 6 и 7, строками 7 и 8 и строками 8 и строками 9, среди которых ряд 6 является самым длинным. Среди 4 и 5 строка 5 самая длинная, поэтому я держу 5. Между строк 5 и 6 самая длинная 6, поэтому я держу 6. Я продолжаю делать это пока я не найду ряд с перекрытием менее 4000 и не выберу самый длинный из них.

Так что мой новый фрейм данных должен выглядеть примерно так:

1 X  12830000  12910000  C
2 X  12960000  13510000  C
3 X  13525000  13675000  C
6 X  13670000  14050000  E3  <- (keeping the longest one)
10 X  15360000  15590000  E3

До сих пор я пытался сделать следующее:

output_4 <- fr_t[NULL,] 
for(i in 1:nrow(fr_t)-1){ 
if(isTRUE(fr_t[i+1,]$V3-fr_t[i,]$V2<40000||fr_t[i,]$V3-fr_t[i+1,]$V2 < 40000)) {
if(isTRUE((fr_t[i+1,]$V3-fr_t[i+1,]$V2)>(fr_t[i+2,]$V3-fr_t[i+2,]$V2))){
next
}
next
}
output_4 <- rbind(output_4, fr_t[i+2,])
}  #fr_t is my original dataframe

Я не могу понять, как выполнить итерацию моего i до тех пор, пока строки не пересекаются, чтобы сохранить самую длинную строку.

Кроме того, как я могу сохранить метаинформацию в столбце 4 в виде матрицы перекрытий? Например, новый фрейм данных следующий:

E3  C  E2  E1  E1  C

Первый столбец - это столбец 4 самой большой строки, а все остальные столбцы - это столбец 4 перекрытий внутри него? Это потребует объекта с различными номерами столбцов. Спасибо за терпение.

Answer 1

Прежде всего, я полностью согласен с комментарием @ zacdav. Посмотрите на пакет R / Bioconductor GenomicRanges; он был разработан именно для такого рода операций. Существуют отличные учебники, например, здесь и здесь .

Что касается вашего вопроса, то ожидаемый результат воспроизводит следующее:

# Your sample data as a data.frame
df <- read.table(text =
    "1 X  12830000  12910000  C
2 X  12960000  13510000  C
3 X  13525000  13675000  C
4 X  13670000  13715000  C
5 X  13670000  13770000  E2
6 X  13670000  14050000  E3
7 X  13765000  14050000  E1
8 X  13910000  14050000  E1
9 X  13940000  14050000  C
10 X  15360000  15590000  E3", header = F, row.names = 1)

# Convert data.frame to GRanges
library(GenomicRanges);
gr <- with(df, GRanges(
    seqnames = df[, 1],
    IRanges(start = df[, 2], end = df[, 3]),
    id = df[, 4]))

# Find overlapping regions within gr
hits <- findOverlaps(gr, gr, minoverlap = 40000);

# Remove self-overlapping hits
hits <- hits[queryHits(hits) != subjectHits(hits)];

# Determine features that are shorter than the overlapping feature
mcols(hits)$queryWidth = width(gr[queryHits(hits)]);
mcols(hits)$subjectWidth = width(gr[subjectHits(hits)]);
mcols(hits)$hit <- ifelse(
    mcols(hits)$queryWidth < mcols(hits)$subjectWidth, 
    queryHits(hits), 
    subjectHits(hits));

# Remove those shorter overlapping features
gr.final <- gr[-unique(mcols(hits)$hit)];
gr.final;
#GRanges object with 5 ranges and 1 metadata column:
#      seqnames               ranges strand |       id
#         <Rle>            <IRanges>  <Rle> | <factor>
#  [1]        X [12830000, 12910000]      * |        C
#  [2]        X [12960000, 13510000]      * |        C
#  [3]        X [13525000, 13675000]      * |        C
#  [4]        X [13670000, 14050000]      * |       E3
#  [5]        X [15360000, 15590000]      * |       E3

Если вы хотите преобразовать gr.final обратно в data.frame, вы можете использовать as.data.frame(gr.final).