###########################################################################
## DOWNLOAD C4 REFERENCE PANEL AND REFERENCE PANEL HG19 POSITIONS        ##
###########################################################################

wget http://mccarrolllab.com/wp-content/uploads/2014/12/MHC_haplotypes_CEU_HapMap3_ref_panel.bgl
wget https://personal.broadinstitute.org/giulio/c4/MHC_haplotypes_CEU_HapMap3_ref_panel.hg19

###########################################################################
## LIFTOVER HG19 POSITIONS TO HG38                                       ##
###########################################################################

wget -q http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/liftOver
chmod a+x liftOver
wget -q http://hgdownload.cse.ucsc.edu/goldenPath/hg19/liftOver/hg19ToHg38.over.chain.gz
awk 'NR>1 {print "chr6\t"$2-1"\t"$2"\t"$1}' MHC_haplotypes_CEU_HapMap3_ref_panel.hg19 |
  ./liftOver /dev/stdin hg19ToHg38.over.chain.gz /dev/stdout /dev/stderr |
  awk 'BEGIN {print "id\tpos"} {print $4"\t"$3}' > MHC_haplotypes_CEU_HapMap3_ref_panel.hg38

###########################################################################
## DOWNLOAD HUMAN GENOME REFERENCE                                       ##
###########################################################################

wget -O- ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.gz |
  gzip -d > human_g1k_v37.fasta
samtools faidx human_g1k_v37.fasta

wget -O- ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz |
  gzip -d > GCA_000001405.15_GRCh38_no_alt_analysis_set.fna
samtools faidx GCA_000001405.15_GRCh38_no_alt_analysis_set.fna

###########################################################################
## INSTALL BCFTOOLS (IF NOT INSTALLED ALREADY)                           ##
###########################################################################

mkdir -p $HOME/bin
git clone --branch=develop git://github.com/samtools/htslib.git
git clone --branch=develop git://github.com/samtools/bcftools.git
cd htslib && autoheader && autoconf && autoconf && ./configure --disable-bz2 --disable-lzma && make && cd ..
cd bcftools && make && cd ..
/bin/cp bcftools/{bcftools,plugins/fixref.so} $HOME/bin/

###########################################################################
## GENERATE VCF FILES FOR HG19 AND HG38                                  ##
###########################################################################

declare -A fasta=( ["hg19"]="human_g1k_v37.fasta" ["hg38"]="GCA_000001405.15_GRCh38_no_alt_analysis_set.fna" )
declare -A chrom=( ["hg19"]="6" ["hg38"]="chr6" )
declare -A length=( ["hg19"]="171115067" ["hg38"]="170805979" )

for ref in hg19 hg38; do
  (echo "##fileformat=VCFv4.1"
  echo "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">"
  echo "##contig=<ID=${chrom[$ref]},length=${length[$ref]}>"
  tr '\r' '\n' < MHC_haplotypes_CEU_HapMap3_ref_panel.bgl |
    grep C4 | cut -f3- | tr '\t' '\n' | sort | uniq | awk '{print "##ALT=<ID="$0">"}'
  echo "##reference=${fasta[$ref]}"
  echo -en "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT"
  tr '\r' '\n' < MHC_haplotypes_CEU_HapMap3_ref_panel.bgl |
    paste MHC_haplotypes_CEU_HapMap3_ref_panel.$ref - |
    awk -v chr="${chrom[$ref]}" 'NR==1 {for (i=5; i<NF; i+=2) printf "\t"$i; printf "\n"}
    NR>1 {printf chr"\t"$2"\t"$4; delete x; delete y; j=0
    for (i=5; i<=NF; i++) if (!($i in x)) {x[$i]=j; y[j++]=$i}
    if ($4=="C4") {printf "\tG\t<"y[0]">"; for (i=1; i<j; i++) printf ",<"y[i]">" }
    else { printf "\t"y[0]"\t"y[1]; for (i=2; i<j; i++) printf ","y[i] }
    printf "\t.\t.\t.\tGT"
    if ($4=="C4") for (i=5; i<NF; i+=2) printf "\t"1+x[$i]"|"1+x[$(i+1)]
    else for (i=5; i<NF; i+=2) printf "\t"x[$i]"|"x[$(i+1)]
    printf "\n"}') |
    $HOME/bin/bcftools +$HOME/bin/fixref.so --no-version -Ov \
    -o MHC_haplotypes_CEU_HapMap3_ref_panel.$ref.vcf -- \
    --fasta-ref ${fasta[$ref]} --mode flip

  $HOME/bin/bcftools norm --check-ref w --fasta-ref ${fasta[$ref]} MHC_haplotypes_CEU_HapMap3_ref_panel.$ref.vcf -o /dev/null 2>&1 |
    grep REF_MISMATCH | cut -f2,3 | awk -F"\t" -v OFS="\t" 'NR==FNR {x[$2]++} NR>FNR && $2 in x {alt=$4; ref=$5; $4=ref; $5=alt
    for (i=10; i<=NF; i++) $i=1-substr($i,1,1)"|"1-substr($i,3,1)} NR>FNR' - MHC_haplotypes_CEU_HapMap3_ref_panel.$ref.vcf |
    bgzip > MHC_haplotypes_CEU_HapMap3_ref_panel.$ref.vcf.gz
