☰

按每组的多数值估算数据

tag-icon • tag-icon awk perl
按每组的多数值估算数据

我在 col1 中有仪器，在 col2 中有变量，在 col3 中有调用。在此示例中，有 17 个 INS 和 3 个 M 类型。可能会有 17 x 3 行。但有些失踪了。

INS1    M1  AA
INS2    M1  AA
INS3    M1  AA
INS4    M1  GG
INS5    M1  GG
INS6    M1  GG
INS7    M1  AA
INS8    M1  GG
INS9    M1  GG
INS10   M1  AA
INS11   M1  AA
INS12   M1  GG
INS13   M1  AA
INS14   M1  AA
INS15   M1  AA
INS1    M2  GG
INS3    M2  TT
INS4    M2  GG
INS5    M2  GG
INS6    M2  TT
INS7    M2  TT
INS8    M2  TT
INS9    M2  TT
INS10   M2  GG
INS11   M2  GG
INS14   M2  GG
INS15   M2  TT
INS1    M3  AA
INS2    M3  TT
INS3    M3  AA
INS4    M3  TT
INS5    M3  TT
INS7    M3  AA
INS8    M3  TT
INS9    M3  AA
INS10   M3  TT
INS15   M3  TT

我有另一个查找文件，其中包含这些仪器所属的组。例如 INS(1,2,3,7,14,16 和 17) 属于组 1。

GR1 INS1
GR1 INS2
GR1 INS3
GR1 INS7
GR1 INS14
GR1 INS16
GR1 INS17
GR2 INS5
GR2 INS6
GR2 INS8
GR2 INS9
GR2 INS15
GR3 INS4
GR3 INS10
GR3 INS11
GR3 INS12
GR3 INS13

我正在尝试根据高于 70% 阈值的多数群体通话来估算丢失的通话。 （估算的行已加星号，但输出中不需要）。

如果在 group1 内，相同的 M 值具有 AA=80% (3/4) 的调用频率和 GG=20% (1/4) 的调用频率，那么我们可以将 group1 内所有缺失的 INS 归因于 80% 的 AA大多数人达到我们 70% 的门槛。

如果 AA 的频率为 66% (2/3)，GG 的频率为 33% (1/3)，我们不会将其估算为 AA，因为 66% 不符合 70% 阈值。

INS1    M1  AA
INS2    M1  AA
INS3    M1  AA
INS4    M1  GG
INS5    M1  GG
INS6    M1  GG
INS7    M1  AA
INS8    M1  GG
INS9    M1  GG
INS10   M1  AA
INS11   M1  AA
INS12   M1  GG
INS13   M1  AA
INS14   M1  AA
INS15   M1  AA
**INS16 M1  AA
INS17   M1  AA**
INS1    M2  GG
INS3    M2  TT
INS4    M2  GG
INS5    M2  GG
INS6    M2  TT
INS7    M2  TT
INS8    M2  TT
INS9    M2  TT
INS10   M2  GG
INS11   M2  GG
**INS12 M2  GG
INS13   M2  GG**
INS14   M2  GG
INS15   M2  TT
INS1    M3  AA
INS2    M3  TT
INS3    M3  AA
INS4    M3  TT
INS5    M3  TT
**INS6  M3  TT**
INS7    M3  AA
INS8    M3  TT
INS9    M3  AA
INS10   M3  TT
**INS11 M3  TT
INS12   M3  TT
INS13   M3  TT
INS14   M3  AA**
INS15   M3  TT
**INS16 M3  AA
INS17   M3  AA**

例如 Grp1 M1 (INS1,2,3,7,14) 有 5 次 AA 调用，频率为 100%。因此，我们可以将 M1 的 INS16 和 INS17（因为它们缺失且属于 Grp 1）推算为 AA，因为调用频率高于 70%。

对于 Grp1 M2，GG(2/4) 和 TT(2/4) 调用均以 50% 进行，因此无法在 70% 以上进行置信插补。对于 Grp1 M2 ，INS2、16 和 17 仍然缺失，因为没有明确的大多数调用（高于 70%）。

