Произвольно генерировать геномные позиции с ограничением на количество мутаций на ген и образец

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

Идея будет состоять в том, чтобы выбирать только позиции в этих интервалах для генерации имитированной мутацииТаблица.Размер таблицы мутаций ~ 50К;генерируемый ~ 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