#!/bin/bash
## Excute from the Mac or Linux command terminal with :
## caffeinate -i sh code_drive_basic_gwas.sh > log_drive_basic_gwas.txt

## DESCRIPTION
## Bash shell script for genome-wide association study of Drosophila melanogaster reproductive fitness data for females and males, using Plink v1.9.
## William P. Gilks wpgilks@gmail.com University of Sussex 2017 http://www.sussex.ac.uk/lifesci/morrowlab/

## Assumes dm6 'bim' file chromosome names have been modified for plink compatibility.
## Downloads essential phenotype-genotype files, and annotation information files.
## Performs GWAS using Plink linear regression.
## Formats and annotates results with SNP function and gene information.
## Generate basic diagnostic graphs.

## Log all, quit on error, state environment, list files, state time.
set -ex ; uname -a ; ls -la ; echo "Start time is `date +%H:%M:%S``date +%d/%m/%y`"


## Input data download:

## Genotype data
curl -O "https://zenodo.org/record/804008/files/lhm_hc_2017.fam"
curl -O "https://zenodo.org/record/804008/files/lhm_hc_2017.bim"
curl -O "https://zenodo.org/record/804008/files/lhm_hc_2017.bed"

## Phenotypes (Female and Male fitness).
curl -O "https://zenodo.org/record/807698/files/lhm_Female_Male_plink.pheno"

## Allele frequencies.
curl -O "https://zenodo.org/record/806849/files/lhm_hc_2017.frq"

## Gene positions.
curl -O "https://zenodo.org/record/815472/files/gene_positions_plink.txt"

## NCBI-dbSNP identifiers.
curl -O "https://zenodo.org/record/815472/files/lhm_ncbi_snp_key.txt"

## SNP functions.
curl -O "https://zenodo.org/record/815472/files/lhm_snp_functions_plink.txt"


## Set variables.
geno=lhm_hc_2017
pheno=lhm_Female_Male_plink.pheno
var_funcs=lhm_snp_functions_plink.txt
genes=gene_positions_plink.txt
ncbi_ids=lhm_ncbi_snp_key.txt



## GWAS: Linear regression, qc-filtered variants from haplotypecaller, thresholds as below.
## (chr 6 and 7 refer to dm6 chromosomes 4 and mitochondria)
plink1.9 --allow-no-sex \
    --bfile ${geno} \
    --linear no-x-sex -ci 0.95 \
    --geno 0 \
    --maf 0.05 \
    --pheno ${pheno} \
    --all-pheno \
    --not-chr 6,7 \
        --out ${geno}


## Remove excess whitespace from GWAS results, and sort by SNP ID.
for i in *.assoc.linear ; do tr -s " " < ${i} | sed 's/^ *//g' | sort -k2 > ${i}.s ; done




## Edit the names of uncharacterised genes so they are deprioritized when assigning one gene per SNP.
## Remove three problematic snRNA genes with grossly extended boundaries.
sed -e 's/CG/yyCG/g' -e 's/CR/zzCR/g' ${genes} | sed -n '/snRNA:U1:21D/!p' | sed -n '/snRNA:U1:95Ca/!p' | sed -n '/snRNA:U1:95Cb/!p' > ${genes}.mod




## ANNOTATION OF SNPS.
for i in *.assoc.linear.s ; do

## Assign genes to SNPs.
plink1.9 --annotate ${i} ranges=${genes}.mod minimal --border 0 --out ${i}.g

## Edit column title. Keep name of first gene overlapping SNP. Edit CG and CR gene names back to original.
sed -e 's/ANNOT/chr_gene/g' ${i}.g.annot | cut -d\| -f1 | sed -e 's/yyCG/CG/g' -e 's/zzCR/CR/g' > ${i}.g.annot.tmp
rm ${i}.g.annot ; rm ${i}.g.log

## Assign functions to SNPs.
plink1.9 --annotate ${i}.g.annot.tmp attrib=${var_funcs} --out ${i}.gf
rm ${i}.g.annot.tmp

## Edit column title.
sed -e 's/ANNOT/snp_funct/g' ${i}.gf.annot > ${i}.gf.tmp
rm ${i}.gf.annot ; rm ${i}.gf.log ; done

rm *.linear.s




## Assign NCBI SNP IDs (Join by 1st column).
for i in *.linear.s.gf.tmp ; do join -1 2 -2 1 ${i} ${ncbi_ids} > ${i}.ncbi ; done
rm *.linear.s.gf.tmp





## For allele frequency file, remove excess whitespace and unnecessary columns, then sort by SNP ID.
tr -s " " < ${geno}.frq | awk '{$1=$1}1' | cut -d" " -f2- | sort > ${geno}.frq.mod

## Join the GWAS and allele frequency results using bash-join command. Assumes both are sorted by ID.
for i in *.ncbi ; do join ${i} ${geno}.frq.mod > ${i}.af ; done
# rm *.assoc.linear.s.gf.tmp.ncbi




## Edit chromosome names to dm6 genome assembly.
for i in *.assoc.linear.s.gf.tmp.ncbi.af ; do
awk '{ sub(/1$/,"chr2L",$2); sub(/2$/,"chr2R",$2); sub(/3$/,"chr3L",$2); sub(/4$/,"chr3R",$2); sub(/5$/,"chrX",$2) }1' < ${i} | \
cut -d" " -f1-3,6-19 > ${i%.assoc.linear.s.gf.tmp.ncbi.af}.gwas_res.txt ; done


rm *.assoc.linear.s.gf.tmp.ncbi.af
rm *.nosex
rm *.assoc.linear.s.gf.tmp.ncbi



### Do initial post-analysis of GWAS results.
R --vanilla < code_plot_basic_gwas.R | tee log_sh_plot_basic_gwas.txt



## Remove leftover files
rm ${geno}.fam ; rm ${geno}.bim
rm ${geno}.bed ; rm ${pheno}
rm ${var_funcs} ; rm ${genes}
rm ${ncbi_ids}



## Notes
## Morrow lab, Sussex Uni http://www.sussex.ac.uk/lifesci/morrowlab/
## Plink 1.9 https://www.cog-genomics.org/plink2/
## Because of RAM limit, using Plink to assign dbSNP IDs couldn't be done, so I just used 'join' command.
## NCBI dbSNP project link https://www.ncbi.nlm.nih.gov/projects/SNP/snp_viewBatch.cgi?sbid=1062461
## Genotype generation and release described here https://f1000research.com/articles/5-2644/v3

## Get file md5 checksums.
find . -type f -name "*" -exec md5 {} \; > log_md5_checksums_lhm_gwas.txt

ls -la ; echo "End time is `date +%H:%M:%S``date +%d/%m/%y`"
echo "William P. Gilks wpgilks@gmail.com University of Sussex, September 2017 http://www.sussex.ac.uk/lifesci/morrowlab/"
exit
###
##
#