done

###########################################################################
## INSTALL BEAGLE AND RUN THE IMPUTATION                                 ##
###########################################################################

declare -A chrom=( ["hg19"]="6" ["hg38"]="chr6" )
wget https://faculty.washington.edu/browning/beagle/beagle.04Jun18.a80.jar
tabix -h in.vcf.gz 6:24894177-33890574 chr6:24893949-33922797 | \
  java -Xmx8g -jar beagle.04Jun18.a80.jar gt=/dev/stdin ref=MHC_haplotypes_CEU_HapMap3_ref_panel.$ref.vcf.gz out=out \
  map=<(awk -v chr="${chrom[$ref]}" 'NR>1 {print chr"\t.\t"$2/1e7"\t"$2}' MHC_haplotypes_CEU_HapMap3_ref_panel.$ref)

###########################################################################
## INSTALL PLINK (IF NOT INSTALLED ALREADY)                              ##
###########################################################################

mkdir -p $HOME/bin
wget https://www.cog-genomics.org/static/bin/plink/plink_linux_x86_64_dev.zip
unzip -od $HOME/bin/ plink_linux_x86_64_dev.zip plink

###########################################################################
## CODE TO COMPARE WITH 1000 GENOMES PHASE 1                             ##
###########################################################################

$HOME/bin/bcftools annotate -x ID -I +'%CHROM:%POS:%REF:%ALT' MHC_haplotypes_CEU_HapMap3_ref_panel.hg19.vcf.gz |
  $HOME/bin/plink --vcf /dev/stdin --biallelic-only --keep-allele-order --const-fid --make-bed --out MHC_haplotypes_CEU_HapMap3_ref_panel.hg19
tabix -h ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase1/analysis_results/integrated_call_sets/ALL.chr6.integrated_phase1_v3.20101123.snps_indels_svs.genotypes.vcf.gz 6:24894177-33890574 |
  awk 'NF==2 {print "##contig=<ID=6,length=171115067>"} {print}' | $HOME/bin/bcftools view -v snps | $HOME/bin/bcftools annotate -x ID -I +'%CHROM:%POS:%REF:%ALT' |
  $HOME/bin/plink --vcf /dev/stdin --keep-allele-order  --const-fid --make-bed --out ALL.chr6.integrated_phase1_v3.20101123.snps_indels_svs.genotypes
$HOME/bin/plink --bfile MHC_haplotypes_CEU_HapMap3_ref_panel.hg19 --bmerge ALL.chr6.integrated_phase1_v3.20101123.snps_indels_svs.genotypes --merge-mode 6

###########################################################################
## CHECK FOR MARKERS THAT MIGHT HAVE FLIPPED                             ##
###########################################################################

awk 'NR>1 {print $1}' plink.diff | uniq -c | sort -k1,1n
awk 'NR>1 && substr($4,1,1)==substr($4,3,1) && substr($5,1,1)==substr($5,3,1) {print $1}' plink.diff | uniq -c | sort -k1,1n

###########################################################################
## CHECK A SINGLE MARKER                                                 ##
###########################################################################

tabix -f MHC_haplotypes_CEU_HapMap3_ref_panel.hg19.vcf.gz
url="ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase1/analysis_results/integrated_call_sets/ALL.chr6.integrated_phase1_v3.20101123.snps_indels_svs.genotypes.vcf.gz"
pos="31293240"
pos="33036549"
fmt="[%SAMPLE\t%POS\t%ID\t%REF\t%ALT\t%GT\n]"
bcftools query -f "$fmt" MHC_haplotypes_CEU_HapMap3_ref_panel.hg19.vcf.gz -r 6:${pos}-${pos} | sort > aswin.txt
bcftools query -f "$fmt" "$url" -r 6:${pos}-${pos} | grep -v esv | sort > 1kg.txt
join 1kg.txt aswin.txt