#!/bin/bash

#SBATCH --job-name=Filter_SNPs_RegioDiv_01
#SBATCH --mail-user=walter.durka@ufz.de
#SBATCH --mail-type=BEGIN,END,FAIL,TIME_LIMIT

#SBATCH --time=24:00:00
#SBATCH --mem-per-cpu=37G

## PARALLEL !! use same number of processors in Configuration file !!
#SBATCH --cpus-per-task=1

# output files
#SBATCH --output=/work/durka/%x-%j.out
#SBATCH --error=/work/durka/%x-%j.err

###########################################################################################################
# Paths and filenames for this analysis
# WORK_DIR: here are the correctly named outputfiles of process_radtags
WORK_DIR=/work/durka/dDocent/RegioDiv 

specs="ACH.MIL AGR.CAP AGR.EUP ANT.ODO ARR.ELA BIS.OFF BRO.ERE CAM.ROT CEN.JAC CYN.CRI EUP.CYP FES.RUB HYP.RAD LAT.PRA" # " TRA.AGG" #CAL.VUL KNA.ARV 

###########################################################################################################


module load Anaconda3
source $(dirname $(which conda))/../etc/profile.d/conda.sh 
conda activate /data/Popgen/programs/ddocent-2.7.8/

for species in $1 #specs
do
  cd "$WORK_DIR/$species"

# decompose complex variants into SNPs and others
vcfallelicprimitives -kg TotalRawSNPs.vcf > "$species".TotalRawSNPsPrim.vcf
##mawk '!/#/' TotalRawSNPsPrim.vcf | wc -l

# remove indels
# keep variants (SNPs), i.e. remove INDELs; remove loci not genotyped in >50% if individuals, min quality score =30 minor allele count of 3 
vcftools --gzvcf "$species".TotalRawSNPsPrim.vcf --remove-indels --max-missing 0.5 --mac 3 --minQ 30 --recode --recode-INFO-all --out "$species".SNPs.g5mac3 

# keep genotypes with >= 3 reads and only biallelic SNPs
vcftools --vcf "$species".SNPs.g5mac3.recode.vcf --minDP 3 --min-alleles 2 --max-alleles 2 --recode --recode-INFO-all --out "$species".SNPs.g5mac3minDP3bi

# output observed and expected heterozygosities and statistics
vcftools --gzvcf "$species".SNPs.g5mac3minDP3bi.recode.vcf --hardy
# search for loci that have an excess of heterozygotes P_HET_EXCESS <1e-15
cat 'out.hwe' | awk '{if($8<1e-15){print $0}}' | cut -f1-2 > het_excess15.loci 
# remove them
vcftools --gzvcf "$species".SNPs.g5mac3minDP3bi.recode.vcf --exclude-positions het_excess15.loci --recode --recode-INFO-all --out "$species".SNPs.g5mac3minDP3birl15

# output a list of all individuals with the percentage of missing genotypes (default filename "out.imiss")
vcftools --gzvcf "$species".SNPs.g5mac3minDP3birl15.recode.vcf --missing-indv 
# search for individuals that have > X missing genotypes
cat 'out.imiss' | awk '{if($5>0.75){print $0}}' | cut -f1 > lowDP75.indv 
# remove them
vcftools --gzvcf "$species".SNPs.g5mac3minDP3birl15.recode.vcf --remove lowDP75.indv --recode --recode-INFO-all --out "$species".SNPs.g5mac3minDP3birl15rm75

# filter loci by allele balance
vcffilter -s -f  "AB > 0.2 & AB < 0.8 | AB < 0.01 | AB > 0.99" "$species".SNPs.g5mac3minDP3birl15rm75.recode.vcf > "$species".SNPs.g5mac3minDP3birl15rm75.recode.F1.vcf
# filter out SNPs that have reads from both stands
vcffilter -f "SAF / SAR > 100 & SRF / SRR > 100 | SAR / SAF > 100 & SRR / SRF > 100" -s "$species".SNPs.g5mac3minDP3birl15rm75.recode.F1.vcf > "$species".SNPs.g5mac3minDP3birl15rm75.recode.F2.vcf
# filter out loci hat have strongly different mapping qualities of two alleles
vcffilter -f "MQM / MQMR > 0.9 & MQM / MQMR < 1.05" "$species".SNPs.g5mac3minDP3birl15rm75.recode.F2.vcf > "$species".SNPs.g5mac3minDP3birl15rm75.recode.F3.vcf
# filter out loci if ref and alt allels are differently paired (ref = paired; alt = non paired)
vcffilter -f "PAIRED > 0.05 & PAIREDR > 0.05 & PAIREDR / PAIRED < 1.75 & PAIREDR / PAIRED > 0.25 | PAIRED < 0.05 & PAIREDR < 0.05" -s "$species".SNPs.g5mac3minDP3birl15rm75.recode.F3.vcf > "$species".SNPs.g5mac3minDP3birl15rm75.recode.F4.vcf

