#!/bin/sh
#$ -N log_flylifter
#$ -pe openmp 1
#$ -S /bin/sh
#$ -cwd
#$ -j y
#$ -q bioinf.q
. /etc/profile.d/modules.sh


## State date & time, set timer, state environment, full log, terminate on bash error.

echo "`date +%y/%m/%d_%H:%M:%S`"
START=$(date +%s.%N)
uname -sa
set -ex
find . -type f -exec ls -lt "{}" \;



## DESCRIPTION
## Shell script for conversion of D.melanogaster DRGP vcf file from reference genome build 3 to build 6.
## Also replaces DGRP variant identifiers with those from NBCI dbSNP.
## For use on a Linux Sun Grid Engine.
## Potentially applicable across species and data sizes. (~3000 genotypes per second?)
## Some DGRP variants ('mostly' indels) do not seem to have dbSNP identifiers.
## Some DGRP variants have two dbSNP rs IDs. First ID is used, second discarded.
# wpgilks@gmail.com


## Source module environment

. /etc/profile.d/modules.sh
module load jre/1.7.0_25
module load gatk/3.4-0




## Set variables.

refO=dmel-all-chromosome-r5.9.fasta
dictO=dmel-all-chromosome-r5.9.dict

refT=dm6.fa
dictT=dm6.dict

chain=mod_dm3ToDm6.over.chain
vcfO=dgrp2.vcf
vcfT=dgrp2_dm6_dbSNP.vcf




## Repair chromosome name error in chain file.

sed 's/chr//1' ${chain} > rep.${chain}




## Decompress input vcf.

gunzip ${vcfO}.gz




## Combine vcf files from NCBI, which contain the new variant IDs.

GenomeAnalysisTK \
    -R ${refO} \
    -T CombineVariants \
    -genotypeMergeOptions UNIQUIFY \
    -V vcf_chr_2L.vcf \
    -V vcf_chr_2R.vcf \
    -V vcf_chr_3L.vcf \
    -V vcf_chr_3R.vcf \
    -V vcf_chr_X.vcf \
    -V vcf_chr_Multi.vcf \
    -V vcf_chr_Un.vcf \
        -o ncbi.vcf \




## Sort NCBI and DGRP vcfs by genomic position.

for i in ncbi.vcf ${vcfO} ; do
perl vcfsorter.pl ${dictO} ${i} > sorted.${i} ; done





## Annotate sorted DGRP variant IDs with those from dbSNP.

GenomeAnalysisTK \
    -R ${refO} \
    -T VariantAnnotator \
    -V sorted.${vcfO} \
    --dbsnp sorted.ncbi.vcf \
        -o annotated.sorted.${vcfO}




## Lift over variants !!

GenomeAnalysisTK \
    -T LiftoverVariants \
    -R ${refO} \
    -V annotated.sorted.${vcfO} \
    -chain rep.${chain} \
    -dict ${dictT} \
        -o lifted.annotated.${vcfO}




## Sort lifted variants.

perl vcfsorter.pl ${dictT} lifted.annotated.${vcfO} > resorted.lifted.annotated.${vcfO}




## Filter lifted variants.

GenomeAnalysisTK \
    -R ${refT} \
    -T FilterLiftedVariants \
    -V resorted.lifted.annotated.${vcfO} \
        -o filtered.rla.${vcfO}




## Split head and body of vcf.

grep -v '^#' filtered.rla.${vcfO} > body.filtered.rla.${vcfO}
grep '^#' filtered.rla.${vcfO} > head.filtered.rla.${vcfO}




## Create file containing only variants that have dbSNP IDs.

grep 'rs' body.filtered.rla.${vcfO} > db.body.filtered.rla.${vcfO}




## Remove the DGRP IDs leaving the dbSNP IDs  (column $3).

awk 'BEGIN{FS=OFS="\t"} { gsub(/^.*;|,.*$/,"",$3) }1' \
db.body.filtered.rla.${vcfO} > clean.db.body.filtered.rla.${vcfO}




## Put header lines back on the modifed variant data. The new vcf file is now complete !

cat head.filtered.rla.${vcfO} clean.db.body.filtered.rla.${vcfO} > ${vcfT}




## Count variants in each vcf.

grep -cv ^\# *.vcf > variant_counts.txt





## Tidy-up a bit.

rm ncbi.vcf*
rm sorted.ncbi.vcf*
rm lifted.annotated.${vcfO}

rm resorted.lifted.annotated.${vcfO}*
rm filtered.rla.${vcfO}*
rm sorted.${vcfO}*

rm annotated.sorted.${vcfO}*
rm body.filtered.rla.${vcfO}
rm db.body.filtered.rla.${vcfO}

rm head.filtered.rla.${vcfO}
rm clean.db.body.filtered.rla.${vcfO}






## Compress vcfs, collect file metrics.

gzip *.vcf
ls -la *
ls -la ${refdir}
find . -type f -exec md5sum "{}" \;





## Data sources.

# dbSNP rs identifiers ftp://ftp.ncbi.nlm.nih.gov/snp/organisms/fruitfly_7227/VCF/
# dgrp2.vcf from http://dgrp2.gnets.ncsu.edu/
# dm6.fa reference genome ftp://hgdownload.cse.ucsc.edu/goldenPath/dm6/bigZips/
# GATK software https://www.broadinstitute.org/gatk/
# dm3ToDm6.over.chain from http://hgdownload-test.cse.ucsc.edu/goldenPath/dm3/liftOver/
# Publicly-available .chain not recognised, ~due to OS conflict.
# Eventually modified by reading into Python, & copy+paste. Failed methods: iconv, dos2unix, tr, perl.
# Dictionary reference built with picard/1.77
# vcf sorter perl script from https://code.google.com/p/vcfsorter/
# http://www.sussex.ac.uk/lifesci/morrowlab
## William P Gilks 26th September 2016 wpgilks@gmail.com



## End stuff.

END=$(date +%s.%N)
Duration=$(echo "$END - $START" | bc)
echo "`date +%y/%m/%d_%H:%M:%S`"
echo $Duration
exit
####
####