1

大きなデータセットがあり、applyで使用するカスタムマージ関数を作成したいのですが、特定の問題を解決できません。時間がかかりすぎるのでループは使えません。データは大まかに次のようになります。

#     [ Name, Strand, Start, End ]

R1 = c( 'GeneA', '+', 1000, 1500 )
R2 = c( 'GeneA', '+', 1510, 2000 )
R3 = c( 'GeneA', '+', 2001, 2500 )
R4 = c( 'GeneB', '-', 3100, 4000 )

データは、行R1：R4のdata.frameです。

これまでのところ、RiとRj（j = i +1）を比較し、Nameが同じで、Strandが同じで、それらの間のギャップが100未満の場合、それらをマージする関数を取得できます。

GAP = Ri[End] - Rj[Start]

関数を各行に適用して出力を作成するとします。次に、関数の出力で作成する必要があります

R1 = c( 'GeneA', '+', 1000, 2500 )
R4 = c( 'GeneB', '-', 3100, 4000 )

applyを使用して2つの連続する要素を1つにマージできる作業関数を取得できますが、3つの連続する要素をマージする方法がわかりません。醜い解決策の1つは、追加のマージが行われなくなるまで2つの連続する要素のマージ関数を実行することでしたが、これは効率的ではありません。私は少しアイデアにこだわっています、そしてどんな洞察もいただければ幸いです。

編集：明確にするために、データは連続する染色体上の開始位置（図示せず）の順に並べられ、各遺伝子名はデータセット内の異なる位置で複数回出現します（つまり、GeneAは1000〜2500で、その後4000〜5000である可能性があります。そして、私はそれらの2つの遺伝子をマージしたくはなく、連続したものだけをマージします）。

編集2：私は以下のTimPソリューションを使用しました。マージの効率に問題が発生しました。また、テキストファイルをスタックオーバーフローに投稿して、データが実際にどのように見えるかを示してから、これまでのスクリプトを投稿する方法はありますか？

 # Let ValidMerge be a logical vector of MergeDown corresponding to the rows in the data
 # as defined by TimP below
 # The data.frame is called RMDB 
 # TeMerge function is previously defined to merge two rows Ri and Rj into one entry
 # which spans the start of Ri to the end of Rj (with same name and strand)

 COUNT = 0
 RMDB.OUT = 0

 while (COUNT < nrow(RMDB)) { # Cycle through all rows of RMDB
  COUNT = COUNT + 1

  # Is this position a merger start?
    # If yes, then returns position in startpt
    # which will be the end position in endpt
    # If no, returns 0

  Merge = match(COUNT,startpt,nomatch=0) 

  if ( Merge == 0 ){
    # No merge starts at this position
    RMDB.OUT = rbind(RMDB.OUT,RMDB[COUNT,]) # Append COUNT Row to output

  }else{
    # Merge COUNT row in RMDB to its endpoint 

      RMDB.OUT = rbind(RMDB.OUT,
                     TeMerge(RMDB[COUNT,],RMDB[endpt[Merge],]))

    #print(paste('Merging Row',COUNT,' and Row',endpt[Merge]))

    COUNT = endpt[Merge] # Move Count to the endpoint      
  }
}

このスクリプトは、私が興味を持っている結果を正確に示します。唯一の問題は、data.frameに5 000 000エントリがあり、GAPサイズのいくつかのパラメーターを使用して分析を実行して結果を比較したいことです。コードのこの部分をより効率的に書き直す方法はありますか？この時点までのすべては、妥当な時間（〜数分）で実行されます。のこの部分は、データの700000サブセットで3時間以上続いています。

編集3：すべてのケースが上部にあるスパイクデータ（MIR3）。列1：4、8、11：15を無視します。

