【问题标题】:Custom Merge Function in RR中的自定义合并函数
【发布时间】:2012-06-11 05:17:23
【问题描述】:

我有一个大型数据集,我想编写一个自定义合并函数以与 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 )

数据是一个 data.frame 行 R1:R4

到目前为止,我可以得到一个比较 Ri 和 Rj (j = i +1) 的函数,如果 Name 相同,Strand 相同,并且它们之间的差距小于 100,则合并它们。

GAP = Ri[End] - Rj[Start]

如果我将函数应用于每一行并建立输出。然后应该创建函数输出

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

我可以获得一个可以使用 apply 将两个连续元素合并为一个的工作函数,但我不知道如何合并三个连续元素。一个丑陋的解决方案是运行两个连续的元素合并功能,直到没有进行额外的合并,但这不是有效的。我有点想不通,任何见解都将不胜感激。

编辑:澄清一下,数据按连续染色体上的起始位置排序(未显示),每个基因名称在数据集中的不同位置多次出现(即 GeneA 可以是 1000-2500,然后再次出现在 4000 -5000,我不想合并这两个基因,只合并连续的)。

编辑 2:我在下面使用了 Tim P 解决方案。合并的效率出现了问题。还有一种方法可以将文本文件发布到堆栈溢出,以便我可以显示数据的真实样子,然后发布我的脚本吗?

 # 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 大小的几个参数来运行分析以比较结果。有没有办法可以重写这部分代码以提高效率?到目前为止,一切都在合理的时间内运行(大约几分钟)。这部分已经在 700 000 个数据子集上运行了 3 多个小时。

编辑 3: 所有案例都位于顶部的尖峰数据 (MIR3)。忽略列 1:4,8,11:15。

输入(RMDB) 结构(列表(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、35217L、35367L、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, 20200L, 24200L, 24400L、24600L、0L、30952L、31131L、31435L、31733L、31754L、 33037L、33456L、33509L、34041L、34108L、34560L、34921L、35366L、 35499L, 37431L, 37861L, 38191L, 39464L, 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, 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=“因子”),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=“因子”), V11 = 结构(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”、“卫星/telo”、“Simple_repeat”、 “SINE/Alu”、“SINE/MIR”、“XXXXXXXX”),class= “因子”),V12 = 结构(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= "系数"), 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=“因子”), 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))

这是应用 ValidMerge 后的输出。

