2013-06-01 1 views
1

Я хочу определить точность вменения программы с использованием данных генотипа SNP, поэтому мне нужно замаскировать часть вызовов SNP для имитации отсутствующих данных.Имитация отсутствующих данных (т. Е. Данных маски) в R для проверки точности вменения

Я тестировал свой код на этом подмножестве данных маркера (см. Ниже). Имена столбцов - это имена отдельных лиц. Имена строк - это идентификаторы SNP. Набор данных содержит отсутствующие данные (отмечен как NA).

 SNPID  AR124 AR124 AR144 AR144 
[1,] "S10_28619" "G" "A" "A" "A" 
[2,] "S10_33499" "A" "A" "G" "G" 
[3,] "S10_47747" "T" "T" NA NA 

Я хочу, чтобы определить точность вменения, используя 10-кратный кросс проверки, так что мне нужно R, чтобы:

  1. маска 10% известных SNPs в общей сложности 10 различных времен (т.е. 10 раундов маскировки).

  2. Каждый раунд должен маскировать другой набор SNP.

  3. Каждый СНП должны быть замаскированы только один раз в течение этих 10 раундов (например, SNP S10_28619 будет отображаться как «NA» только один раз в течение 10 раундов маскировки).

Это код, который я использую:

##the function will return the marker matrix with an additional 10% missing data 

CV10NA= function (input, seed, foldno) {#input is the SNP matrix, seed is the random number seed, fold number is a number indicating which cross validation fold you are on 
set.seed(seed) 
a = unlist(input) 
b = is.na(a) #matrix b where TRUE indicates missing SNP and FALSE indicates known SNP 
pres = grep(FALSE, b) #finds cases of FALSE in matrix b and gives an integer vector containing the FALSEs' index numbers 
sets= sample(rep(1:10, c(length(pres)/10)), replace = FALSE) #repeat numbers 1 through 10 a total of length(pres)/10) times then randomly sample from values with no replacement 
a[which(sets==foldno)] = NA #find where sets==foldno in matrix a and replace it with NA 
a = matrix(a, ncol = ncol(input)) 
return(a) 
} 

функция, кажется, работает для foldno = 1 до 9, но не работает, когда foldno = 10. Не появляется сообщение об ошибке. ПРИМЕЧАНИЕ. Я удалил имена столбцов и строк перед выполнением функции, чтобы функция не рассматривала их как «маскируемые» элементы.

Вот выход для foldno = 1, 2, 3, и 10, соответственно:

> CV10NA(beagle.subset, 1, 1) 
    [,1] [,2] [,3] [,4] 
[1,] "G" "A" "A" NA 
[2,] "A" "A" "G" "G" 
[3,] "T" "T" NA NA 

> CV10NA(beagle.subset, 1, 2) 
    [,1] [,2] [,3] [,4] 
[1,] "G" "A" "A" "A" 
[2,] "A" NA "G" "G" 
[3,] "T" "T" NA NA 

> CV10NA(beagle.subset, 1, 3) 
    [,1] [,2] [,3] [,4] 
[1,] NA "A" "A" "A" 
[2,] "A" "A" "G" "G" 
[3,] "T" "T" NA NA 

> CV10NA(beagle.subset, 1, 10) 
    [,1] [,2] [,3] [,4] 
[1,] "G" "A" "A" "A" 
[2,] "A" "A" "G" "G" 
[3,] "T" "T" NA NA 

foldno = 10 не маскирует любую SNP в наборе данных.

Любые предложения/отзывы будут оценены! У меня нет опыта программирования, поэтому, пожалуйста, простите меня, если я сделаю очевидную ошибку или спрошу «глупый» вопрос.

Дополнительные попытки/мысли: Я попытался отладить код, запустив его по строкам, но ничего не вышло из него. Я запускал код с другим случайным числом семян, и проблема, похоже, не связана с тем, какое значение я присваиваю foldno. SNP в [2,4] матрицы просто не маскируется, независимо от складчатости и количества семян.


Для тех, кто заинтересован, вот пересмотренный код, который я использовал для маскировки:

CV10NA= function (input, seed, foldno) { 
set.seed(seed) 
a = unlist(input) 
b = is.na(a) 
pres = grep(FALSE, b) 
pres = sample(pres) 
sets= sample(rep(1:10, length(pres)/10), replace = FALSE) 
a[pres[which(sets==foldno)]] = NA 
a = matrix(a, ncol = ncol(input)) 
enter code here 
return(a) 
} 

ответ

1

Так sets перестановка чисел 1-10, например:

3 4 5 7 2 8 9 6 10 1 

И хотите узнать индекс значения, равный номеру сгиба:

which(sets==foldno) 

Проблема в том, что это приведет только к возвращению чисел в диапазоне от 1 до 10.Таким образом, здесь значения в beagle.subset, которые потенциально могут быть установлены в НС:

> beagle.subset[1:10] 
[1] "G" "A" "T" "A" "A" "T" "A" "G" NA "A" 

Обратите внимание, как один из них уже NA! Нет никакого способа, чтобы beagle.subset[2,4] никогда не пострадали, так как он находится в индексе 11. Вместо этого beagle.subset[3,3] будет установлен только на NA --again.

Я думаю, что вместо перетасовки 1-10 вы хотите перетасовать индексы, которые не являются NA. Затем вы можете назначить каждому из них значение 1-10, по существу, помещая их в десять ящиков.

pres = sample(pres) 
bins = seq_along(pres) %% 10 + 1 
a[pres[bins == foldno]] = NA 
+0

Спасибо за ваш отзыв @Peyton! Я пробовал ваш код, и он работает для beagle.subset. К сожалению, я попробовал код в большем наборе данных, и он не генерирует результаты, которые я хочу. Я хочу замаскировать 10% известных SNP с каждым раундом маскировки (всего 10 раундов). В конце всех 10 запусков все известные SNP должны были маскироваться ровно 1 раз. Поэтому каждый раунд должен маскировать другой набор SNP (и один набор когда-либо маскирует тот же SNP, что и другой набор). Это то, чего я надеялся достичь (но не удалось) с помощью: sets = sample (rep (1:10, length (pres)/10), replace = FALSE) – user2441436

+0

(Извините, закончились символы.) Код, который вы предоставили, маскирует только 1 SNP за раз (вместо 10%). К сожалению, я должен был предоставить больший набор данных в моем вопросе, чтобы избежать путаницы. Есть ли еще рекомендации? Еще раз спасибо за вашу помощь @Peyton! – user2441436

+0

Хорошая добыча! Я исправил это. Теперь код присваивает каждому индексу значение до 1-10. – Peyton

Смежные вопросы