dput（RMDB）構造（list（V1 = c（3612L、318L、318L、318L、318L、318L、318L、318L、318L、318L、318L、318L、318L、318L、741L、444L、741L、407L、2059L、 407L、656L、2058L、257L、4051L、456L、351L、850L、335L、1000L、1566L、236L、588L、3877L、750L、2292L、783L、747L、666L、260L、1118L、341L、7010L、320L、7010L 249L、458L、24L、631L、631L、875L）、V2 = c（11.4、23、23、23、23、23、23、23、23、23、23、23、23、23、18.7、11.6、24.9 、28.8、14.1、28.8、23.6、18.3、25、7、8.2、23.6、24.9、29.5、13.5、19.4、34.8、17.4、22.9、27.6、12、26.6、30.4、12.9、38.5、35.4、27.8、19.2 、0、17.3、21.2、19.3、0、3.9、26.6、22.6）、V3 = c（27、3.8、3.8、3.8、3.8、3.8、3.8、3.8、3.8、3.8、3.8、3.8、3.8、3.8、 4.5、0、0、7.3、0.3、7.3、12.2、4.9、0、0.4、0、1.8、15.1、12.7、0、3.6、3、7.5、14、9.4、0、14.1、4.1、4、1.4、 9.4、5.9、2.7、0、9.3、1、0、0、0、1.8、9.5）、V4 = c（1.3、0、0、0、0、0、0、0、0、0、0、0、0、0、0、3.1、1.1 、3、0.3、3、3、0、0、0、0、3.6、0.8、2.7、0、2.8、0.4、0.8、1.6、6.4、0、3.8、3.7、0、0、2.6、1.5、2.5 、2.4、1.4、2、5、0、0、0.6、0.5）、V5 =構造（c（1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、 1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、 1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L）、. Label = "chr1"、class = "factor"）、V6 = c（10469L、20001L、20210L、21000L 、21201L、22000L、22201L、23000L、23201L、20000L、20001L、24001L、24205L、24405L、0L、30855L、30953L、31293L、31436L、31734L、32841L、33048L、33466L、33529L、34048L、34451L、34565L、35217 、37045L、37733L、38059L、38256L、39465L、39624L、39925L、40333L、40629L、40736L、41380L、42370L、43243L、44836L、44877L、45887L、46079L、46217L、46416L、46553L、46893L）、V7 = c（11447L、20200L、20400L、21200L、21400L、22200L、22400L、23200L、23400L、20001L、 、24200L、24400L、24600L、0L、30952L、31131L、31435L、31733L、31754L、33037L、33456L、33509L、34041L、34108L、34560L、34921L、35366L、35499L、37431L、37861L、38191L、39464L、39623L 、40626L、40729L、40878L、42285L、42504L、44835L、44876L、45753L、45987L、46198L、46240L、46493L、46722L、47092L）、V8 =構造（c（38L、37L、37L、37L、37L、37L、37L、 37L、37L、37L、37L、37L、37L、37L、36L、35L、34L、33L、32L、31L、30L、29L、28L、27L、26L、25L、24L、23L、22L、21L、20L、19L、 18L、17L、16L、15L、14L、13L、12L、11L、10L、9L、8L、7L、6L、5L、4L、3L、2L、1L）、.ラベル=c（ "（249203529）"、 "（249203899）"、 "（249204128）"、 "（249204381​​）"、 "（249204423）"、 "（249204634）"、 "（249204868）"、 "（249205745） "、"（249205786） "、"（249208117） "、"（249208336） "、"（249209743） "、"（249209892） "、"（249209995） "、"（249210327） "、"（249210697） "、 "（249210998）"、 "（249211157）"、 "（249212430）"、 "（249212760）"、 "（249213190）"、 "（249215122）"、 "（249215255）"、 "（249215700）"、 "（ 249216061） "、"（249216513） "、"（249216580） "、"（249217112） "、"（249217165） "、"（249217584） "、"（249218867） "、"（249218888） "、"（249219186） "、"（249219490） "、"（249219669） "、"（249219773） "、"（249235266） "、"（249239174） "）、class =" factor "）、V9 =構造（c（2L、1L、1L、2L、2L、2L、2L、2L、1L、2L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L 、2L、1L、2L、2L、2L、1L、1L、1L、1L、1L、1L、1L、1L、1L、2L、2L、2L、1L、1L、1L、1L、1L、1L、1L、2L 、1L、2L）、. Label = c（ "+"、 "C"）、class = "factor"）、V10 = structure（c（32L、23L、23L、23L、23L、23L、24L、23L、23L 、23L、23L、23L、23L、23L、34L、33L、27L、26L、1L、26L、22L、12L、5L、14L、13L、16L、30L、25L、2L、7L、16L、15L、29L、28L 、3L、28L、15L、4L、18L、8L、19L、10L、31L、10L、9L、11L、6L、17L、20L、21L）、. Label = c（ "AluJo"、「AluJr」、「AluSx」、「AluSz6」、「AluYc」、「AT_rich」、「Charlie5」、「ERVL-E-int」、「L1M5」、「L1MA8」、「L1MA9」、「L1MB5」、「L1P1 "、" L1PA6 "、" L2a "、" L2c "、" LTR12F "、" LTR16C "、" MamRep1527 "、" MER45A "、" MER58A "、" MIR "、" MIR3 "、" MIR3a "、" MIRb "、 "MIRc"、 "MLT1A"、 "MLT1E1A"、 "MLT1E1A-int"、 "MLT1J2"、 "（TAAA）n"、 "TAR1"、 "（TC）n"、 "XXXXX"）、class = "factor" ）、V11 = structure（c（10L、13L、13L、13L、13L、13L、13L、13L、13L、13L、13L、13L、13L、13L、14L、11L、9L、13L、12L、13L、13L、 3L、12L、3L、3L、4L、9L、13L、12L、1L、4L、4L、9L、9L、12L、9L、4L、12L、8L、8L、6L、3L、11L、3L、3L、3L、5L、 7L、2L、1L）、. Label = c（ "DNA / hAT-Charlie"、 "DNA / hAT-Tip100"、 "LINE / L1"、 "LINE / L2"、 "Low_complexity"、 "LTR"、 "LTR / ERV1 "、" LTR / ERVL "、" LTR / ERVL-MaLR "、" Satellite / telo "、" Simple_repeat "、" SINE / Alu "、" SINE / MIR "、" XXXXXXXX "）、class =" factor " ）、V12 = structure（c（19L、5L、5L、5L、5L、5L、5L、5L、5L、5L、5L、5L、5L、5L、2L、9L、8L、25L、2L、10L、9L、 23L、2L、1L、15L、12L、1L、6L、2L、11L、16L、13L、7L、2L、2L、8L、14L、1L、3L、26L、17L、18L、9L、21L、22L、24L、 2L、4L、27L、20L）、. Label = c（ "（0）"、 "1"、 "（113） "、"（115） "、"（119） "、"（13） "、" 131 "、" 173 "、" 2 "、" 218 "、" 2234 "、"（231） "、" 2705 "、" 2923 "、" 2970 "、" 3242 "、" 359 "、" 3715 "、"（399） "、"（4） "、" 5306 "、" 5334 "、" 5746 "、" 6167 " 、"67"、 "（685）"、 "7"）、class = "factor"）、V13 = c（1712L、143L、143L、143L、143L、143L、143L、143L、143L、143L、143L、143L 、143L、143L、162L、96L、349L、217L、298L、238L、216L、6174L、44L、6154L、3030L、3156L、448L、255L、133L、2623L、3375L、2846L、1489L、172L、301L、666L、3217L 、312L、376L、4982L、499L、5305L、41L、6290L、5433L、6280L、24L、404L、178L、220L）、V14 =構造（c（23L、24L、24L、24L、24L、24L、24L、24L、24L、24L、24L、24L、24L、24L、8L、1L、10L、25L、4L、13L、21L、1L、 11L、26L、15L、14L、20L、30L、5L、2L、1L、27L、1L、18L、3L、4L、7L、6L、9L、19L、22L、29L、1L、2L、28L、16L、1L、 17L、1L、12L）、. Label = c（ "（0）"、 "（1）"、 "（11）"、 "（14）"、 "（179）"、 "208"、 "（209） "、"（212） "、" 232 "、"（25） "、"（255） "、" 3 "、"（30） "、" 3049 "、"（3116） "、"（32） "、 "327"、 "（388）"、 "4016"、 "41"、 "（46）"、 "（470）"、 "483"、 "49"、 "（51）"、 "5640"、 "（ 573） "、"（691） "、"（838） "、" 91 "）、class = "factor"）、V15 = c（2L、5L、5L、5L、5L、5L、5L、5L、5L、5L、5L、5L、5L、5L、22L、23L、22L、24L、25L、24L 、26L、27L、28L、29L、30L、31L、32L、33L、34L、35L、36L、37L、38L、38L、39L、38L、37L、40L、41L、42L、43L、44L、45L、44L、46L 、47L、48L、49L、50L、51L））、. Names = c（ "V1"、 "V2"、 "V3"、 "V4"、 "V5"、 "V6"、 "V7"、 "V8"、 "V9"、 "V10"、 "V11"、 "V12"、 "V13"、 "V14"、 "V15"）、class = "data.frame"、row.names = c（NA、-50L））40L、41L、42L、43L、44L、45L、44L、46L、47L、48L、49L、50L、51L））、. Names = c（ "V1"、 "V2"、 "V3"、 "V4"、 " V5 "、" V6 "、" V7 "、" V8 "、" V9 "、" V10 "、" V11 "、" V12 "、" V13 "、" V14 "、" V15 "）、class =" data.frame "、row.names = c（NA、-50L））40L、41L、42L、43L、44L、45L、44L、46L、47L、48L、49L、50L、51L））、. Names = c（ "V1"、 "V2"、 "V3"、 "V4"、 " V5 "、" V6 "、" V7 "、" V8 "、" V9 "、" V10 "、" V11 "、" V12 "、" V13 "、" V14 "、" V15 "）、class =" data.frame "、row.names = c（NA、-50L））

