从后验平均基因型计算等位基因频率

Question

在以下文件 (input_file.txt) 中：

col1, col2, col3, NA, NA, NA, NA, 0, 0.1, 0, 0.002, 1, 0, NA, NA, NA, 1, 1.9, 0, 0.1, 2, 1, 0, 0.8

后验平均基因型 (PMG) 为 0-2，或 NA（前 3 列不包含 PMG）。

因此，input_file.txt 有 21 个人的 PMG 信息（=字段数 (NF) -3）。

要根据input_file.txt中包含的信息计算等位基因频率，我们首先需要做一些预处理并将一些值保存在内存中。

如前所述，PMG 的范围为 0-2，但这些值中的每一个都代表三种基因型之一：

0 至 <0.5 = 纯合参考
0.5 至 <1.5 = 杂合子
1.5 到 2 = 纯合子交替

因此我们接下来需要将 PMG 转换为 input_file.txt

中的基因型

类似于：

for NF >3 && column value !=NA 
if "column value" = 0 to <0.5, replace value with 0 (= homozygous reference)
if "column value" = 0.5 to <1.5, replace value with 1 (= heterozygous)
if "column value" = 1.5 to 2, replace value with 2 (= homozygous alternate)

中间输出：

col1, col2, col3, NA, NA, NA, NA, 0, 0, 0, 0, 1, 0, NA, NA, NA, 1, 2, 0, 0, 2, 1, 0, 1

所以我们有 21 个人中的 14 个人的非 NA PMG。

我可以使用以下方法确定计数：

awk -F, ' NF>3 {print (split([=13=],a,"NA")-1) }' input_file.txt > NA_count
awk -F, ' NF>3 {print (split([=13=],a,"0")-1) }' input_file.txt > ref_count
awk -F, ' NF>3 {print (split([=13=],a,"1")-1) }' input_file.txt > het_count
awk -F, ' NF>3 {print (split([=13=],a,"2")-1) }' input_file.txt > alt_count

要从此数据中获取等位基因频率： (ref_count + het_count + 2alt_count)/2(ref_count + het_count + alt_count )

这与： 0+0+0+0+1+0+1+2+0+0+2+1+0+1/2(14) = 8/28 = 0.286 = 等位基因频率

期望的输出：

"col1", 8/28, 0.286

请注意，输入文件可以有可变数量的列，其中将有可变数量的 0-2 和 NA 值。

上面显示的中间输出本身并不是必需的（只需要所需的输出）

首选 awk 代码但不是必需的

Answer 1

awk 救援！

$ awk -F', *' -v OFS=', ' '
       {for(i=4;i<=NF;i++) if($i!="NA") a[$i=int($i+0.5)]++}
   1;
   END {print "\"""\"", (n=a[1]+2*a[2])"/"(d=2*(a[0]+a[1]+a[2])),n/d}' file

col1, col2, col3, NA, NA, NA, NA, 0, 0, 0, 0, 1, 0, NA, NA, NA, 1, 2, 0, 0, 2, 1, 0, 1
"col1", 8/28, 0.285714

要过滤中间结果，请从代码中删除 1;。

从后验平均基因型计算等位基因频率

Calculate allele frequencies from posterior mean genotypes

arrays

awk

bioinformatics