# keep SNPs genotyped in >66% of individuals; mean read depth for genotype call = 5/10/20
vcftools --gzvcf "$species".SNPs.g5mac3minDP3birl15rm75.recode.vcf --max-missing 0.66 --maf 0.05 --min-meanDP  5 --max-meanDP 1000 --recode --recode-INFO-all --out "$species".SNPs.g66mac3minDP3birl15rm75.meanDP05maxDP1000maf05
vcftools --gzvcf "$species".SNPs.g5mac3minDP3birl15rm75.recode.vcf --max-missing 0.66 --maf 0.05 --min-meanDP 10 --max-meanDP 1000 --recode --recode-INFO-all --out "$species".SNPs.g66mac3minDP3birl15rm75.meanDP10maxDP1000maf05
vcftools --gzvcf "$species".SNPs.g5mac3minDP3birl15rm75.recode.vcf --max-missing 0.66 --maf 0.05 --min-meanDP 20 --max-meanDP 1000 --recode --recode-INFO-all --out "$species".SNPs.g66mac3minDP3birl15rm75.meanDP20maxDP1000maf05

vcftools --gzvcf "$species".SNPs.g5mac3minDP3birl15rm75.recode.F4.vcf --max-missing 0.66 --maf 0.05 --min-meanDP  5 --max-meanDP 1000 --recode --recode-INFO-all --out "$species".SNPs.g66mac3minDP3birl15rm75.F4.meanDP05maxDP1000maf05
vcftools --gzvcf "$species".SNPs.g5mac3minDP3birl15rm75.recode.F4.vcf --max-missing 0.66 --maf 0.05 --min-meanDP 10 --max-meanDP 1000 --recode --recode-INFO-all --out "$species".SNPs.g66mac3minDP3birl15rm75.F4.meanDP10maxDP1000maf05
vcftools --gzvcf "$species".SNPs.g5mac3minDP3birl15rm75.recode.F4.vcf --max-missing 0.66 --maf 0.05 --min-meanDP 20 --max-meanDP 1000 --recode --recode-INFO-all --out "$species".SNPs.g66mac3minDP3birl15rm75.F4.meanDP20maxDP1000maf05

# keep only one SNP per contig; however, no choice on which one, likely the first
vcftools --gzvcf "$species".SNPs.g66mac3minDP3birl15rm75.meanDP05maxDP1000maf05.recode.vcf --thin 1000 --recode --recode-INFO-all --out "$species".SNPs.g66mac3minDP3birl15rm75.meanDP05maxDP1000maf05thin
vcftools --gzvcf "$species".SNPs.g66mac3minDP3birl15rm75.meanDP10maxDP1000maf05.recode.vcf --thin 1000 --recode --recode-INFO-all --out "$species".SNPs.g66mac3minDP3birl15rm75.meanDP10maxDP1000maf05thin
vcftools --gzvcf "$species".SNPs.g66mac3minDP3birl15rm75.meanDP20maxDP1000maf05.recode.vcf --thin 1000 --recode --recode-INFO-all --out "$species".SNPs.g66mac3minDP3birl15rm75.meanDP20maxDP1000maf05thin

