Shell Programming and Scripting

Hello,
I have a nexus alignment file that looks like this:


bar101_min2covg_binarynex 11001-100111
bar102_min2covg_binarynex 110010010011
bar103_min2covg_binarynex 11101010--11
etc.

There are 41 rows and 28014 characters in each, with 0, 1, and missing data (-) as the three possibilities. Probably 80% of all the sites are invariant, and I would like to remove them from the alignment. So, I'm looking for a way to scan through this alignment file and remove all sites where all rows' values match, or where only 1 row differs, ignoring missing datapoints to make this determination (i.e. if several rows have missing data at a site but all the others match, it gets chopped). A slight complication is that the data come in pairs, so I need to evaluate sites 1/2, 3/4, 5/6, 7/8, etc. etc. in pairs and eliminate them only if both sites are invariant across all rows. I'm kind of stumped at how to approach this, and fairly new to this kind of data manipulation. Does anyone have suggestions for how I might approach this?

The ideal output from the example would be:

bar101_min2covg_binarynex 001001
bar102_min2covg_binarynex 000100
bar103_min2covg_binarynex 1010--

Thanks for the help!

I think you going to need gawk for this as awk can't handle lines with more that about 3000 characters.

Try:

Code:

gawk '
{
   key[NR]=$1;
   for(i=1;i<length($2);i+=2) {
       site[i,NR]=substr($2,i,2)
       if(i>maxi)maxi=i
   }
}
END {
    c=0
    for(i=1;i<=maxi;i+=2) {
           v=""
           for(r=1;r<=NR&&!keep[i];r++) {
              if(v=="" && !(site[i,r]~"-"))v=site[i,r];
              if(length(v)&&site[i,r]!=v) keep[++c]=i;
           }
        }
        for(r=1;r<=NR;r++) {
           printf "%s ", key[r]
           for(i=1;i<c;i++) printf "%s", site[keep[i],r]
           printf "\n"
        }
}' infile

Thanks for the reply, I tried that code and it duplicated every site on each line, so the output is 56024 sites (I was off by 2 before when I wrote 28014).

Funny, it seems to be working fine for me here using the test file you posted with gawk version 4.01.

Can you try it with your 3 line test file and see how it goes?

Alright, it works ok on the example, but when it gets more complicated it starts to not handle the missing characters in the same way, and ultimately breaks down completely with a few full-length individuals, that's when the duplication of every site occurs.

Why are those first 2 pairs included? According to your rules, if only one row differs, a pair should be excluded.

Regards,
Alister

You're right, this is just a short snippet taken from the dataset so I didn't follow that rule in the example. If just one sample of the 41 differs, it's not informative, but it's not the end of the world to leave it in the dataset either.