将基因型转换为 0/1

Question

我有一个文件，看起来像：

1    rs6687776    1020428    T    C    T    C    T    C    C    C    T    C    C    C    T    C

第 4 和第 5 列是该位点的两个不同的可能等位基因。我需要更改第 6 列以显示0是否有T等位基因以及1是否有C等位基因。我的文件是 20805 x 459。所以应该看起来像：

1   rs6687776   1020428 T   C   0   1   0   1   1   1   0   1   1   1   0   1

我试过了：

cat file | while read line
do if [ [,6-] = [,4] ]
then
    echo "0"
    echo "1"
fi
done

不过，我刚刚结束了与交替的一个文件0的和1的是41610名行长。也许AWK更有用？

Answer 1

这是另一种awk方法：

$ awk '{a[$4]=0;a[$5]=1; for(i=6;i<=NF;i++){$i=a[$i]}}1;' file
1 rs6687776 1020428 T C 0 1 0 1 1 1 0 1 1 1 0 1

解释

• a[$4]=0;a[$5]=1;: 创建a具有两个键的数组，$4和$5。的值$4设置为0，的值设置$5为 1。

• for(i=6;i<=NF;i++){$i=a[$i]} ：对于从 6 到最后一个的每个字段编号，将该字段设置为数组中存储的用于找到的核苷酸的任何内容。

• 1; : awk “打印这一行”的简写。

你也可以用 Perl 来做：

$ perl -lane 's/$F[3]/0/ for @F[5..$#F]; s/$F[4]/1/ for @F[5..$#F]; print "@F"' file
1 rs6687776 1020428 T C 0 1 0 1 1 1 0 1 1 1 0 1

这是同样的想法。该-a品牌perl像awk，劈裂上空白每行到阵列中@F。然后，我们将第 4 个字段（$F[3]，数组从 0 开始）中找到的所有核苷酸情况替换为0，将第 5 个 ( $F[4]) 中的所有情况替换为1。这for @F[5..$#F]意味着替换仅适用于字段 6 到最后。最后，打印修改后的数组。