そして、これがValidMergeが適用された後の出力です。

dput（RMDB.OUT）構造（list（V1 = c（3612L、NA、NA、318L、318L、318L、318L、318L、318L、NA、741L、444L、741L、407L、2059L、407L、656L、2058L、 257L、4051L、456L、351L、850L、335L、1000L、1566L、236L、588L、3877L、750L、2292L、783L、747L、666L、260L、1118L、341L、7010L、320L、7010L、249L、458L、24L 631L、631L、875L）、V2 = c（11.4、NA、NA、23、23、23、23、23、23、NA、18.7、11.6、24.9、28.8、14.1、28.8、23.6、18.3、25、7 、8.2、23.6、24.9、29.5、13.5、19.4、34.8、17.4、22.9、27.6、12、26.6、30.4、12.9、38.5、35.4、27.8、19.2、0、17.3、21.2、19.3、0、3.9、26.6 、22.6）、V3 = c（27、NA、NA、3.8、3.8、3.8、3.8、3.8、3.8、NA、4.5、0、0、7.3、0.3、7.3、12.2、4.9、0、0.4、0、 1.8、15.1、12.7、0、3.6、3、7.5、14、9.4、0、14.1、4.1、4、1.4、9.4、5.9、2.7、0、9.3、1、0、0、0、1.8、9.5） 、V4 = c（1.3、NA、NA、0、0、0、0、0、0、NA、0、3.1、1.1、3、0.3、3、3、0、0、0、0、3.6、0.8、2.7、0、2.8、0.4、0.8、1.6、6.4、0、3.8、3.7、0、0、 2.6、1.5、2.5、2.4、1.4、2、5、0、0、0.6、0.5）、V5 =構造（c（1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L 、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L 、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L）、. Label = "chr1"、class = "factor"）、V6 = c（10469L、20001L、21000L、22000L、22201L、 23000L、23201L、20000L、20001L、24001L、0L、30855L、30953L、31293L、31436L、31734L、32841L、33048L、33466L、33529L、34048L、34451L、34565L、35217L、35367L、37045L、37733L、38059L、38256 39624L、39925L、40333L、40629L、40736L、41380L、42370L、43243L、44836L、44877L、45887L、46079L、46217L、46416L、46553L、46893L）、V7 = c（11447L、20400L、21400L、22200L、22400L、23200L、23400L、20001L、20200L、24600L、0L、30952L、31131L、31435L、31733L、31754L、33037L、33456L、33509L、34041L、34108L、34560L、34921L、35366L、35499L、37431L、37861L、38191 39623L、39924L、40294L、40626L、40729L、40878L、42285L、42504L、44835L、44876L、45753L、45987L、46198L、46240L、46493L、46722L、47092L）、V8 =構造（c（38L、37L、37L、37L、37L） 、37L、37L、37L、37L、37L、36L、35L、34L、33L、32L、31L、30L、29L、28L、27L、26L、25L、24L、23L、22L、21L、20L、19L、18L、17L 、16L、15L、14L、13L、12L、11L、10L、9L、8L、7L、6L、5L、4L、3L、2L、1L）、. Label = c（ "（249203529）"、 "（249203899）" 、"（249204128）"、 "（249204381​​）"、 "（249204423）"、 "（249204634）"、 "（249204868）"、 "（249205745）"、 "（249205786）"、 "（249208117） "、"（249208336） "、"（249209743） "、"（249209892） "、"（249209995） "、"（249210327） "、"（249210697） "、"（249210998） "、"（249211157 ）"、"（249212430） "、"（249212760） "、"（249213190） "、"（249215122） "、"（249215255） "、"（249215700） "、"（249216061） "、"（249216513） " 、"（249216580）"、 "（249217112）"、 "（249217165）"、 "（249217584）"、 "（249218867）"、 "（249218888）"、 "（249219186）"、 "（249219490）"、 " （249219669） "、"（249219773） "、"（249235266） "、"（249239174） "）、class =" factor "）、V9 = structure（c（2L、1L、2L、2L、2L、2L、1L 、2L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、2L、1L、2L、2L、2L、1L、1L、1L、1L、1L、1L、1L、1L、1L、2L、2L、2L、1L、1L、1L、1L、1L、1L、 1L、2L、1L、2L）、. Label = c（ "+"、 "C"）、class = "factor"）、V10 = structure（c（32L、23L、23L、23L、24L、23L、23L、 23L、23L、23L、34L、33L、27L、26L、1L、26L、22L、12L、5L、14L、13L、16L、30L、25L、2L、7L、16L、15L、29L、28L、3L、28L、 15L、4L、18L、8L、19L、10L、31L、10L、9L、11L、6L、17L、20L、21L）、. Label = c（ "AluJo"、 "AluJr"、 "AluSx"、 "AluSz6"、 "AluYc"、 "AT_rich"、 "Charlie5"、 "ERVL-E-int"、 "L1M5"、 "L1MA8"、 "L1MA9"、 "L1MB5"、 "L1P1"、 "L1PA6"、 "L2a"、 "L2c "、" LTR12F "、" LTR16C "、" MamRep1527 "、"MER45A "、" MER58A "、" MIR "、" MIR3 "、" MIR3a "、" MIRb "、" MIRc "、" MLT1A "、" MLT1E1A "、" MLT1E1A-int "、" MLT1J2 "、"（TAAA）n "、" TAR1 "、"（TC）n "、" XXXXX "）、class =" factor "）、V11 = structure（c（10L、13L、13L、13L、13L、13L、13L、13L、13L、13L 、14L、11L、9L、13L、12L、13L、13L、3L、12L、3L、3L、4L、9L、13L、12L、1L、4L、4L、9L、9L、12L、9L、4L、12L、8L 、8L、6L、3L、11L、3L、3L、3L、5L、7L、2L、1L）、. Label = c（ "DNA / hAT-Charlie"、 "DNA / hAT-Tip100"、 "LINE / L1" 、"LINE / L2"、 "Low_complexity"、 "LTR"、 "LTR / ERV1"、 "LTR / ERVL"、 "LTR / ERVL-MaLR"、 "Satellite / telo"、 "Simple_repeat"、 "SINE / Alu "、" SINE / MIR "、" XXXXXXXX "）、class =" factor "）、V12 = structure（c（19L、NA、NA、5L、5L、5L、5L、5L、5L、NA、2L 、9L、8L、25L、2L、10L、9L、23L、2L、1L、15L、12L、1L、6L、2L、11L、16L、13L、7L、2L、2L、8L、14L、1L、3L、26L 、17L、18L、9L、21L、22L、24L、2L、4L、27L、20L）、. Label = c（ "（0）"、 "1"、 "（113）"、 "（115）"、 " （119） "、"（13） "、" 131 "、" 173 "、" 2 "、" 218 "、" 2234 "、"（231） "、" 2705 "、" 2923 "、" 2970 "、" 3242 "、" 359 "、" 3715 "、"（399） "、"（4） "、" 5306 "、" 5334 "、" 5746 "、" 6167 "、" 67 "、"（685） "、" 7 "）、class =" factor "）、V13 = c（1712L、NA、NA、143L、143L、143L、143L、143L、143L、NA、162L、96L、349L、217L、298L、238L、216L、6174L、44L、6154L、3030L、3156L、448L、 255L、133L、2623L、3375L、2846L、1489L、172L、301L、666L、3217L、312L、376L、4982L、499L、5305L、41L、6290L、5433L、6280L、24L、404L、178L、220L）、V14=構造（c（23L、NA、NA、24L、24L、24L、24L、24L、24L、NA、8L、1L、10L、25L、4L、13L、21L、1L、11L、26L、15L、14L、20L、30L 、5L、2L、1L、27L、1L、18L、3L、4L、7L、6L、9L、19L、22L、29L、1L、2L、28L、16L、1L、17L、1L、12L）、. Label = c （"（0）"、 "（1）"、 "（11）"、 "（14）"、 "（179）"、 "208"、 "（209）"、 "（212）"、 "232" 、"（25）"、 "（255）"、 "3"、 "（30）"、 "3049"、 "（3116）"、 "（32）"、 "327 "、"（388） "、" 4016 "、" 41 "、"（46） "、"（470） "、" 483 "、" 49 "、"（51） "、" 5640 "、"（573 ）"、"（691） "、"（838） "、" 91 "）、class =" factor "）、V15 = c（2L、5L、5L、5L、5L、5L、5L、5L、5L、5L 、22L、23L、22L、24L、25L、24L、26L、27L、28L、29L、30L、31L、32L、33L、34L、35L、36L、37L、38L、38L、39L、38L、37L、40L、41L 、42L、43L、44L、45L、44L、46L、47L、48L、49L、50L、51L））、. Names = c（ "V1"、 "V2"、 "V3"、 "V4"、 "V5"、 "V6"、 "V7"、 "V8"、 "V9"、 "V10"、 "V11"、 "V12"、 "V13"、 "V14"、 "V15"）、class = "data.frame"、行.names = c（NA、-46L））

