#!/bin/csh -f

# 26.May.2011 0:29:19[yunmli]codon~/projects/ESP/ref/2011-05
# get_union_ped.csh

foreach s (HeartGO WHISP)
  if (-e $s.fam && -e $s.bim && -e $s.bed) then
    echo "$s proceeding ..."
  else
    echo "ERROR: missing $s PLINK input file(s) $s.fam $s.bim $s.bed"
    exit
  endif
  
  if ($s == HeartGO) then
    set sd = ~/projects/HeartGO
  else
    set sd = ~/imputation_GWA/WHI
  endif
  
  mkdir -p esp/$s/ped
  mkdir -p esp/$s/map
  mkdir -p esp/$s/dat
  mkdir -p esp/$s/snps
  
  mkdir -p GWAS/$s/ped
  mkdir -p GWAS/$s/map
  mkdir -p GWAS/$s/dat
  mkdir -p GWAS/$s/snps
  
  mkdir -p concord/$s/log
  mkdir -p concord/$s/bySNP
  mkdir -p concord/$s/minusstrand
  mkdir -p concord/$s/highly_discordant
  
  mkdir -p union/$s/ped
  mkdir -p union/$s/dat
  mkdir -p union/$s/snps
  mkdir -p union/$s/map
  foreach chr (`seq 1 22`)
    ## ESP data ##
    set epf = esp/$s/ped/chr$chr.ped
    if (!(-e $epf.gz)) then
      echo "$s chr$chr generating ESP files"
      # extract information for chr$chr
      plink --bfile $s --chr $chr --make-bed --out chr$chr.plink --noweb >& chr$chr.plink.log
      rm -f chr$chr.plink.log
      
      # rename ESP SNPs that are also in GWAS to rs ID
      mv chr$chr.plink.bim chr$chr.plink.orig.bim
      awk '{print $4}' chr$chr.plink.orig.bim > chr$chr.plink.pos
      if ($s == HeartGO) then
        set mapf = ~/projects/HeartGO/b37/map/chr$chr.map
      else
        set mapf = ~/imputation_GWA/WHI/b37/map/chr$chr.map
      endif
      if (!(-e $mapf)) then
        echo "ERROR: $mapf does not exist"
        exit
      endif
      extractor_field 3 chr$chr.plink.pos < $mapf > chr$chr.overlap.map
      awk '{print $2}' chr$chr.overlap.map > chr$chr.overlap.snps
      awk '{print $3}' chr$chr.overlap.map > chr$chr.overlap.pos
      cat chr$chr.overlap.pos chr$chr.plink.pos | sort | uniq -c | awk '$1 == 1' | awk '{print $2}' > chr$chr.plink_only.pos
      extractor_field 4 chr$chr.plink_only.pos < chr$chr.plink.orig.bim | awk '{print $1,$2,$4}' > chr$chr.plink_only.map
      cat chr$chr.overlap.map chr$chr.plink_only.map | sort -g -k3 | awk '{print $2}' > chr$chr.plink.renamed.pos_sorted.snps
      paste chr$chr.plink.renamed.pos_sorted.snps chr$chr.plink.orig.bim | awk '{print $2,$1,$4,$5,$6,$7}' | sed -e 's/ /\t/g' > chr$chr.plink.bim
      rm -f chr$chr.plink.pos chr$chr.overlap.map chr$chr.overlap.pos chr$chr.plink_only.pos chr$chr.plink_only.map chr$chr.plink.renamed.pos_sorted.snps
      
      plink --bfile chr$chr.plink --recode --out chr$chr --noweb >& chr$chr.plink.log
      rm -f chr$chr.plink.* 
      rm -f chr$chr.nosex
      ped_convert.pl -i chr$chr.recode.ped -o esp/$s/ped/chr$chr.ped -h 5 -col 7 > chr$chr.$s.pedconvert.log
      rm -f chr$chr.$s.pedconvert.log
      rm -f chr$chr.log chr$chr.nosex chr$chr.recode.ped
      gzip esp/$s/ped/chr$chr.ped
      
      awk '{print $1,$2,$4}' chr$chr.recode.map > esp/$s/map/chr$chr.map
      rm -f chr$chr.recode.map
      
      awk '{print $2}' esp/$s/map/chr$chr.map > esp/$s/snps/chr$chr.snps
      awk '{print $2}' esp/$s/map/chr$chr.map | sed -e 's/^/M /g' > esp/$s/dat/chr$chr.dat
    endif
    
    ## GWAS data ##
    set gpf = GWAS/$s/ped/chr$chr.ped
    if (!(-e $gpf.gz)) then
      echo "$s chr$chr generating GWAS files"
      if ($s == HeartGO) then
        set pedf = ~/projects/HeartGO/b37/ped/chr$chr.ped.gz
      else
        set pedf = ~/imputation_GWA/WHI/b37/ped/chr$chr.ped.gz
      endif
      if (!(-e $pedf)) then
        echo "ERROR: $pedf does not exist"
        exit
      endif
      zcat $pedf | extractor_field 2 $s.with_ESP_and_GWAS.GWAS.ID > GWAS/$s/ped/chr$chr.ped
      gzip GWAS/$s/ped/chr$chr.ped
      cp -p $sd/b37/dat/chr$chr.dat GWAS/$s/dat/
      cp -p $sd/b37/map/chr$chr.map GWAS/$s/map/
      cp -p $sd/b37/snps/chr$chr.snps GWAS/$s/snps/
    endif

    ## check concordance ##
    set bySNPf = concord/$s/bySNP/chr$chr.bySNP
    if (!(-e $bySNPf)) then
      echo "$s chr$chr checking concordance"
      set logf = concord/$s/log/chr$chr.log
      CalcMatch --trueped $gpf.gz --truedat GWAS/$s/dat/chr$chr.dat --impped $epf.gz --impdat esp/$s/dat/chr$chr.dat -o chr$chr > $logf
      mv chr$chr.bySNP $bySNPf
      mv chr$chr.minusstrand concord/$s/minusstrand/
    endif
    
    set upf = union/$s/ped/chr$chr.ped
    if (!(-e $upf.gz)) then
      echo "$s chr$chr generating union files"
      ## (1) mask overlapping SNPs in ESP ##
        mask_dat.pl -d esp/$s/dat/chr$chr.dat -s chr$chr.overlap.snps -o chr$chr.esp.mask_in_S2.dat > chr$chr.mask.log
        rm -f chr$chr.mask.log
        ln -s esp/$s/ped/chr$chr.ped.gz chr$chr.esp.mask_in_S2.ped
      ## (2) mask highly discordant (gErr > 5%) SNPs in GWAS ##
        awk '$2 > 0.05' $bySNPf | grep -v SNP | awk '{print $1}' > concord/$s/highly_discordant/chr$chr.highly_discordant.snps
        mask_dat.pl -d GWAS/$s/dat/chr$chr.dat -s concord/$s/highly_discordant/chr$chr.highly_discordant.snps -o chr$chr.gwas.mask_in_S2.dat > chr$chr.mask.log
        rm -f chr$chr.mask.log
        ln -s GWAS/$s/ped/chr$chr.ped.gz chr$chr.gwas.mask_in_S2.ped
      ## (3) prepare order pedigree files (dat in order, ped empty) ##
        cat chr$chr.esp.mask_in_S2.dat chr$chr.gwas.mask_in_S2.dat | grep "M " | awk '{print $2}' > chr$chr.union_to_keep.snps
        cat GWAS/$s/map/chr$chr.map esp/$s/map/chr$chr.map | sort -g -k 3 | uniq > chr$chr.union.all.map
        extractor_field 2 chr$chr.union_to_keep.snps < chr$chr.union.all.map > union/$s/map/chr$chr.map
        awk '{print $2}' union/$s/map/chr$chr.map > union/$s/snps/chr$chr.snps
        rm -f chr$chr.union_to_keep.snps chr$chr.union.all.map
        awk '{print $2}' union/$s/map/chr$chr.map | sed -e 's/^/M /g' > chr$chr.union.raw.pos_ordered.dat
        touch chr$chr.union.raw.pos_ordered.ped
      ## (4) merge ##
        pedmerge_yun chr$chr.union.raw.pos_ordered chr$chr.esp.mask_in_S2 chr$chr.gwas.mask_in_S2 union/$s/ped/chr$chr
        mv union/$s/ped/chr$chr.dat union/$s/dat/
        gzip union/$s/ped/chr$chr.ped
      ## (5) clean ##
      rm -f chr$chr.esp.mask_in_S2.* 
      rm -f chr$chr.gwas.mask_in_S2.*
      rm -f chr$chr.union.raw.pos_ordered.*
      rm -f chr$chr.overlap.snps
    endif

  end
end

