Check group consistencies

hello masters ,

please help here. I have 4 cols, I am looking for consistent 'geno' values within
'line', 'part' combinations. If the geno values are not consistent within a 'line', 'part' block, then we delete that block. One of the complications is that geno values are always 2 character, but AT and TA mean the same thing,similarly CG=GC etc.

I have 2 requests

1) Convert all 'geno' values in the whole dataset like AT and TA to either one or the other. I dont care which one, as long as they are consistent across the whole dataset.
2) Delete inconsistent 'line', 'part' combinations (red) from the data. If it is consistent (green) then just keep any one, in the example I have kept the first one.

line part serial geno
ax211 part1 1234 AA
gf345  part1 1345 TT
gf345  part1  3456 AA
gf345  part1 1346 TT
ax211 part2 1834 AA
gf345  part2 1395 TT
gf345  part2  3656 AA
gf345  part2 13746 TT
ax211 part3 1634 AA
gf345  part3 13345 AT
gf345  part3  34256 TA
gf345  part3 13446 AT

expected output

line part serial geno
ax211 part1 1234 AA
ax211 part2 1834 AA
ax211 part3 1634 AA
gf345  part3 13345 AT

For the first part I can try the following

sed "s/TA/AT/g;s/TC/CT/g;s/TG/GT/g;s/GC/CG/g;s/GA/AG/g;s/CA/AC/g;" file

but is there a smarter way to do this?

for the second part this is what I tried,

awk     'FNR==NR        {a[$1$2]+=$4;  b[$1$2]=$4;next}
                    $1$2 in b  {if (a[$1$2] ==1 ) print $0 }
        ' file file

Being a bit short on time i will deal only with the first part:

No, your solution is good. It is possible to do the same in awk , though, and if you want to use awk anyway for the second part of your requirement you might want to do all in one pass instead of doing something like sed '....' file | awk '....' . This, of course, is only the case if your two requirements are always two parts of one step always done together. If not, this point is moot.

I would make the regexp a bit more robust, though (here just the first rule as example):

sed "s/TA/AT/g" file

This will change every "TA" to "AT". But i fact you are interested only in changing the "TA" at the end of the line. This will also improve your code because accidental changes somewhere in the middle of a line can't happen any more. Therefore:

sed "s/TA$/AT/g" file

For the same reason you can drop the "g" option of the substitution: you want to change only one "TA" per line (the one immediately before the line end), hence:

sed "s/TA$/AT/" file

Finally, as long as it is possible you should enclose sed-expressions in single quotes instead of double quotes. In your special case it doesn't matter but some meta-characters are expanded within double-quotes whereas in single quotes nothing is expanded at all. Therefore, finally:

sed 's/TA$/AT/' file

One last tip: i usually write sed-scripts line per line, it is easier to read this way, IMHO:

sed 's/TA$/AT/
     s/TC$/CT/
     s/TG$/GT/
     s/GC$/CG/
     s/GA$/AG/
     s/CA$/AC/' file

I hope this helps.

bakunin

thanks, very helpful explanation ! .. I was wondering if there is a shortcut to reverse matching the 'geno' values without using specific sed arguments..like converting all generic xy$ and yx$ values to xy ... waiting on some part 2 help

Is your file sorted? How many rows?

---------- Post updated 01-11-15 at 03:50 PM ---------- Previous update was 01-10-15 at 10:40 PM ----------

A very crude solution, maybe the experts can help with a better one, would be problematic if you file is huge..worth a try

awk     'FNR==NR {a[$1" "$2]++;  
          b[$1" "$2" "$4]++;
  next}
                    $1" "$2 in a{
   $4 = ( a[$1" "$2] == b[$1" "$2" "$4] ) ? $4 : "multiple"; 
   delete a[$1" "$2] ;
   delete b[$1" "$2" "$4] ;
   print ;
   }' file file | awk '$4!="multiple"'

Another way:

awk '
  NR==1{
    print
    next
  }
  {
    i=$1 SUBSEP $2
    if(i in A) {
      if(B!=$4 && C!=$4) N
    }
    else {
      A=$0
      B=$4
      C=substr($4,2) substr($4,1,1)
    }
  }
  END {
    for(i in A)
      if(!(i in N))
        print A
  }
' file

Hi senhia83,

thank you, your script may be problematic since the file is being read twice into memory..

Hi Scrutinizer, will your script work on unsorted data? My file is 23Gb and when I try to use

sort -k2,2 -k1,1 

I`m getting

sort: write failed: /tmp/sortmcTb: No space left on device

any alternative suggestions on how to prepare the data?

Also, I do t have a header line , should I just ignore the part?

NR==1{
    print
    next
  }

Yes, Scrutinizer's script works with unsorted data.
Yes, this section prints the 1st line then jumps to the next line.
The script takes lots of memory that is processed in the END section, after reading the whole file.
Senhias script reads the file twice that takes time, but maybe takes less memory.

---------- Post updated at 01:14 PM ---------- Previous update was at 12:24 PM ----------

The following is a little improvement regarding memory consumption and speed:

awk '
{
  i=$1 FS $2
  if (i in A) {
    if (!((i,$4) in B)) {A=""}
  } else {
    A=$0
    B[i,$4]; B[i,substr($4,2,1) substr($4,1,1)]
  }
}
END {for (i in A) if (A!="") print A}
' file

Doesn't need a sorted file.

---------- Post updated at 01:59 PM ---------- Previous update was at 01:14 PM ----------

The following uses minimum memory but is a little slower:

awk '
{
  i=$1 FS $2
  if (i in A) {
    if (B!=$4 && B!=substr($4,2) substr($4,1,1)) A=""
  } else {
    A=$3
    B=$4
  }
}
END {for (i in A) if (A!="") print i,A,B}
' file

This should take less memory:

wak '
  {
    i=$1 SUBSEP $2
  }
  NR==FNR {
    if(i in A) {
      if(A && A!=$4 && A!=substr($4,2) substr($4,1,1)) A=x
    }
    else {
      A=$4
    }
    next
  }
  i in A && A {
    print
    A=x
  }
' file file

--
Regarding

sort: write failed: /tmp/sortmcTb: No space left on device

Have a look at sort's -T option