在R中计算单个值的数量

Question

我有一些RNA-seq数据，需要计算单体数量。我们定义单体为读取（read），其周围100个碱基内没有其他读取。

我有一个包括每个读取的起始坐标和终止坐标的数据框。我正在使用R进行操作。

目前我已编写了下面的代码，但是apply函数使用不正确，因此出现错误。

begin_end <- data.frame(begin_coordinate, final_coordinate)
apply(begin_end, 1, function(x) x[,1]-(x-1)[,2])

数据帧的前几行是：

> head(begin_end)

    begin   final
1   60507   60551
2   60790   60840
3   62004   62051
4   62819   62868
5   65141   65187

第一个似乎是单例，因为下一个读取开始的位置比它结束的位置要超过100个碱基，而且数据集的前几行都是这样。但是数据框很长，我希望不是所有的读取都是单例。

- newa123

Answer 1

以下是使用 data.table 实现与 @jeremycg 使用 dplyr 的 lag 和 lead 相同操作：

library(data.table)
setDT(begin_end)

begin_end[{
  d = begin - shift(final, type="lag")
  pmin(d, shift(d, type="lead"), na.rm=TRUE) > 100
}]

---

注释。 基本的data.table语法是DT[i,j]。其中i用于过滤输入数据，而j则用于修改输出结果。

我们在上面使用了i，但为了探究其工作原理，我们可以将相关向量放入j中：

begin_end[,{
  d       = begin - shift(final, type="lag")
  d_lead  = shift(d, type="lead")
  my_pmin = pmin(d, d_lead, na.rm=TRUE)
  c(.SD, list(d = d, d_lead = d_lead, my_pmin = my_pmin))
}]

#    begin final    d d_lead my_pmin
# 1: 60507 60551   NA    239     239
# 2: 60790 60840  239   1164     239
# 3: 62004 62051 1164    768     768
# 4: 62819 62868  768   2273     768
# 5: 65141 65187 2273     NA    2273

.SD 是一个已在表格中的列向量列表，是 Subset of Data（数据子集）的缩写。

- Frank

Answer 2

你似乎试图使用 (x-1) 从 apply 中获取先前的结束值。不幸的是，在 apply 函数族内部无法这样做。

幸运的是，有一个名为 lag 的函数（有几个，所以我将使用来自 dplyr 的函数）。这让我们可以按给定的条目数来 lag 一列：

begin_end$space <- begin_end$begin - dplyr::lag(begin_end$final)

这是输出结果：

  begin final space
1 60507 60551    NA
2 60790 60840   239
3 62004 62051  1164
4 62819 62868   768
5 65141 65187  2273

然后您可以尝试：

begin_end$issingle <- begin_end$space >= 100

- jeremycg

Answer 3

使用Bioconductor的GenomicRanges，我认为想法是从读取数据使用GenomicAlignments::readGAlignments()或makeGRangesFromDataFrame()创建一个GRanges()，然后使用resize()在每个方向上扩展它们，然后使用findOverlaps()来识别单个重叠自身的读取。大致如此。

library(GenomicRanges)
gr = GRanges(seqnames="chr1",
             IRanges(start=c(1000, 1150, 1500), width=100))
gr100 = resize(gr, width(gr) + 200, fix="center")
hits = findOverlaps(gr100)
gr100[tabulate(queryHits(hits), queryLength(hits)) == 1]

导致

>     gr100[tabulate(queryHits(hits), queryLength(hits)) == 1]
GRanges object with 1 range and 0 metadata columns:
      seqnames       ranges strand
         <Rle>    <IRanges>  <Rle>
  [1]     chr1 [1400, 1699]      *
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

这对于数百万条记录来说将会非常快速。

- Martin Morgan