1

I have instruments in col1, their variables in col2 and calls in col3. In this example, there are 17 INS and 3 M types. Potentially there will be 17 x 3 rows. But some are missing.

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

I have another lookup file which has the groups, that these instruments belong to. e.g. INS(1,2,3,7,14,16 and 17) belong to group1.

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

I am trying to impute the missing calls, based on majority group call above a threshold of 70%. (imputed rows starred, but not desired in output).

If within group1 , the same M value has call frequency of AA=80% (3/4) and call frequency of GG=20% (1/4), then we can impute all the missing INS withing group1 as AA since 80% majority meets our 70% threshold.

If the frequency was 66% (2/3) for AA, and 33% (1/3) for GG, we do not impute as AA since , 66% does not meet the 70% threshold.

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**

For example Grp1 M1 (INS1,2,3,7,14) has 5 calls of AA with frequency 100%. So we can impute INS16 and INS17 (since they are missing and belonging to Grp 1) for M1 as AA because the call frequency is above 70%.

For Grp1 M2, both GG(2/4) and TT(2/4) calls are made at 50% , so a confident imputation can not be made above 70%. For Grp1 M2 , INS2,16 and 17 remain missing because there is no clear majority of a call (above 70%).

Please guide on how to achieve this in awk or perl. My attempt of a solution is to add the group to the data and then find the frequency of the highest call, to check with the threshold, I`m getting a little lost with hash arrays.

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
    How do we know which calls are missing? You've added INS16 and INS17 for M1, why not INS13 as well? For the other groups, you seem to be adding only the missing values but how can we know? Also, what's the relationship between the two files exactly? You say GR1 has INS 1-7, 14,15 and 17 but M1 has many more. Please edit and clarify, it is very hard to understand what you're asking here. – terdon Jun 26 '15 at 21:51
  • 1
    I understand this is about statistics, however is it feasible to explain this to someone who has no clue about what a "call of AA with frequency 100%" is without having them opening a statistics book? What you're describing is very technical, perhaps explaining it in simpler terms (if possible) will help you receiving a better response. – kos Jun 26 '15 at 21:58
  • Missing calls are INS x M combinations which are not there in the data. If you look at the Group lookup table there are 17 INS types and there are 3 M types in the data. So there can be 51 INS x M combinations. INS16 and INS17 was imputed for M1 because they belong to GR1 which had 4 out of 4 (INS 1,2,7,15) calls of AA at a frequency of 100%. INS13 for M1 is not missing, so it doesnt need imputation. Perhaps I couldnt explain the groups clearly, GR1 (group 1) has INS (1,2,4,7,14,16,17) GR2 has INS(5,6,8,9,15).. – Hia Sen Jun 26 '15 at 23:42
  • there is not much statistics here, its all computation. Within GR2 M1, INS (5,6,8,9) have GG, so the frequency is 4 out of a total 5 that is 80% while AA has only 1 call from INS15 with a frequency of 1/5 or 20%. – Hia Sen Jun 26 '15 at 23:48
  • 1
    Please edit the question to add more information, it's harder to read and easy to miss in the comments. So, if I understand correctly, you're also missing INS 2, 16 and 17 for M2, right? Why aren't they in your output? Or did I misunderstand? – terdon Jun 27 '15 at 0:48
2

This is a bit too complex to be legible as a one liner so here's a commented gawk script:

#!/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
                        }
                    }
                }
            }
        }
    }
}

Save the file as foo.awk, make it executable chmod +x foo.awk and run it on your files:

$ ./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

Note that this approach requires loading the entire dataset (both files) into memory. I don't really see a way around that though since you need to read the entire thing before knowing whether there are cases of >=70%. The only other approach I can think of would involve processing the file multiple times. Let me know if loading into memory is a problem and I'll see if I can come up with a different option.

  • 1
    To sort the output: ./foo.awk data groups | sort -b -k 2,2 -k 1,3V – A.B. Jun 27 '15 at 15:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.