編集4：申し訳ありませんが、ここに簡略化されたバージョンがあります：Initial Data.frame

dput（RMDB.cut）

structure（list（Chromosome = structure（c（1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L） 、.Label = "chr1"、class = "factor"）、Start = c（10469L、20001L、20210L、21000L、21201L、22000L、22201L、23000L、23201L、20000L、20001L、24001L、24205L、24405L、0L、30855L 、30953L、31293L、31436L、31734L）、エンド= c（11447L、20200L、20400L、21200L、21400L、22200L、22400L、23200L、23400L、20001L、20200L、24200L、24400L、24600L、0L、30952L、31131L、31435L、 31733L、31754L）、(Left)=構造（c（38L、37L、37L、37L、37L、37L、37L、37L、37L、37L、37L、37L、37L、37L、36L、35L、34L、33L、32L、31L）、. Label = c （"（249203529）"、 "（249203899）"、 "（249204128）"、 "（249204381​​）"、 "（249204423）"、 "（249204634）"、 "（249204868）"、 "（249205745）"、 " （249205786） "、"（249208117） "、"（249208336） "、"（249209743） "、"（249209892） "、"（249209995） "、"（249210327） "、"（249210697） "、"（249210998 ）"、"（249211157） "、"（249212430） "、"（249212760） "、"（249213190） "、"（249215122） "、"（249215255） "、"（249215700） "、"（249216061） " 、"（249216513）"、 "（249216580） "、"（249217112） "、"（249217165） "、"（249217584） "、"（249218867） "、"（249218888） "、"（249219186） "、"（249219490） "、"（249219669 ）"、"（249219773） "、"（249235266） "、"（249239174） "）、class =" factor "）、Strand = structure（c（2L、1L、1L、2L、2L、2L、2L、2L 、1L、2L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L）、. Label = c（ "+"、 "C"）、class = "factor"）、repName = structure （c（32L、23L、23L、23L、23L、23L、24L、23L、23L、23L、23L、23L、23L、23L、34L、33L、27L、26L、1L、26L）、. Label = c（ " AluJo "、" AluJr "、" AluSx "、" AluSz6 "、" AluYc "、" AT_rich "、" Charlie5 "、" ERVL-E-int "、「L1M5」、「L1MA8」、「L1MA9」、「L1MB5」、「L1P1」、「L1PA6」、「L2a」、「L2c」、「LTR12F」、「LTR16C」、「MamRep1527」、「MER45A」、「MER58A」 "、" MIR "、" MIR3 "、" MIR3a "、" MIRb "、" MIRc "、" MLT1A "、" MLT1E1A "、" MLT1E1A-int "、" MLT1J2 "、"（TAAA）n "、" TAR1 、"（TC）n"、 "XXXXX"）、class = "factor"）、ValidMerge = c（FALSE、TRUE、FALSE、TRUE、FALSE、FALSE、FALSE、FALSE、FALSE、FALSE、FALSE、TRUE、TRUE、 FALSE、FALSE、FALSE、FALSE、FALSE、FALSE、FALSE））、. Names = c（ "Chromosome"、 "Start"、 "End"、 "（Left）"、 "Strand"、 "repName "、" ValidMerge "）、row.names = c（NA、20L）、class =" data.frame "）

