随机生成基因组位置，限制每个基因和样本的突变数量

Question

假设我有一个包含突变的基因组位置列表。每个突变都与一个基因和一个样本相关联。我的数据集如下所示：

chrom    pos   gene sample
chr1    1000    ABC S1
chr1    1500    ABC S2
chr1    1000    ABC S3
chr2    5000    XYZ S1
chr2    5000    XYZ S2
chr2    6000    XYZ S1
chr3    500     MNO S1

我的目标是生成一个类似的模拟突变表，其中每个样本的突变数量（跨所有基因）和每个基因的突变数量（跨所有样本）与参考突变表（上面那个）相同.在这种情况下：

Gene :
ABC : 3
XYZ : 3
MNO : 1

Sample :
S1 : 4
S2 : 2
S3 : 1

除此之外，我还有一个基因表：

gene    chrom   start   end
ABC      chr1   500     1100
ABC      chr1   1300    1600
ABC      chr1   2000    2500
XYZ      chr2   4000    5500
XYZ      chr2   5800    6500
MNO      chr3   200     300
MNO      chr3   400     600
MNO      chr3   800     1000

我们的想法是仅在这些间隔中选择位置以生成模拟突变表。变异表的大小 ~50K ；可生成~200K

模拟变异表示例：

chrom   pos    gene sample
chr1    600     ABC S1
chr1    1400    ABC S1
chr1    1500    ABC S2
chr2    4500    XYZ S1
chr2    6200    XYZ S1
chr2    6400    XYZ S2
chr3    900     MNO S3

您观察到每个基因和样本的突变数量与参考突变表中的相同。

我的第一个想法是首先使用基因表选择基因中的 X_i 个随机位置；其中 X_i = 基因 i 参考突变表的突变数。然后根据参考突变表中突变样本的数量为这些位置中的每一个分配一个样本。

在 R 中：

res <- 
    refmut %>% 
    group_by(gene) %>% 
    summarise(nmut=n()) %>% # compte number of mutations per gene
    right_join(gene.table) %>% # right join with gene table
    mutate(size = end-start + 1) %>% # compute size of each gene interval
    group_by(gene) %>% 
    sample_n(size=nmut,replace = T,weight = size) %>% # Randomly sample rows, proportional to the length of each range
    rowwise() %>% # for each row
    mutate(pos=sample(start:end,size=1)) %>% # Randomly sample uniformly within each chosen range
    ungroup() %>% # globally
    mutate(sample=sample(refmut$sample)) %>% # permute samples across positions
    select(-nmut,-start,-end,-size,chrom,pos,gene,sample) # format result 

在我的代码中，可能会在模拟之前计算一些行，例如间隔大小和 right_join ；走得更快。

还有什么想法吗？

可重现的数据集：

structure(list(chrom = c("chr1", "chr1", "chr1", "chr2", "chr2", 
"chr2", "chr3"), pos = c(1000L, 1500L, 1000L, 5000L, 5000L, 6000L, 
500L), gene = c("ABC", "ABC", "ABC", "XYZ", "XYZ", "XYZ", "MNO"
), sample = c("S1", "S2", "S3", "S1", "S2", "S1", "S1")), class = "data.frame", row.names = c(NA, 
-7L))

structure(list(gene = c("ABC", "ABC", "ABC", "XYZ", "XYZ", "MNO", 
"MNO", "MNO"), chrom = c("chr1", "chr1", "chr1", "chr2", "chr2", 
"chr3", "chr3", "chr3"), start = c(500L, 1300L, 2000L, 4000L, 
5800L, 200L, 400L, 800L), end = c(1100L, 1600L, 2500L, 5500L, 
6500L, 300L, 600L, 1000L)), class = "data.frame", row.names = c(NA, 
-8L))

Answer 1

我看到的主要内容是rowwise。您的原始表中将有 50,000 个组，每个组调用 sample。这是很多循环。

另一种方法是使用runif() 一次生成随机数并进行规范化。具体来说：

start+as.integer(runif(n()) * (size-1))

完整代码：

refmut %>% 
  count(gene, name = 'nmut') %>% #different - no faster
  right_join(gene.table)%>%
  mutate(size = end-start + 1) %>% 
  group_by(gene) %>% 
  sample_n(size=nmut,replace = T,weight = size)%>%
  ungroup()%>%
  mutate(pos = start + as.integer(runif(n()) * (size-1)), #different - should be faster
         sample = sample(refmut$sample))

# A tibble: 7 x 8
#  gene   nmut chrom start   end  size   pos sample
#  <chr> <int> <chr> <int> <int> <dbl> <int> <chr> 
#1 ABC       3 chr1   2000  2500   501  2176 S1    
#2 ABC       3 chr1    500  1100   601   966 S3    
#3 ABC       3 chr1    500  1100   601   807 S2    
#4 MNO       1 chr3    200   300   101   200 S2    
#5 XYZ       3 chr2   5800  6500   701  6368 S1    
#6 XYZ       3 chr2   5800  6500   701  5871 S1    
#7 XYZ       3 chr2   4000  5500  1501  5309 S1