【发布时间】:2019-09-27 04:54:35
【问题描述】:
我有一个 gds 文件,描述了与参考基因组相关的许多个体的 SNP 变体。我使用了 R 中 SeqVarTools 包中的 hwe() 函数。这给了我每个变体的参考等位基因频率。我想获得次要等位基因频率,但我不知道如何解决这个问题,因为许多软件包需要将数据转换为对进一步分析无用的模糊矩阵分类。
我的主要问题:在给定参考等位基因频率的情况下,如何获得次要等位基因频率?
下面是一个小例子,可以帮助我直观地了解我的问题。
# Allele frequencies
af <- c(0.082, 0.765, 0.125, 0.986)
# Desired outcome
maf <- c(0.082, 0.235, 0.125, 0.014)
# List for outcome
maf <- c()
# Loop to take 1-af
for (i in 1:length(af)) {
if (af[i] > 0.501) {
maf[i] <- 1-af[i]
} else {maf[i] <- af[i] }
}
我正在开发的一个解决方案是一个 for 循环减去 (1 -af) if (i > 0.5) else {pass}。
我的数据集非常大,包含超过 30,000 个变量,因此 for 循环并不理想。
【问题讨论】:
-
见
maf <- af; ind <- seq(2,length(maf), by = 2);maf[ind] <- 1-maf[ind]。您似乎只是在更新偶数索引。 -
基本上,如果值超过 0.5,我需要它取 1-af。有什么方法可以实现您提供的解决方案吗?
-
我想我刚刚找到了解决这个问题的方法,我发布了一个更新作为 for 循环。如果有一种方法可以将其变成具有相同结果的更高效的代码行,那就太好了。