そして合併後のアウトプット

dput（RMDB.out.cut）

structure（list（Chromosome = structure（c（1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L、1L）、.Label = "chr1"、 class = "factor"）、Start = c（ "10469"、 "20001"、 "21000"、 "22000"、 "22201"、 "23000"、 "23201"、 "20000"、 "20001"、 "24001" 、"0"、 "30855"、 "30953"、 "31293"、 "31436"、 "31734"）、End = c（ "11447"、 "20400"、 "21400"、 "22200"、 "22400"、 「23200」、「23400」、「20001」、「20200」、「24600」、「0」、「30952」、「31131」、「31435」、「31733」、「31754」）、(Left)= structure（c（38L、37L、37L、37L、37L、37L、37L、37L、37L、37L、36L、35L、34L、33L、32L、31L）、. Label = c（ "（249203529）"、 " （249203899） "、"（249204128） "、"（249204381​​） "、"（249204423） "、"（249204634） "、"（249204868） "、"（249205745） "、"（249205786） "、"（249208117 ）"、"（249208336） "、"（249209743） "、"（249209892） "、"（249209995） "、"（249210327） "、"（249210697） "、"（249210998） "、"（249211157） " 、"（249212430）"、 "（249212760）"、 "（249213190）"、 "（249215122）"、 "（249215255）"、 "（249215700）"、 "（249216061）"、 "（249216513）"、 " （249216580） "、"（249217112） "、"（249217165） "、"（249217584） "、"（249218867） "、"（249218888） "、"（249219186） "、"（249219490） "、"（249219669） "、"（249219773 ）"、"（249235266） "、"（249239174） "）、class =" factor "）、Strand = structure（c（2L、1L、2L、2L、2L、2L、1L、2L、1L、1L、1L 、1L、1L、1L、1L、1L）、. Label = c（ "+"、 "C"）、class = "factor"）、repName = structure（c（32L、23L、23L、23L、24L、23L 、23L、23L、23L、23L、34L、33L、27L、26L、1L、26L）、. Label = c（ "AluJo"、 "AluJr"、 "AluSx"、 "AluSz6"、 "AluYc"、 "AT_rich" 、"Charlie5"、 "ERVL-E-int"、 "L1M5"、 "L1MA8"、 "L1MA9"、 "L1MB5」、「L1P1」、「L1PA6」、「L2a」、「L2c」、「LTR12F」、「LTR16C」、「MamRep1527」、「MER45A」、「MER58A」、「MIR」、「MIR3」、「MIR3a」 、"MIRb"、 "MIRc"、 "MLT1A"、 "MLT1E1A"、 "MLT1E1A-int"、 "MLT1J2"、 "（TAAA）n"、 "TAR1"、 "（TC）n"、 "XXXXX"）、 class = "factor"）、ValidMerge = c（FALSE、TRUE、FALSE、TRUE、FALSE、FALSE、FALSE、FALSE、FALSE、FALSE、FALSE、TRUE、TRUE、FALSE、FALSE、FALSE））、. Names = c（ "Chromosome"、 "Start"、 "End"、 "（Left）"、 "Strand"、 "repName"、 "ValidMerge"）、row.names = c（ "2"、 "3"、 "4"、「6」、「7」、「8」、「9」、「10」、「11」、「111」、「15」、「16」、「17」、「18」、「19」、「20」 "）、class =" data.frame "）