vcftools --gzvcf "$species".SNPs.g66mac3minDP3birl15rm75.F4.meanDP05maxDP1000maf05.recode.vcf --thin 1000 --recode --recode-INFO-all --out "$species".SNPs.g66mac3minDP3birl15rm75.F4.meanDP05maxDP1000maf05thin
vcftools --gzvcf "$species".SNPs.g66mac3minDP3birl15rm75.F4.meanDP10maxDP1000maf05.recode.vcf --thin 1000 --recode --recode-INFO-all --out "$species".SNPs.g66mac3minDP3birl15rm75.F4.meanDP10maxDP1000maf05thin
vcftools --gzvcf "$species".SNPs.g66mac3minDP3birl15rm75.F4.meanDP20maxDP1000maf05.recode.vcf --thin 1000 --recode --recode-INFO-all --out "$species".SNPs.g66mac3minDP3birl15rm75.F4.meanDP20maxDP1000maf05thin

cp "$species".SNPs.g66mac3minDP3birl15rm75.F4.meanDP10maxDP1000maf05.recode.vcf /gpfs1/data/Popgen/2021_RegioDiv/data.vcf/
cp "$species".SNPs.g66mac3minDP3birl15rm75.F4.meanDP05maxDP1000maf05thin.recode.vcf /gpfs1/data/Popgen/2021_RegioDiv/data.vcf/
cp "$species".SNPs.g66mac3minDP3birl15rm75.F4.meanDP10maxDP1000maf05thin.recode.vcf /gpfs1/data/Popgen/2021_RegioDiv/data.vcf/
cp "$species".SNPs.g66mac3minDP3birl15rm75.F4.meanDP20maxDP1000maf05thin.recode.vcf /gpfs1/data/Popgen/2021_RegioDiv/data.vcf/
cp TotalRawSNPs.vcf /gpfs1/data/Popgen/2021_RegioDiv/data.vcf/"$species".TotalRawSNPs.vcf
cp "$species".SNPs.g5mac3minDP3birl15rm75.recode.F4.vcf /gpfs1/data/Popgen/2021_RegioDiv/data.vcf/

cp lowDP75.indv /gpfs1/data/Popgen/2021_RegioDiv/Logfiles/dDocent_filter_logs/"$species".lowDP75.indv
cp out.imiss    /gpfs1/data/Popgen/2021_RegioDiv/Logfiles/dDocent_filter_logs/"$species".out.imiss
cp het_excess15.loci /gpfs1/data/Popgen/2021_RegioDiv/Logfiles/dDocent_filter_logs/"$species".het_excess15.loci
cp dDocent_main.LOG /gpfs1/data/Popgen/2021_RegioDiv/Logfiles/dDocent_main_logs/"$species".dDocent_main.LOG
cp namelist         /gpfs1/data/Popgen/2021_RegioDiv/Logfiles/dDocent_filter_logs/"$species".namelist
mkdir /gpfs1/data/Popgen/2021_RegioDiv/reference/"$species"/
cp reference.fasta /gpfs1/data/Popgen/2021_RegioDiv/reference/"$species"/

files=(TotalRawSNPs.vcf
    ${species}.TotalRawSNPsPrim.vcf
    ${species}.SNPs.g5mac3.recode.vcf
    ${species}.SNPs.g5mac3minDP3bi.recode.vcf
    ${species}.SNPs.g5mac3minDP3birl15.recode.vcf
    ${species}.SNPs.g5mac3minDP3birl15rm75.recode.vcf
    ${species}.SNPs.g5mac3minDP3birl15rm75.recode.F1.vcf
    ${species}.SNPs.g5mac3minDP3birl15rm75.recode.F2.vcf
    ${species}.SNPs.g5mac3minDP3birl15rm75.recode.F3.vcf
    ${species}.SNPs.g5mac3minDP3birl15rm75.recode.F4.vcf
    ${species}.SNPs.g66mac3minDP3birl15rm75.F4.meanDP05maxDP1000maf05.recode.vcf
    ${species}.SNPs.g66mac3minDP3birl15rm75.F4.meanDP10maxDP1000maf05.recode.vcf
    ${species}.SNPs.g66mac3minDP3birl15rm75.F4.meanDP20maxDP1000maf05.recode.vcf
    ${species}.SNPs.g66mac3minDP3birl15rm75.F4.meanDP05maxDP1000maf05thin.recode.vcf
    ${species}.SNPs.g66mac3minDP3birl15rm75.F4.meanDP10maxDP1000maf05thin.recode.vcf
    ${species}.SNPs.g66mac3minDP3birl15rm75.F4.meanDP20maxDP1000maf05thin.recode.vcf
    )

  for file in ${files[*]}
  do
    echo "Logging loci and stacks in ${file}."
    echo "${file} $(grep -Eo '^dDocent\w*' ${file} | wc -l) $(grep -Eo '^dDocent\w*' ${file} | uniq | wc -l)" >> vcf.files_N.loci.stacks.txt
  done

  cp vcf.files_N.loci.stacks.txt /gpfs1/data/Popgen/2021_RegioDiv/Logfiles/dDocent_filter_logs/${species}.vcf.files_N.loci.stacks.txt