请指导如何在 awk 或 perl 中实现此目的。我尝试的解决方案是将组添加到数据中，然后找到最高调用的频率，以检查阈值，我对哈希数组有点迷失。

awk 'NR==FNR{a[$2]=$1;next} $1 in a { print $0 FS a[$1]}' groups data > tmp
awk '{count[$4 FS $1]++}END{for(j in count) print j":"count[j]}' tmp > tmp2
awk -F, '{if (a[$2]< $3)a[$2]=$3;}END{for(i in a){print i,a[i];}}' tmp2 > tmp3

答案1

这有点太复杂了，无法作为单行代码阅读，因此这里有一个带注释的gawk脚本：

#!/usr/bin/gawk -f
## Save the data in array data: data[M][INS]=dinucleotide
NR==FNR{
    data[$2][$1]=$3;
    next
}
## Save the groups in array groups: groups[GRN][INS]
{
    groups[$1][$2]++
}
## Now that everything is stored in memory, analyze
END{
    ## Get averages: for each group
    for(group in groups){
        ## For each INS in this group
        for(ins in groups[group]){
            ## For each MN in the data file
            for(m in data){
                ## If this INS had a value for this M
                if(data[m][ins]){
                    ## This counts the number of times this dinucleotide
                    ## (data[m][ins]) was found in this M among the INSs 
                    ## of this group.
                    num[group][m][data[m][ins]]++
                    ## My version of gawk doesn't seem to support
                    ## length for multidimensional arrays, so this array
                    ## only exists to count the number of Ms of this group.
                    len[group][m]++;
                }
            }
        }
    }
    ## Foreach group of the groups file
    for(group in num){
        ## For each M of this group 
        for(m in num[group]){
            ## For each INS of this group
            for(ins in groups[group]){
                ## If this INS has a value for this m in
                ## the data file, print it. 
                if(data[m][ins]){
                    printf "%-5s %s %s\n", ins,m,data[m][ins]
                }
                ## If it doesn't, check if there's an nt at
                ## >=70% for this group and print that
                else{
                    for(nt in num[group][m]){
                        if(num[group][m][nt]*100/len[group][m] >= 70){
                            printf "%-5s %s %s\n", ins,m,nt
                        }
                    }
                }
            }
        }
    }
}

将文件另存为foo.awk，使其可执行chmod +x foo.awk并在您的文件上运行：

$ ./foo.awk data groups 
INS1  M1 AA
INS2  M1 AA
INS14 M1 AA
INS3  M1 AA
INS16 M1 AA
INS17 M1 AA
INS7  M1 AA
INS1  M2 GG
INS14 M2 GG
INS3  M2 TT
INS7  M2 TT
INS1  M3 AA
INS2  M3 TT
INS14 M3 AA
INS3  M3 AA
INS16 M3 AA
INS17 M3 AA
INS7  M3 AA
INS9  M1 GG
INS15 M1 AA
INS5  M1 GG
INS6  M1 GG
INS8  M1 GG
INS9  M2 TT
INS15 M2 TT
INS5  M2 GG
INS6  M2 TT
INS8  M2 TT
INS9  M3 AA
INS15 M3 TT
INS5  M3 TT
INS6  M3 TT
INS8  M3 TT
INS10 M1 AA
INS11 M1 AA
INS12 M1 GG
INS13 M1 AA
INS4  M1 GG
INS10 M2 GG
INS11 M2 GG
INS12 M2 GG
INS13 M2 GG
INS4  M2 GG
INS10 M3 TT
INS11 M3 TT
INS12 M3 TT
INS13 M3 TT
INS4  M3 TT

请注意，此方法需要将整个数据集（两个文件）加载到内存中。不过，我确实没有找到解决方法，因为您需要先阅读整篇文章，然后才能知道是否存在 >=70% 的情况。我能想到的唯一其他方法是多次处理文件。让我知道加载到内存中是否有问题，我会看看是否可以提出不同的选项。

相关内容