4

4 に答える 4

1

戦略はDoMergeと呼ばれる別の列を生成することだと思います-各行R_i（iの範囲は1からn-1）について、名前とストランドがR_iとR_ {i + 1}の間で一致する場合、DoMergeはTRUEであり、End for R_iは、R_ {i + 1}の開始に十分に近いです（正しい値の場合は100以内）。行nのDoMergeは、慣例によりFALSEです。直感的には、DoMergeがTRUEであるということは、その行を次の行とマージすることが有効であることを意味します。

次に、TRUEの連続した文字列があるすべての行をマージします。それが最善の戦略であることに同意すれば、そのための簡単なコードをいくつかノックアップできます。:)

アップデート：

mydfがName、Strand、Start、Endの各列を持つ情報のデータフレームであると仮定した場合のタスクのコードは次のとおりです...以下の出力は、マージする必要のある開始点と終了点です。マージする必要がありますそれは簡単なはずです:)

distanceThresh = 100
isSameName=(mydf$Name==c(mydf$Name[-1],"void"))
isSameStrand=(mydf$Strand==c(mydf$Strand[-1],"void"))
isWithinDistance=(c(mydf$Start[-1],max(mydf$End)+2*distanceThresh)
                  -mydf$End) <= distanceThresh
validMerge = isSameName & isSameStrand & isWithinDistance