#grep '^dDocent' TotalRawSNPs.vcf | wc -l >vcf.files_N.loci.txt
# grep '^dDocent' "$species".TotalRawSNPsPrim.vcf | wc -l >>vcf.files_N.loci.txt
 # keep variants (SNPs), i.e. remove INDELs; remove loci not genotyped in >50% if individuals, min quality score =30 minor allele count of 3 
# grep '^dDocent' "$species".SNPs.g5mac3.recode.vcf  | wc -l >>vcf.files_N.loci.txt
 # keep genotypes with >= 3 reads and only biallelic SNPs
# grep '^dDocent' "$species".SNPs.g5mac3minDP3bi.recode.vcf  | wc -l >>vcf.files_N.loci.txt
 # search for loci that have an excess of heterozygotes P_HET_EXCESS <1e-15 # remove them
# grep '^dDocent' "$species".SNPs.g5mac3minDP3birl15.recode.vcf  | wc -l >>vcf.files_N.loci.txt
 # search for individuals that have > X missing genotypes
# grep '^dDocent' "$species".SNPs.g5mac3minDP3birl15rm75.recode.vcf  | wc -l >>vcf.files_N.loci.txt
 # filter loci by allele balance
# grep '^dDocent' "$species".SNPs.g5mac3minDP3birl15rm75.recode.F1.vcf | wc -l >>vcf.files_N.loci.txt
 # filter out SNPs that have reads from both stands
# grep '^dDocent' "$species".SNPs.g5mac3minDP3birl15rm75.recode.F2.vcf | wc -l >>vcf.files_N.loci.txt
 # filter out loci hat have strongly different mapping qualities of two alleles
# grep '^dDocent' "$species".SNPs.g5mac3minDP3birl15rm75.recode.F3.vcf | wc -l >>vcf.files_N.loci.txt
 # filter out loci if ref and alt allels are differently paired (ref = paired; alt = non paired)
# grep '^dDocent' "$species".SNPs.g5mac3minDP3birl15rm75.recode.F4.vcf | wc -l >>vcf.files_N.loci.txt

 # keep SNPs genotyped in >66% of individuals; mean read depth for genotype call = 5/10/20
# grep '^dDocent' "$species".SNPs.g66mac3minDP3birl15rm75.F4.meanDP05maxDP1000maf05.recode.vcf | wc -l >>vcf.files_N.loci.txt
# grep '^dDocent' "$species".SNPs.g66mac3minDP3birl15rm75.F4.meanDP10maxDP1000maf05.recode.vcf | wc -l >>vcf.files_N.loci.txt
# grep '^dDocent' "$species".SNPs.g66mac3minDP3birl15rm75.F4.meanDP20maxDP1000maf05.recode.vcf | wc -l >>vcf.files_N.loci.txt

 # keep only one SNP per contig; however, no choice on which one, likely the first
# grep '^dDocent' "$species".SNPs.g66mac3minDP3birl15rm75.F4.meanDP05maxDP1000maf05thin.recode.vcf | wc -l >>vcf.files_N.loci.txt
# grep '^dDocent' "$species".SNPs.g66mac3minDP3birl15rm75.F4.meanDP10maxDP1000maf05thin.recode.vcf | wc -l >>vcf.files_N.loci.txt
# grep '^dDocent' "$species".SNPs.g66mac3minDP3birl15rm75.F4.meanDP20maxDP1000maf05thin.recode.vcf | wc -l >>vcf.files_N.loci.txt

# cp vcf.files_N.loci.txt /gpfs1/data/Popgen/2021_RegioDiv/Logfiles/dDocent_filter_logs/"$species".vcf.files_N.loci.txt




 #     for i in $f/*gz
 #     do


done
#####################################################################################################################
conda deactivate

