R："label" 一行基于另一个 data.table 的条件

Question

我有一个超过 100,000 行的 data.table (A)。有 3 列。

铬开始结束 1: 字符 1 6484847 6484896 2: chr1 6484896 6484945 3：chr1 6484945 6484994 4：chr1 6484994 6485043 5: chr1 6485043 6485092 ---<br> 183569: chrX 106893605 106893654 183570: chrX 106893654 106893703 183571: chrX 106893703 106893752 183572: chrX 106893752 106893801 183573: chrX 106893801 106894256

我想生成一个名为 "gene" 的新列，它为来自另一个 data.table 的每行注释提供一个标签，该注释具有约 90 行 (B)。如下所示：

chrom start end gene 1: chr1 6484847 6521004 ESPN 2: chr1 41249683 41306124 KCNQ4 3: chr1 55464616 55474465 BSND 42: chrX 82763268 82764775 POU3F4 43: chrX 100600643 100603957 TIMM8A 44: chrX 106871653 106894256 PRPS1

如果 data.table A 中的行起始值在 data.table B 的行起始值和结束值之内，我需要相应地用正确的基因标记 A 中的行。

例如，生成的完整 data.table A 将是

染色体起始结束基因 1: chr1 6484847 6484896 ESPN 2: chr1 6484896 6484945 ESPN 3: chr1 6484945 6484994 ESPN 4: chr1 6484994 6485043 ESPN 5: chr1 6485043 6485092 ESPN ---<br> 183569: chrX 106893605 106893654 TIMM8A 183570: chrX 106893654 106893703 TIMM8A 183571: chrX 106893703 106893752 TIMM8A 183572: chrX 106893752 106893801 TIMM8A 183573: chrX 106893801 106894256 TIMM8A

我尝试了一些嵌套循环来执行此操作，但这似乎会花费太长时间。我认为必须有一种方法可以使用 data.table 包来做到这一点，但我似乎无法弄清楚。

如有任何建议，我们将不胜感激。

Answer 1

虽然可以在 base R 中执行此操作（或可能使用 data.table），但我强烈建议使用 GenomicRanges；这是一个非常强大和灵活的 R/Bioconductor 库，专为此类任务而设计。

这是一个使用 GenomicRanges::findOverlaps 的例子：

# Sample data
df1 <- read.table(text =
    "chrom     start      end
     chr1   6484847   6484896
     chr1   6484896   6484945
     chr1   6484945   6484994
     chr1   6484994   6485043
     chr1   6485043   6485092", sep = "", header = T, stringsAsFactors = F);

df2 <- read.table(text =
    "chrom     start       end     gene
     chr1   6484847   6521004     ESPN
     chr1  41249683  41306124     KCNQ4
     chr1  55464616  55474465     BSND
     chrX  82763268  82764775     POU3F4
     chrX 100600643 100603957     TIMM8A
     chrX 106871653 106894256     PRPS1", sep = "", header = TRUE, stringsAsFactors = F);

# Convert to GRanges objects
gr1 <- with(df1, GRanges(chrom, IRanges(start = start, end = end)));
gr2 <- with(df2, GRanges(chrom, IRanges(start = start, end = end), gene = gene));

# Find features from gr1 that overlap with gr2
m <- findOverlaps(gr1, gr2);

# Add gene annotation as metadata to gr1
mcols(gr1)$gene[queryHits(m)] <- mcols(gr2)$gene[subjectHits(m)];
gr1;
#GRanges object with 5 ranges and 1 metadata column:
#      seqnames             ranges strand |        gene
#         <Rle>          <IRanges>  <Rle> | <character>
#  [1]     chr1 [6484847, 6484896]      * |        ESPN
#  [2]     chr1 [6484896, 6484945]      * |        ESPN
#  [3]     chr1 [6484945, 6484994]      * |        ESPN
#  [4]     chr1 [6484994, 6485043]      * |        ESPN
#  [5]     chr1 [6485043, 6485092]      * |        ESPN
#  -------
#  seqinfo: 1 sequence from an unspecified genome; no seqlengths

Answer 2

除了之外，还有另一种 data.table 方法使用 non-equi join 和 update on join.

A[B, on = .(chrom, start >= start, start <= end), gene := i.gene][]

    chrom     start       end  gene
 1:  chr1   6484847   6484896  ESPN
 2:  chr1   6484896   6484945  ESPN
 3:  chr1   6484945   6484994  ESPN
 4:  chr1   6484994   6485043  ESPN
 5:  chr1   6485043   6485092  ESPN
 6:  chrX 106893605 106893654 PRPS1
 7:  chrX 106893654 106893703 PRPS1
 8:  chrX 106893703 106893752 PRPS1
 9:  chrX 106893752 106893801 PRPS1
10:  chrX 106893801 106894256 PRPS1

根据 OP，A 和 B 已经是 data.table 个对象。因此，这种方法避免了对 GRanges 对象的强制转换。

可重现的数据

library(data.table)
A <- fread("rn         chrom     start      end 
     1:  chr1   6484847   6484896 
     2:  chr1   6484896   6484945 
     3:  chr1   6484945   6484994 
     4:  chr1   6484994   6485043 
     5:  chr1   6485043   6485092
183569:  chrX 106893605 106893654
183570:  chrX 106893654 106893703
183571:  chrX 106893703 106893752
183572:  chrX 106893752 106893801
183573:  chrX 106893801 106894256", drop = 1L)

B <- fread("rn    chrom     start       end     gene
 1:  chr1   6484847   6521004     ESPN
 2:  chr1  41249683  41306124     KCNQ4
 3:  chr1  55464616  55474465     BSND
42:  chrX  82763268  82764775     POU3F4
43:  chrX 100600643 100603957     TIMM8A
44:  chrX 106871653 106894256     PRPS1", drop = 1L)

R："label" 一行基于另一个 data.table 的条件

R: "label" a row based on conditions from another data.table

r

bioinformatics

data.table

可重现的数据