从其他 data.table 应用一个 data.table 列函数

Question

我是 data.table 的新手，正在尝试切换。我有 2 个 data.tables（variable_sites 和 dt_bam）并想使用 variable_sites$POS（称之为 refPOS）使用来自 [=16] 的数据执行一个函数=].要在摘要 table 中获取变量 read_base，我想在 dt_bam 中找到一行 refPOS 小于 pos + qwidth 并从中提取一个字符字符串 dt_bam$seq 基于 refPOS 和 pos

之间的差异

我让它为 refPOS 的一个值工作，但我真的不知道如何在 data.table 语法中 sapply 一个 refPOS 的向量。感谢任何帮助。

这是我的代码：

dt_bam<-data.table(qname=lst[[1]],rname=lst[[2]],strand=lst[[3]],pos=lst[[4]],qwidth=lst[[5]],cigar=lst[[6]],
                   seq=as.character(lst[[7]]))
refPOS<-1000140 # renamed POS so not to confuse with pos
summ_tab <-  dt_bam[refPOS < pos +qwidth & refPOS >pos,
                    .(locus_pos=refPOS,read_base = substr(seq,abs(refPOS-pos),abs(refPOS-pos)))] 

# sapply(variable_sites[,POS],) then the individual values from variable_sites[POS] become refPOS

预期输出，如下所示，但 dt1 variable_sites[POS]:

中的每一行对应一行

    refPOS read_base
1: 1000140         C

这是一些示例数据：

> head(variable_sites)
    CHR     POS REF
1: chr1 1013855   G
2: chr1 1045080   G
3: chr1 1051873   C
4: chr1 1083795   C
5: chr1 1091327   C
6: chr1 1091421   T    

> head(dt_bam)
                qname rname strand     pos qwidth  cigar
1: SRR709972.27609810  chr1      + 1000135    101   101M
2: SRR709972.27609810  chr1      - 1000145    101   101M
3: SRR709972.23678227  chr1      + 1000545    101 91M10S
4: SRR709972.23678227  chr1      - 1000632    101   101M
5: SRR709972.11643848  chr1      + 1000651    101   101M
6: SRR709972.18299955  chr1      + 1000669    101   101M
                                                                                                     seq
1: GCCGCGGGGTGTGTGAACCCGGCTCCGCATTCTTTCCCACACTCGCCCCAGCCAATCGACGGCCGCGCTCCTCCCCCGCTCGCTGTCAGTCACGCCTCGGC
2: GTGTGAACCCGGCTCCGCATTCTTTCCCACACTCGCCCCAGCCAATCGACGGCCGCGCTCCTCCCCCGCTCGCTGTCAGTCACGCCTCGGCTCCGGGCGCG
3: CGAGCCTCGGTCTCGAGCCTCTTGGCTTCCTCCGCCCTTCCCCACTCCGGTCCCGGTTTGGGCCCTGCTCTGTCTCCGAGTTTGATCCGACCCCGCCTCGC
4: CGACACCGGCTCGGCCTCCGGGGGTCCCCCCCTCAGGTGTGCGGCCTGGAGCACGGAGGGCTGCAGAAAGCCTTGGGAGCGACAGAGCCGGGGGAAGGTTG
5: GGGGGTCCCACCCTCAGGTGTGCGGCCTGGAGCACGGAGGGCTGCAGAAAGCCTTGGGAGCGACAGAGCCGGGGGAAGGTTGGCGGCCGGGTCGGCAGGCG
6: TGTGCGGCCTGGAGCACGGAGGGCTGCAGAAAGCCTTGGGAGCGACAGAGCCGGGGGAAGGTTGGCTGCCGGGTCGGCAGGCGGGAGGGCGGAGTCAGCGG

> dput(head(variable_sites))
setDT(structure(list(CHR = c("chr1", "chr1", "chr1", "chr1", "chr1", 
"chr1"), POS = c(1013855L, 1045080L, 1051873L, 1083795L, 1091327L, 
1091421L), REF = c("G", "G", "C", "C", "C", "T")), row.names = c(NA, 
-6L), class = c("data.table", "data.frame")))

Answer 1

这就是您正在寻找的 data.table 方法。我们在 dt_bam 中创建一个临时变量 end ，然后执行非相等连接。请注意，在执行连接时，您必须使用 x.POS 来引用 variable_sites$POS。 POS 会给你错误的变量。 i.pos/pos/POS 都指的是 dt_bam$pos，因为默认情况下，您加入的变量（POS 在这种情况下）被第一个替换data.table中对应的变量（本例中为pos）加入.

library(data.table)

