r - создать возможный исход генотипов потомства

Question

Пусть есть два бинарных вектора p1 и p2 длиной 20, которые представляют генотипы двух родительских особей.

p1 <- sample(c(0,1), 20, replace=T)
p2 <- sample(c(0,1), 20, replace=T)

Тогда первые 10 чисел представляют аллели на одной хромосоме иВторые 10 чисел представляют аллели на второй хромосоме, т.е. мы рассматриваем 10 двуаллельных локусов.Теперь я хотел бы создать все возможные генотипы, которые могут иметь потомство от обоих родителей.Как я могу создать эти последовательности в предположении рекомбинации (то есть, если один из аллелей в локусе i выбран из 2 аллелей p2, второй аллель в положении i + 10 должен быть взят из p1 и наоборот)

Answer 1

(Примечание: у меня была неправильная индексация в подсчете потомства в предыдущем проекте.)

Вот решение, основанное на моих комментариях выше.

library(compositions)

p1 <- matrix(sample(0:1, 20, replace = TRUE), ncol = 2)
p2 <- matrix(sample(0:1, 20, replace = TRUE), ncol = 2)

for (choice1 in 0:1023) {
  p1choices <- bit(choice1, 0:9) + 1
  for (choice2 in 0:1023) {
    p2choices <- bit(choice2, 0:9) + 1
    offspring <- cbind(p1[cbind(1:10, p1choices)], p2[cbind(1:10, p2choices)])
    # record this somehow
  }
}

Я пропустил шагзаписи всех генотипов потомства.Вы можете преобразовать столбцы offspring в два числа в 0: 1023, используя

apply(offspring, 2, function(x) sum(x*2^(0:9)))

, но вам решать, что с ними делать.

Отредактировано, чтобы добавить:

Вышеуказанные циклы генерируют около миллиона потомков, но во многих случаях это не нужно.Если p1 или p2 является гомозиготным (равные значения в обоих столбцах), не имеет значения, какой из них вы выберете.При использовании простой модели в среднем около половины локусов будут гомозиготными у каждого родителя, поэтому на самом деле требуется всего около тысячи вариантов.Эта версия кода учитывает это.Это сложнее (и, следовательно, более вероятно, содержит ошибки!), Но примерно в тысячу раз быстрее:

library(compositions)

p1 <- matrix(sample(0:1, 20, replace = TRUE), ncol = 2)
hetero1 <- p1[,1] != p1[,2]
count1 <- sum(hetero1)
p1choices <- rep(1, 10)
p2 <- matrix(sample(0:1, 20, replace = TRUE), ncol = 2)
hetero2 <- p2[,1] != p2[,2]
count2 <- sum(hetero2)
p2choices <- rep(1, 10)

for (choice1 in 0:(2^count1 - 1)) {
  p1choices[hetero1] <- bit(choice1, 0:(count1 - 1)) + 1
  for (choice2 in 0:(2^count2 - 1)) {
    p2choices[hetero2] <- bit(choice2, 0:(count2 - 1)) + 1
    offspring <- cbind(p1[cbind(1:10, p1choices)], p2[cbind(1:10, p2choices)]) 
    # record this somehow
  }
}