fthent=which(!validMerge & c(validMerge[-1],FALSE))
tthenf=which(validMerge & !c(validMerge[-1],TRUE))
startpt = fthent+1; if (validMerge[1]) {startpt=c(1,startpt)}
endpt = tthenf+1

instructions=NULL
for (kk in seq_along(startpt)) {
instructions = c(instructions,
               paste("Merge",startpt[kk],"to",endpt[kk],"inclusive",sep=" "))
}

これがすべて理にかなっているかどうか教えてください！:)

マージ方法（6月8日）：

このようなものはどうですか（いくつかのテストはありましたが、実際のデータではありません）...

doMerge = unlist(mapply(function(x,y) {seq(x,y,1)},startpt,endpt))
doNotMerge = setdiff(seq_along(validMerge),doMerge)
dataMerged=data.frame(Name=RMDB$Name[startpt], Strand=RMDB$Strand[startpt],
                      Start=RMDB$Start[startpt], End=RMDB$End[endpt],
                      stringsAsFactors=FALSE)
dataUnmerged=RMDB[doNotMerge,]

基本的に、マージが必要な行を示します（それぞれから対応するまでdoMerge実行されるシーケンスを設定するだけです）。そして、他のすべての行です（の長さがデータの長さであると仮定します）。startptendptdoNotMergevalidMerge

次にdataMerged、マージされたデータがどのように見えるかを直接構築します-明らかに、Name行から継承され、行から継承されます。（他に関心のある列がある場合は、それらがどこから来ているかを決定する必要があります。もちろん...）の行数は、もちろん、との長さに一致します。最後に、マージに不適格だったすべての行です。StrandStartstartptEndendptdataMergedstartptendptdataUnmerged

うまくいけば、上記のすべてが理にかなっており、すべてを元のシーケンスに戻すために組み合わせdataMergedてdataUnmerged並べ替えると（おそらく、これに使用できるインデックス列があります）、望ましい結果が得られることは明らかです。

そして、私は上記が非常に、非常に速く実行されることを期待しています:)

于 2012-06-07T07:45:57.237 に答える
1

これがGenomicRangesソリューションです。最初のステップは、1回限りで、パッケージとその依存関係をインストールすることです。

source("http://bioconductor.org/biocLite.R")
biocLite("GenomicRanges")

次に、パッケージをロードしGRanges、データからインスタンスを作成します。これは、というデータフレームに配置しました。df

library(GenomicRanges)
gr = with(df, GRanges(Chromosome, IRanges(Start, End), Strand, repName=repName))

データは実際にはもう少し複雑です。「GRangesのリスト」では、各リスト要素が遺伝子によって定義されているため、

grl = split(gr, values(gr)$repName)

reduceこれは要素ごとに行い、隣接する要素間の最小ギャップ幅が100の場合に縮小できるようにします。

merged = reduce(grl, min.gapwidth=100L)

これを強制的にadata.frameに戻すことができas(merged, "data.frame")ます。強制前の結果は次のようになります

> merged
GRangesList of length 8:
$AluJo
GRanges with 1 range and 0 elementMetadata cols:
      seqnames         ranges strand
         <Rle>      <IRanges>  <Rle>
  [1]     chr1 [31436, 31733]      +

$MIR3
GRanges with 7 ranges and 0 elementMetadata cols:
      seqnames         ranges strand
  [1]     chr1 [20001, 20400]      +
  [2]     chr1 [23201, 23400]      +
  [3]     chr1 [24001, 24600]      +
  [4]     chr1 [20000, 20001]      -
  [5]     chr1 [21000, 21400]      -
  [6]     chr1 [22000, 22200]      -
  [7]     chr1 [23000, 23200]      -

ととしてdata.frame

> as(merged, "data.frame")
   element seqnames start   end width strand
1    AluJo     chr1 31436 31733   298      +
2     MIR3     chr1 20001 20400   400      +
3     MIR3     chr1 23201 23400   200      +
4     MIR3     chr1 24001 24600   600      +
5     MIR3     chr1 20000 20001     2      -
6     MIR3     chr1 21000 21400   401      -
7     MIR3     chr1 22000 22200   201      -
8     MIR3     chr1 23000 23200   201      -

100,000の遺伝子に配置された100万行の場合、

> length(grl)
[1] 100000
> table(elementLengths(grl))
    10
100000
> system.time(reduce(grl, min.gapwidth=100))
   user  system elapsed
  9.468   0.064   9.553
于 2012-06-08T03:38:09.750 に答える
0

まず、データフレーム（）からさまざまな列を再ベクトル化しますdf。以下のように必要な操作を行ってから、data-frame

