У меня есть набор пар точек x для рисования сегментов вдоль оси x для создания пользовательской карты чтения в R:
Половина задачи построения этих сегментовопределяет их позиции y так, чтобы не было двух перекрывающихся сегментов на одном уровне y.Для каждого сегмента я перебираю уровни y из первой позиции, пока не доберусь до позиции, которая еще не содержит сегмент, который будет перекрывать текущий.Затем я записываю конечную позицию текущего сегмента и перехожу к следующему.
Фактический код представляет собой следующую функцию:
# Dummy data
# A list of start and end positions for each segment along the X axis. Sorted by start.
# Passing the function few.reads draws a map in half a second. Passing it many.reads takes about half an hour to complete.
few.reads <- data.frame( start=c(rep(10,150), rep(16,100), rep(43,50)), end=c(rep(30,150), rep(34,100), rep(57,50)) );
many.reads <- data.frame( start=c(rep(10,15000), rep(16,10000), rep(43,5000)), end=c(rep(30,15000), rep(34,10000), rep(57,5000)) );
#---
# A function to draw a series of overlapping segments (or "reads" in my along
# The x-axis. Where reads overlap, they are "stacked" down the y axis
#---
drawReads <- function(reads){
# sort the reads by their start positions
reads <- reads[order(reads$start),];
# minimum and maximum for x axis
minstart <- min(reads$start);
maxend <- max(reads$end);
# initialise yread: a list to keep track of used y levels
yread <- c(minstart - 1);
ypos <- c(); #holds the y position of the ith segment
#---
# This iteration step is the bottleneck. Worst case, when all reads are stacked on top
# of each other, it has to iterate over many y levels to find the correct position for
# the later reads
#---
# iterate over segments
for (r in 1:nrow(reads)){
read <- reads[r,];
start <- read$start;
placed <- FALSE;
# iterate through yread to find the next availible
# y pos at this x pos (start)
y <- 1;
while(!placed){
if(yread[y] < start){
ypos[r] <- y;
yread[y] <- read$end;
placed <- TRUE;
}
# current y pos is used by another segment, increment
y <- y + 1;
# initialize another y pos if we're at the end of the list
if(y > length(yread)){
yread[y] <- minstart-1;
}
}
}
#---
# This is the plotting step
# Once we are here the rest of the process is very quick
#---
# find the maximum y pos that is used to size up the plot
maxy <- length(yread);
miny = 1;
reads$ypos <- ypos + miny;
print("New Plot...")
# Now we have all the information, start the plot
plot.new();
plot.window(xlim=c(minstart, maxend+((maxend-minstart)/10)), ylim=c(1,maxy));
axis(3,xaxp=c(minstart,maxend,(maxend-minstart)/10));
axis(2, yaxp=c(miny,maxy,3),tick=FALSE,labels=FALSE);
print("Draw the reads...");
maxy <- max(reads$ypos);
segments(reads$start, maxy-reads$ypos, reads$end, maxy-reads$ypos, col="blue");
}
Мой фактический набор данных очень большой и содержит области, которые, насколько я могу судить, могут иметь до 600000 операций чтения.Чтения будут естественным образом накладываться друг на друга, поэтому очень легко реализовать сценарий наихудшего случая, когда все чтения накладываются друг на друга.Время, необходимое для построения большого количества операций чтения, для меня неприемлемо, поэтому я ищу способ сделать этот процесс более эффективным.Могу ли я заменить свои петли чем-то более быстрым?Есть ли алгоритм, который может организовать чтение быстрее?Я действительно не могу придумать лучшего способа сделать это в данный момент.
Спасибо за вашу помощь.