1

我想使用 SNP 基因型数据确定程序的插补准确性，因此我需要屏蔽一部分 SNP 调用以模拟缺失数据。

我一直在这个标记数据子集上测试我的代码（见下文）。列名是个人的名字。行名称是 SNP ID。数据集包含缺失数据（标记为 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% 的已知SNP，总共 10 次不同的时间（即 10 轮掩蔽）。

  2. 每一轮都需要掩盖一组不同的 SNP。

  3. 每个SNP 在这 10 轮中只应被掩蔽一次（例如，SNP S10_28619 将在 10 轮掩蔽中仅显示为“NA”一次）。

这是我一直在使用的代码：

##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 的值。矩阵的 [2,4] 中的 SNP 不会屏蔽，无论 foldno 和种子数如何。

对于任何有兴趣的人，这是我用于屏蔽的修改后的代码：

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)
}
4

1 回答 1

1

sets数字 1-10 的排列也是如此，例如：

3  4  5  7  2  8  9  6 10  1

并且您想找出等于折叠数的值的索引：

which(sets==foldno)

问题是，这只会返回 1-10 范围内的数字。因此，以下是beagle.subset可能设置为 NA 的值：

> 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 的值，基本上将它们放入 10 个箱中。

pres = sample(pres)
bins = seq_along(pres) %% 10 + 1 
a[pres[bins == foldno]] = NA
于 2013-06-01T04:13:56.680 回答