2

我有一个基因行数据集，每个基因行都有它们的基因长度，我希望使用拒绝抽样通过基因长度分布从这些基因中取样 - 因为我在这个数据集中有太多基因太小而无法进入进一步分析（但是我不能自己设置一个截止点来删除它们）。我有一个基因长度的基因数据集可供采样，还有另一个基因长度的提议分布，我想使用它来对第一个数据集进行拒绝采样。

我的数据示例如下所示：

#df1 data to sample from:
Gene  Length
Gene1  5
Gene2  6
Gene3  400000
Gene4  1000
Gene5  25000
Gene6  10
Gene7  50
Gene8  4
Gene9  100
Gene10 2000

我的提案数据集：

#df2
Gene  Length
Gene1  5000
Gene2  60000
Gene3  400000
Gene4  1000
Gene5  25000
Gene6  10000
Gene7  50000
Gene8  4000
Gene9  1000
Gene10 2000

我没有任何统计背景，我正在尝试进行拒绝抽样（我的总体目标是获取长度极小基因较少的基因样本以进行进一步分析）。

要进行拒绝抽样，我正在从我在这里找到的教程中尝试这个：

X = df1$Length
U = df2$Length

accept = c()
count = 1

pi_x <- function(x) {
  new_x = (3/2)*(x^3)+(11/8)*(x^2)+(1/6)*(x)+(1/12)
  return(new_x)
}


while(count <= 50 & length(accept) < 50){
  test_u = U[count]
  test_x = pi_x(X[count])/(3.125*dunif(X[count],0,1))
  if (test_u <= test_x){
    accept = rbind(accept, X[count])
    count = count + 1
  }
  count = count + 1
}

我的问题是它只选择了 25 个基因（我进一步分析的理想采样范围是选择 50-100 个基因），而这 25 个基因中的大多数在采样后仍然很小。在运行此拒绝采样代码之前，我是否需要以X某种方式进行转换？我的实际数据df1是 800 个基因长度呈偏态/β 分布的基因（大多数都非常小）。还是我完全错过了我理解的其他东西？任何指导将不胜感激。

输入数据：

df1 <- structure(list(Gene = c("Gene1", "Gene2", "Gene3", "Gene4", "Gene5", 
"Gene6", "Gene7", "Gene8", "Gene9", "Gene10"), Length = c(5L, 
6L, 400000L, 1000L, 25000L, 10L, 50L, 4L, 100L, 2000L)), row.names = c(NA, 
-10L), class = c("data.table", "data.frame"))

df2 <- structure(list(Gene = c("Gene1", "Gene2", "Gene3", "Gene4", "Gene5", 
"Gene6", "Gene7", "Gene8", "Gene9", "Gene10"), Length = c(5000L, 
60000L, 400000L, 1000L, 25000L, 10000L, 50000L, 40000L, 1000L, 2000L)), row.names = c(NA, 
-10L), class = c("data.table", "data.frame"))

编辑：

我也试过：

sampled <- data.frame(proposal = df2$Length)
sampled$targetDensity <- dbeta(sampled$proposal, 3,6)

maxDens = max(sampled$targetDensity, na.rm = T)
sampled$accepted = ifelse(df1$Length < sampled$targetDensity / maxDens, TRUE, FALSE)

hist(sampled$proposal[sampled$accepted], freq = F, col = "grey", breaks = 100)
curve(dbeta(x, 3,6),0,1, add =T, col = "red")

但我确定我没有dbeta()正确使用，因为sampled$targetDensity输出全为零 - 有没有办法解决这个问题？我尝试过更改值，dbeta()但没有任何成功。

4

1 回答 1

0

如果您知道要采样的基因数量，则 sample 函数应该可以正常工作：

sampled = sample(df$genes, size = n, prob = df$length)

如果您想进一步降低对长度较小的基因进行采样的概率，您可以对prob参数的长度进行平方。参数 prob 将采样概率与每个元素相关联（此处基于长度）

如果你不知道你想要获取的基因数量，那么你可以定义自己的概率函数：

rejection.prob = function(x){
  if (x<too.small) {return(0)} # all genes smaller than too.small won't be sampled
  if (x > it.is.ok) {return(1)} # all these ones will be sampled
  if (x>too.small & (x<it.is.ok){
    # here make any function that is equal to 0 when x == too.small
    # and 1 when x == it.is.ok
    # it can be a simple linear function
}

请注意，您也可以将函数的输出rejection.prob用于sample函数。

根据您的期望，您可能希望拒绝功能更加连续（在 too.small 和 it.is.ok 处没有这些中断）。如果是这种情况，我会使用 sigmoid 函数，您可以在其中根据所需的输出调整参数。

于 2021-05-26T14:57:04.763 回答