name_strand=paste(df$Name,df$Strand,sep=",")  #for grouping
st=df$Start
en=df$End
unq=unique(name_strand) #get the unique Name+Strand tags
ls1=list()
for(i in 1:length(unq)){
  index=which(name_strand %in% unq[i]) #get the index ids or group by unique Name+Strand 
  un_ls=unlist(strsplit(unq[i],split=","))
  ls1[[i]]=list(Name=un_ls[1],Strand=un_ls[2],Start=min(st[index]),End=max(en[index]))
}
df=as.data.frame(do.call(rbind, ls1))

Nameしたがって、これにより、が他の場所でStrand発生している場合の状況も解決されdata-frameます。

編集：

あなたの質問には疑問があるので、100未満にしたい場合や、行が順序付けられていない可能性があることgap_distanceを考慮して、解決策を変更しました。data-frameしたがって、私のソリューションでは、これらの両方を考慮に入れています（行がすでにとの降順で並べられている場合は、変更できStartますEnd：

for(i in 1:length(unq)){
  index=which(name_strand %in% unq[i]) #get the index ids or group by unique Name+Strand 
  un_ls=unlist(strsplit(unq[i],split=","))
  #previous sol: ls1[[i]]=list(Name=un_ls[1],Strand=un_ls[2],Start=min(st[index]),End=max(en[index]))

  new_st=st[index] #new Start vector from the "st" vector for this group by name_strand
  new_en=en[index]
  new_st=new_st[order(new_st)] #when data-frame rows are not ordered:
  new_en=new_en[order(new_en)]  #"order" method gives the index in ascending order of START

  #using embed method to lag and taking the diff for checking with gap distance considerations
  gap_diff=embed(new_st,2)[,1]-new_en[-length(new_en)]  #subtract from End vector except the last element of End  
  a_ind_st=1 #index of new_st (Start) vector
  b_ind_en=1 #index of new_en (End) vector     
  ls_ind=1  #index for list
  for (j in 1:(length(new_en)-1)){  #also acts as the index for gap_diff
      if (gap_diff[j] < gap_distance){ #where gap_distance = 100 (consideration #2) < or <=
          b_ind_en=j+1
           if (j + 1 > length(new_en)-1){  #special case for last element in vector
              ls1[[ls_ind]]=list(Name=un_ls[1],Strand=un_ls[2],Start=new_st[a_ind_st],End=new_en[b_ind_en]))
              ls_ind=ls_ind+1 #increment list index
           }
       }else{
          ls1[[ls_ind]]=list(Name=un_ls[1],Strand=un_ls[2],Start=new_st[a_ind_st],End=new_en[b_ind_en]))
          ls_ind=ls_ind+1 #increment list index
          a_ind_st=b_ind_en+1
          b_ind_en=b_ind_en+1
      }      
   }
}
df=as.data.frame(do.call(rbind, ls1))

この変更では、すべてが考慮されます。制限が必要ない場合は、ソリューションを変更できます。

于 2012-06-07T10:54:12.213 に答える
0

を使用してこれを行うことも可能です。これdplyrは、他の人がより簡単に理解でき、ループを使用しません。

  1. Name同じでStrand使用しているグループに行を分割しますgroup_by
  2. 各潜在的なグループ内のどの行が最小ギャップ範囲内にあるかを使用して決定しますmutate
  3. group_by次に、とを使用して隣接する行を「マージ」しますsummarize

dput動作させることができなかったため、出力にタイプミスがあるはずなので、代わりに小さいサンプルデータセットを使用しました。

library(dplyr)

X <- tibble(
    id = c("R1", "R2", "R3", "R4"),
    Name = c("GeneA", "GeneA", "GeneA", "GeneB"),
    Strand = c("+", "+", "+", "-"),
    Start = c(1000, 1510, 2001, 3100),
    End = c(1500, 2000, 2500, 4000)
)
print(X)

投稿からのサンプルデータ

# A tibble: 4 x 5
  id    Name  Strand Start   End
  <chr> <chr> <chr>  <dbl> <dbl>
1 R1    GeneA +       1000  1500
2 R2    GeneA +       1510  2000
3 R3    GeneA +       2001  2500
4 R4    GeneB -       3100  4000

dplyr解決：

# Set max distance between strands
max_gap <- 100

X %>%

    # consider row groups with the same Name and Strand
    group_by(Name, Strand) %>%

    # logically identify rows that meet the closeness criteria 
    mutate(group_id = (!is.na(lag(End)) & (lag(End) - Start) < max_gap)) %>%
    

    # Sum over the logical identifier to create a group identifier on every row
    ungroup() %>%
    mutate(group_id = cumsum(!group_id)) %>%

    # Merge rows that are in the same group
    group_by(group_id, Name, Strand) %>%
    summarize(id = min(id), Start = min(Start), End = max(End)) %>%
    ungroup() %>%

    # Keep only the desired columns
    select(id, Name, Strand, Start, End)

結果

# A tibble: 2 x 5
  id    Name  Strand Start   End
  <chr> <chr> <chr>  <dbl> <dbl>
1 R1    GeneA +       1000  2500
2 R4    GeneB -       3100  4000
于 2021-04-13T21:50:08.807 に答える