#
# Read inher.out produced by Allegro
#
AWK=/usr/bin/awk

if [ $# = 0 ] || [ "$1" = "-h" ] || [ "$1" = "--help" ]
then
        echo "Usage: readinh [-012zd[123]] [inheritance_file]"
        echo "       Reads a Allegro inheritance file (defaulting to inher.out)"
        echo "       and a corresponding Linkage style locus file."
        echo
        echo "       -0  Print all haplotype information."
        echo "       -1  Print haplotype information if 1 or more recombinants." 
        echo "       -2  Print haplotype information if 2 or more recombinants." 
        echo "       -dn Print list of doubly recombinant loci for deletion."
        echo "           as likely typing error, where n is the number of"                            
        echo "           consecutive markers between the recombination sites,"                         
        echo "           defaulting to 2."                            
        echo "       -zn Print list of individuals carrying more"
        echo "           than n recombinants for removal from data."
        exit
fi
INHERITANCE=inher.out
RECCNT=0
for i
do case $i in
  --double|-d) RECCNT=-d2 ;;
  -[0-9dz]*)   RECCNT=$i  ;;
          *)   if [ -f "$i" ]
               then
                 INHERITANCE=$i
               fi
   esac
done

if ! [ -f $INHERITANCE ]
then
  echo "$INHERITANCE not found!"
  exit
fi

$AWK -v reccnt=$RECCNT '
BEGIN { nloci=6
        width=0
        while(getline && NF>0) {
          if (width<length($0)) width=length($0)
          for(i=1;i<=length($0);i++) {
            col[i]=substr($0,i,1) col[i] 
          }
        }
        for(i=1;i<=width;i++) {
          gsub(/[ ]+/,"",col[i])
          if (length(col[i])>0) {
            nloci++
            locus[nloci]=col[i]
          }
          delete col[i]
        }

        reccnt=substr(reccnt,2)
        edit_script=0
        del_script=0
        long_error=3
        if (substr(reccnt,1,1)=="z") {
          del_script=1
          long_error=substr(reccnt,2,1)
        }
        if (substr(reccnt,1,1)=="d") {
          edit_script=1
          long_error=substr(reccnt,2,1)+1
        }
        name[2]="double "
        name[3]="triple "
        name[4]="quadruple "
        name[5]="quintuple "
        name[6]="6-" 
        name[7]="7-"
        name[8]="8-"
        name[9]="9-"
        name[10]="10-"

        first=7
        start=first+1
}
substr($0,1,1)!=" " && NF>0 {
        ped=$1
        cid=$2
        fid=$3
        mid=$4
        plist=""
        plast=""
        last_rec_pos=""
        prec=0
        origin=$first
        pchrom=origin
        for(i=start;i<=NF;i++) {
          if ($i!=origin) {
            prec++
            origin=$i
            if (last_rec_pos!="") gap["m",i-last_rec_pos]++
            last_rec_pos=i

            if (locus[i-1]==plast) {
              plist=plist "//" locus[i]
            }else{
              plist=plist " " locus[i-1] "//" locus[i]
            }
            plast=locus[i]
            ploc[prec]=i
          }
          pchrom=pchrom $i
        }
        getline 
        mlist=""
        mlast=""
        last_rec_pos=""
        mrec=0
        origin=$first
        mchrom=origin
        for(i=start;i<=NF;i++) {
          if ($i!=origin) {
            mrec++
            origin=$i
            if (last_rec_pos!="") gap["f",i-last_rec_pos]++
            last_rec_pos=i
            if (locus[i-1]==mlast) {
              mlist=mlist "//" locus[i]
            }else{
              mlist=mlist " " locus[i-1] "//" locus[i]
            }
            mlast=locus[i]
            mloc[mrec]=i
          }
          mchrom=mchrom $i
        }
        dummy=prec; if (dummy>9) dummy=9;  pcount[dummy]++
        dummy=mrec; if (dummy>9) dummy=9;  mcount[dummy]++
# Write output
        if (del_script && (mrec>long_error || prec>long_error)) {
          print "!"
          if (mrec>long_error) {
            print "! Maternal " name[mrec] "recombinant in " ped "-" cid 
          }
          if (prec>long_error) {
            print "! Paternal " name[prec] "recombinant in " ped "-" cid 
          }
          print "!"
          print "delete",ped,cid
        }else if (edit_script) {
          if (mrec>1) {
            for(maxrec=2;maxrec<=mrec;maxrec++) {
            if ((mloc[maxrec]-mloc[maxrec-1])<long_error) {
              printf "!\n! Maternal " name[mrec] "recombinant in " ped "-" cid "\n!\n"
              for(ll=mloc[maxrec-1];ll<mloc[maxrec];ll++) {
                print "edit",ped,cid,locus[ll],"to x x"
              }
            }}
          }
          if (prec>1) {
            for(maxrec=2;maxrec<=prec;maxrec++) {
            if ((ploc[maxrec]-ploc[maxrec-1])<long_error) {
              printf "!\n! Paternal " name[prec] "recombinant in " ped "-" cid "\n!\n"
              for(ll=ploc[maxrec-1];ll<ploc[maxrec];ll++) {
                print "edit",ped,cid,locus[ll],"to x x"
              }
            }}
          }
        }else if (prec>=reccnt || mrec>=reccnt) {
          if (ped!=lastped) {
            print "Pedigree " ped ":"
            lastped=ped
          }
          printf "%10s %10s %10s %10s %3d %s %s\n%43s %3d %s %s\n", \
                 ped, cid, fid, mid, prec, pchrom, plist,
                                " ", mrec, mchrom, mlist
        }
}
END   { print " "
        print "!"
        print "! Number of recombination events in haplotype file"
        print "!"
        print "!----------------------------------------------------------"
        print "! Number:   0    1    2    3    4    5    6    7    8    9+"
        print "!----------------------------------------------------------"

        printf "!  Male:"
        for(i=0;i<=9;i++) {
          if (pcount[i]=="") pcount[i]=0
          printf " %4d",pcount[i]
        }
        printf "\n!Female:"
        for(i=0;i<=9;i++) {
          if (mcount[i]=="") mcount[i]=0
          printf " %4d",mcount[i]
        }
        printf "\n"

        print "!"
        print "! Number of markers between sequential recombination events"
        print "!"
        print "!----------------------------------------------------------"
        print "! Number:   1    2    3    4    5    6    7    8    9    10"
        print "!----------------------------------------------------------"
        printf "!  Male:"
        for(i=1;i<=10;i++) {
          if (gap["m",i]=="") gap["m",i]=0
          printf " %4d",gap["m",i]
        }
        printf "\n!Female:"
        for(i=1;i<=10;i++) {
          if (gap["f",i]=="") gap["f",i]=0
          printf " %4d",gap["f",i]
        }
        printf "\n!\n\n"
}' $INHERITANCE