输入(RMDB.OUT) 结构(列表(V1 = c(3612L,NA,NA,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, 不适用, 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、38256L、39465L、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、38191L、 39464L, 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=“因子”),V9 =结构(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= "因子"), V10 = 结构(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= "因子"), V11 = 结构(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-查理", “DNA/hAT-Tip100”、“LINE/L1”、“LINE/L2”、“Low_complexity”、“LTR”、 “LTR/ERV1”、“LTR/ERVL”、“LTR/ERVL-MaLR”、“卫星/telo”、“Simple_repeat”、 “SINE/Alu”、“SINE/MIR”、“XXXXXXXX”),class= “因子”),V12 = 结构(c(19L, 不适用,不适用,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= “因子”)、V13 = c(1712L、 北美,北美,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, 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= “因子”)、 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”,row.names = c(NA, -46L))

编辑 4:抱歉,这是简化版: 初始数据帧

输入(RMDB.cut)

结构(列表(染色体=结构(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"), 开始 = 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=“因子”),链=结构(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= "因子"), 代表名称=结构(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=“因子”), ValidMerge = c(假,真,假,真,假,假,假, 假,假,假,假,真,真,假,假,假, 错误,错误,错误,错误)),.Names = c(“染色体”,“开始”, "End", "(Left)", "Strand", "repName", "ValidMerge"), row.names = c(NA, 20L), class= "data.frame")

以及合并后的输出

输入(RMDB.out.cut)

结构(列表(染色体=结构(c(1L,1L,1L,1L,1L,1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = "chr1", class= "factor"), 开始 = c("10469", "20001", "21000", "22000", "22201", "23000", “23201”、“20000”、“20001”、“24001”、“0”、“30855”、“30953”、 "31293", "31436", "31734"), 结束 = c("11447", "20400", "21400", “22200”、“22400”、“23200”、“23400”、“20001”、“20200”、“24600”、 “0”、“30952”、“31131”、“31435”、“31733”、“31754”)、(Left)=结构(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= "因子"), 链 = 结构 (c(2L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("+", "C"), class= "因子"), repName = 结构(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= “因子”)、 ValidMerge = c(假,真,假,真,假,假,假, 假,假,假,假,真,真,假,假,假 )), .Names = c("染色体", "开始", "结束", "(左)", "链", "repName", "ValidMerge"), row.names = c("2", "3", "4", "6", "7", “8”、“9”、“10”、“11”、“111”、“15”、“16”、“17”、“18”、“19”、“20” ), class= "data.frame")

【问题讨论】:

  • 请贴出您目前使用的代码。否则我们必须重新发明它。
  • 数据是否已经按名称、链和开始排序?如果他们不是,他们不应该是吗?
  • 看起来很有趣!假设一切都正确排序(正如 BenBarnes 指出的那样),我们能否确定正确的想法基本上是尽可能地扩展合并的“链”,然后全部执行?所以如果我们发现 i 和 i+1 可以合并,我们看看 i+2 是否也可以合并,如果可以,我们继续直到链结束?就您的问题领域而言,这是否具有生物学意义?
  • @BenBarnes 不需要对数据进行排序。看看我的解决方案
  • Bioconductor 包 GenomicRanges 的价值使这个和许多其他范围内的操作变得微不足道(而且快速);请参阅 ?GRanges 和,例如 ?reduce,以及此和 IRanges 包中的小插曲(可从网页登录页面获得)。

标签: r merge bioinformatics


【解决方案1】:

我认为策略应该是生成另一个名为 DoMerge 的列 - 对于每一行 R_i(其中 i 的范围为 1 到 n-1),如果 Name 和 Strand 在 R_i 和 R_{i+1} 之间匹配,则 DoMerge 为 TRUE,并且 R_i 的 End 与 R_{i+1} 的 Start 足够接近(如果这是正确的值,则在 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 告诉您需要合并的行(只需设置一个从每个startpt 运行到相应endpt 的序列)。而doNotMerge 是所有其他行(假设validMerge 的长度是数据的长度)。

然后dataMerged 只是直接构造合并后的数据应该是什么样子——显然NameStrandStart 是从行startpt 继承的,End 是从行endpt 继承的。 (如果还有其他感兴趣的列,您必须确定它们来自哪里,显然......)dataMerged 中的行数当然与startptendpt 的长度相匹配。最后,dataUnmerged 是所有不符合合并条件的行。

希望以上所有内容都有意义,并且很明显,如果您将 dataMergeddataUnmerged 组合并重新排序以使所有内容恢复到原始顺序(大概有一个索引列可以用于此),那么您得到想要的结果。

而且我希望上面的内容确实会运行得非常非常快:)

【讨论】:

  • 这是解决问题的好方法。实施我的数据需要几个小时,但如果它在合理的时间内有效,我会通知您。
  • 感谢 Artem,让我知道这是怎么回事 - 如果您需要,很乐意帮助您完成最后的合并阶段(尽管我认为这很容易:合并后的结果显然具有相同的名称和链值作为您要合并的行,并且它从行 startpt 继承 Start 并从 row endpt 继承 End(如果有意义的话)。问:我可以仔细检查开始列下的数字是否严格增加? End 列也一样吗?并且无法想象开始/结束间隔可以重叠(例如,在您的示例中,R2 永远不能从 1490 而不是 1510 开始)?
  • Ri$End和Rj$Start之间的数字不增加的情况有两种;第一个是小重叠,第二个是另一个染色体开始,在这种情况下,差异将在 10^6 的数量级上。为了说明 isWithinDistance 的原因,我将其写为差异的绝对值。我目前遇到了合并问题,请参阅上面的编辑 2。
  • 随着您的合并方法更新,我认为以这种方式使用它更有意义。如果我理解正确,输出将是两个数据帧,所有合并元素在一帧中沿 startpt 排序,第二帧中来自 RMDB 的所有未合并元素。然后,我将新合并的元素附加到未合并的元素,并按唯一 ID 列重新排序它们(幸运地包含在 V15 中)。快速查看 ValidMerge 以获取完整数据表明有 20 000 个 Merge 行和大约 700 000 个未合并的行。我将在接下来的几个小时内实现它,并告诉你它是如何运行的。
  • 是的,您的总结是正确的。听起来不错——希望一切顺利。我想您会发现代码几乎立即运行。我想从我的角度来看,一个突出的问题是您建议对所有其他字段(列)做什么。如果它们在合并的“块”中采用相同的值,则没有问题,但如果不是,您需要决定正确的值以进行合并... :)
【解决方案2】:

这是GenomicRanges 解决方案。第一步,只有一次,是安装包及其依赖项

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)

您可以使用as(merged, "data.frame") 将其强制回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      -

对于一百万行,排列成 100000 个基因,我们有

> 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

【讨论】:

  • 我想我要去度过周末,详细学习一下生物导体。对于一般的第二代数据操作,您会推荐哪些软件包?特别是我正在处理 RNA-seq 数据。
  • GenomicRanges、GenomicFeatures、Biostrings、IRanges(所有这些都有小插曲,这是一个很好的信息来源)。还可以查看最近的 coursesRNASeq biocViews。
【解决方案3】:

我将首先重新矢量化数据框 (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))

因此,当您在 data-frame 的其他位置出现 NameStrand 时,这也可以解决这种情况。

编辑

由于您的问题存在疑问,我有一个修改后的解决方案,以防您想要 gap_distance data-frame 行可能未排序。所以我的解决方案同时考虑了这两个(如果行已经按StartEnd 的降序排列,您可以修改它:

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))

此修改将考虑所有因素。如果您不需要任何限制,您可以修改解决方案。

【讨论】:

  • minmax 方法对分组后获得起始和结束位置有很大帮助。
  • 另一个不错的解决方案。我在 cmets 中的问题是数据是否可以按 Name、Strand 和 Start 排序,或者是否需要将这些变量的不相交组分开。如果需要将不相交的组分开,这在原始帖子中并不清楚,那么您的解决方案将产生不正确的结果。
  • @BenBarnes,OP 从未提及任何需要注意的不相交的注意事项。从示例中 - 1500 和 1510 不是连续的,而是组合在一起的。但是他确实说间隙距离小于100。我没看清楚。需要为此考虑编辑我的解决方案。但是OP对此并不清楚,让他消除疑虑。 OP 说“到目前为止..”但不是“必需的”
【解决方案4】:

也可以使用dplyr 来做到这一点,这可能更容易被其他人理解并且仍然不使用循环。

  1. 使用group_by 将行分成具有相同NameStrand 的组
  2. 使用mutate 确定每个潜在组中的哪些行在最小间隙范围内
  3. 然后使用group_bysummarize“合并”相邻的行

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-05-13
    • 2019-04-05
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2018-04-16
    • 2020-03-15
    相关资源
    最近更新 更多