variable_sites[
  dt_bam[, end:=pos+qwidth], read_base:=substr(seq, x.POS - i.pos, x.POS - i.pos), 
  on = .(POS > pos, POS < end)
]
dt_bam[, end:=NULL]

输出

> variable_sites
    CHR     POS REF read_base
1: chr1 1013855   G      <NA>
2: chr1 1045080   G      <NA>
3: chr1 1051873   C      <NA>
4: chr1 1083795   C      <NA>
5: chr1 1091327   C      <NA>
6: chr1 1091421   T      <NA>
7: chr1 1000140   ?         C

数据

variable_sites <- data.table::setDT(structure(list(CHR = c("chr1", "chr1", "chr1", "chr1", "chr1", 
"chr1", "chr1"), POS = c(1013855L, 1045080L, 1051873L, 1083795L, 
1091327L, 1091421L, 1000140L), REF = c("G", "G", "C", "C", "C", 
"T", "?")), row.names = c(NA, -7L), class = c("data.table", "data.frame")))

dt_bam <- data.table::setDT(structure(list(qname = c("SRR709972.27609810", "SRR709972.27609810", 
"SRR709972.23678227", "SRR709972.23678227", "SRR709972.11643848", 
"SRR709972.18299955"), rname = c("chr1", "chr1", "chr1", "chr1", 
"chr1", "chr1"), strand = c("+", "-", "+", "-", "+", "+"), pos = c(1000135L, 
1000145L, 1000545L, 1000632L, 1000651L, 1000669L), qwidth = c(101L, 
101L, 101L, 101L, 101L, 101L), cigar = c("101M", "101M", "91M10S", 
"101M", "101M", "101M"), seq = c("GCCGCGGGGTGTGTGAACCCGGCTCCGCATTCTTTCCCACACTCGCCCCAGCCAATCGACGGCCGCGCTCCTCCCCCGCTCGCTGTCAGTCACGCCTCGGC", 
"GTGTGAACCCGGCTCCGCATTCTTTCCCACACTCGCCCCAGCCAATCGACGGCCGCGCTCCTCCCCCGCTCGCTGTCAGTCACGCCTCGGCTCCGGGCGCG", 
"CGAGCCTCGGTCTCGAGCCTCTTGGCTTCCTCCGCCCTTCCCCACTCCGGTCCCGGTTTGGGCCCTGCTCTGTCTCCGAGTTTGATCCGACCCCGCCTCGC", 
"CGACACCGGCTCGGCCTCCGGGGGTCCCCCCCTCAGGTGTGCGGCCTGGAGCACGGAGGGCTGCAGAAAGCCTTGGGAGCGACAGAGCCGGGGGAAGGTTG", 
"GGGGGTCCCACCCTCAGGTGTGCGGCCTGGAGCACGGAGGGCTGCAGAAAGCCTTGGGAGCGACAGAGCCGGGGGAAGGTTGGCGGCCGGGTCGGCAGGCG", 
"TGTGCGGCCTGGAGCACGGAGGGCTGCAGAAAGCCTTGGGAGCGACAGAGCCGGGGGAAGGTTGGCTGCCGGGTCGGCAGGCGGGAGGGCGGAGTCAGCGG"
)), row.names = c(NA, -6L), class = c("data.table", "data.frame")))

这些联接涉及隐式复制。

summ_tab <- variable_sites[dt_bam[, end:=pos+qwidth], .(
  refPOS = x.POS, 
  read_base = substr(seq, x.POS - i.pos, x.POS - i.pos)
), on = .(POS > pos, POS < end), nomatch=NULL]

或

summ_tab <- dt_bam[, end:=pos+qwidth][variable_sites, .(
  refPOS = i.POS, 
  read_base = substr(seq, i.POS - x.pos, i.POS - x.pos)
), on = .(pos < POS, end > POS), nomatch=NULL]

nomatch=NULL 删除在 variable_sites 中无法匹配的所有行。删除此开关然后上面的两种方法会给你不同的行为。选择你想要的那个。

根据我们的讨论对您的问题的扩展解决方案

dt_bam[, c("start", "end") := .(
  pos - qwidth * (strand == "-"), 
  pos + qwidth * (strand == "+")
)]

variable_sites[dt_bam, .(
  refPOS = x.POS, 
  read_base = substr(seq, x.POS - i.start, x.POS - i.start)
), on = .(POS > start, POS < end), nomatch=NULL]

dt_bam[, c("start", "end") := NULL]

从其他 data.table 应用一个 data.table 列函数

sapply over one data.table column function from other data.table

r

data.table