4

Rでカスタム読み取りマップを作成するためにx軸に沿ってセグメントを描画するためのxポイントのペアのセットがあります。

地図を読む例

これらのセグメントをプロットするタスクの半分は、オーバーラップする2つのセグメントが同じ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回の読み取りが可能な領域が含まれています。読み取りは自然に互いに重なり合うため、すべての読み取りが互いに重なり合うという最悪のシナリオを簡単に実現できます。大量の読み取りをプロットするのにかかる時間は私には受け入れられないので、プロセスをより効率的にする方法を探しています。ループをもっと速いものに置き換えることはできますか？読み取りをより迅速に調整できるアルゴリズムはありますか？現時点では、これを行うためのより良い方法を考えることはできません。

ご協力いただきありがとうございます。

4

2 に答える 2

2

貪欲な方法で各yレベルを埋めます。レベルが一杯になった後、1レベル下に移動し、決して上に戻らないでください。

擬似コード：

 y <- 1
 while segment-list.not-empty
   i <- 1
   current <- segment-list[i]
   current.plot(y)
   segment-list.remove(i)
   i <- segment-list.find_first_greater(current.end)
   while (i > 0)
     current <- segment-list[i]
     current.plot(y)
     segment-list.remove(i)
   y <- y + 1

これは、必ずしも「最適な」プロットを生成するわけではありませんが、少なくともO（n log n）です。

于 2012-03-26T13:20:23.867 に答える
1

開始値でソートできませんか？次に、リストを前から後ろに移動します。各アイテムについて、それをプロットしてから、リストの残りの部分で、プロットしたばかりのアイテムの終了座標よりも大きい最初のアイテムを2分探索します。何も見つからない場合は、Yを増やします。プロットするときに各アイテムを削除します。

ソートはO（N lg N）であり、各アイテムの二分探索はO（lg N）であるため、合計はO（N lg N）になります。

于 2012-03-26T17:00:14.277 に答える