# for convinience
data=target
inpath=/data/module5/sim/

wc -l ${inpath}/${data}.ped
wc -l ${inpath}/${data}.map

# check where you are, and set your working directory
pwd
mkdir tmp
cd tmp/

plink \
--file ${inpath}/${data}  \
--make-bed --out ${data}

plink  --bfile ${data}  --missing   --out   ${data}
plink  --bfile ${data}  --freq      --out   ${data}
plink  --bfile ${data}  --hardy     --out   ${data}

filename="target"

imiss <- read.table(paste0(filename,".imiss"), header=T, check.names=F)
lmiss <- read.table(paste0(filename,".lmiss"), header=T, check.names=F)
freq <- read.table(paste0(filename,".frq"), header=T,    check.names=F)
hwe <- read.table(paste0(filename,".hwe"), header=T,   check.names=F)

# You can either run the plotting function in interactive Jupyter R. 
# Or un-mute the png and dev.off command lines to save the picture. 

#png(paste0(filename,"_Quality_of_genotype.png"), type="cairo")
par(mfrow=c(2,2))
hist(1-imiss$F_MISS, breaks="sturges",main="Individuals",col="tan", 
     xlab="Genotyping Rate", ylab="Number of Individuals")
hist(1-lmiss$F_MISS, breaks="sturges", main="SNPs", col="tan", 
     xlab="Genotyping Rate", ylab="Number of SNPs")
hist(hwe$P, breaks="sturges", main="HWE P-Value", col="tan", 
     xlab="HWE P-value", ylab="Number of SNPs")
hist(freq$MAF, breaks="sturges", main="MAF", col="tan", 
     xlab="MAF", ylab="Number of SNPs")
#dev.off()

# don't run this chunk with practical data. it does not have X chromosome
plink  --bfile ${data}  --check-sex --out  ${data}

## don't need to run this chunk. It's just for showing you the R package.
BiocManager::install("SNPlocs.Hsapiens.dbSNP144.GRCh37")
library(SNPlocs.Hsapiens.dbSNP144.GRCh37) 
ref <- SNPlocs.Hsapiens.dbSNP144.GRCh37

chr=22

# Pull out data for relevant chromosome and convert to VCF. 
plink --bfile ${data}_chr${chr}  --recode vcf --out  ${data}_chr${chr}

# Sort and compress the VCF file
vcf-sort ${data}_chr${chr}.vcf | bgzip -c > ${data}_chr${chr}.vcf.gz

# Fix the reference allele to match the GRCh37 reference fasta (human_g1k_v37.fasta). 
ref2fix=${refpath}/human_g1k_v37.fasta
BCFTOOLS_PLUGINS=/software/bin/
bcftools \
  +fixref \
  ${data}_chr${chr}.vcf.gz \
  -Oz \
  -o fixed_${data}_chr${chr}.vcf.gz  \
  -- -d \
  -f ${ref2fix} \
  -m flip

zcat fixed_${data}_chr${chr}.vcf.gz | bgzip -c > indexed_fixed_${data}_chr${chr}.vcf.gz

# create index file. 
tabix indexed_fixed_${data}_chr${chr}.vcf.gz

geneticmap1=${refpath}/genetic_map_chr${chr}_combined_b37_modified.txt
reference=${refpath}/ALL.chr${chr}.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz

# Use EAGLE to generate phased haplotypes
# it takes about 4 minutes.
eagle  \
    --vcfRef=$reference  \
    --vcfTarget=indexed_fixed_${data}_chr${chr}.vcf.gz  \
    --geneticMapFile=$geneticmap1  \
    --vcfOutFormat=z \
    --outPrefix=phased_chr${chr} > phasing.log

# index the vcf.gz file
tabix -p vcf  phased_chr${chr}.vcf.gz

zless phased_chr${chr}.vcf.gz | less -S

# less -S will align the columns for you, so the file is clearer to visualize. 
# press q to exit from the file

# CAUTION: genetic map file is different from the one used in phasing!
# impute5 doesn't want the chr column in genetic map, so we removed that column

# imputing chr22 takes around 8min. If you are in a hurry, try with next chunk. 

geneticmap2=${refpath}/genetic_map_chr${chr}_combined_b37.txt
reference=${refpath}/ALL.chr${chr}.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz

impute5_1.1.5_static \
        --m  $geneticmap2 \
        --h  $reference  \
        --g  phased_chr${chr}.vcf.gz  \
        --r  ${chr} \
        --ne 20000   \
        --threads 1 \
        --o  imputed_chr${chr}.vcf.gz \
        --l  imputed_chr${chr}.log

# set the boundary of the chunk to impute
intstart=40000001 
intend=50000000

impute5_1.1.5_static \
        --m  $geneticmap2 \
        --h  $reference  \
        --g  phased_chr${chr}.vcf.gz  \
        --r  ${chr}:${intstart}-${intend} \
        --ne 20000   \
        --threads 1 \
        --o  imputed_chr${chr}_chunk.vcf.gz \
        --l  imputed_chr${chr}_chunk.log

## convert the format using plink
plink --vcf  imputed_chr${chr}.vcf.gz   \
            --id-delim '_'  \
            --keep-allele-order \
            --make-bed \
            --out  imputed_chr${chr}

## get info score, need bcftools
tabix -p vcf   imputed_chr${chr}.vcf.gz
bcftools query -f '%CHROM\t%ID\t%QUAL\t%POS\t%REF\t%ALT\t%INFO/AF\t%INFO/INFO\n' \
imputed_chr${chr}.vcf.gz > imputed_chr${chr}.info

## get the list of SNPs with info > 0.3
awk ' $8 > 0.3 ' imputed_chr${chr}.info | awk '{print $2}' > snps_with_invo_over_0.3.info

## QC with info score 
plink2 --bfile imputed_chr${chr}  \
--extract snps_with_invo_over_0.3.info  \
--rm-dup force-first  \
--make-bed \
--out imputed_chr${chr}_QCed

## note that we are using PLINK2 here, since PLINK does not have function --rm-dup. 
## PLINK and PLINK2 have slight difference, and each of them have a few specific functions. We are more used to PLINK, so I kept to PLINK for all the other steps. 

## don't run this chunk. It's just an example, without usable inputs. 

# in vcf format
filelist=$(ls  imputed_chr*vcf.gz | tr '\n'   '\t')  
bcftools concat -Oz -o ${data}_imputed.vcf.gz  ${filelist}

# in PLINK format
plink \
--bfile imputed_chr1_QCed  \
--merge-list file_names_of_imputed_chr_2_to_22.txt \
--allow-no-sex  \
--make-bed  \
--out  ${data}_imputed_qced_autosomes