#########################################
###### From aligned reads to SFSs #######  
#########################################

######	Set up useful environment variables

REF="/home/DATA10/PAOLO/Turbot_2b/Assembly/GCA_003186165.1_ASM318616v1/GCA_003186165.1_ASM318616v1_genomic.fna"
MOMENTS="/home/DATA10/PAOLO/Turbot_2b/All_Trimmed/MOMENTS"
ANGSD="/home/DATA10/PAOLO/Turbot_2b/All_Trimmed/ANGSD"
SCRIPTS="/home/egru/PAOLO/2bRAD_GATK/"
###### With the varibale AFS i link to where i have some stored some useful scripts from Mikhail Matz (https://github.com/z0on/AFS-analysis-with-moments/tree/master/multimodel_inference)
AFS="/home/egru/PAOLO/AFS-analysis-with-moments/multimodel_inference/"
######	Filter to remove bad reads, triallelic snps, and based on mapping and phred quality as well as > 60% Individuals
FILTERS="-uniqueOnly 1 -remove_bads 1  -skipTriallelic 1 -minMapQ 20 -minQ 20 -doHWE 1 -minInd 100 "
###### ANGSD command:
TODO="-doMajorMinor 1 -doMaf 1 -dosnpstat 1 -doPost 2 -doGeno 11"
HOME="/home/egru/PAOLO/"
TODO2="-doSaf 1 -anc $REF -ref $REF"



######	ANGSD command to create MAF files
cat  $MOMENTS/NS.bamlist  $MOMENTS/BS.bamlist > $MOMENTS/NS_BS.bamlist
angsd -b $MOMENTS/NS_BS.bamlist -GL 1 -P 10 $FILTERS $TODO -out $MOMENTS/NS_BS/NS_BS_sfilt

######	Remove sites with het > 0.5 (likely paralogs) and make a site index
zcat $MOMENTS/NS_BS/NS_BS_sfilt.snpStat.gz | awk '($3+$4+$5+$6)>0' | awk '($12+$13+$14+$15)/($3+$4+$5+$6)<0.7' | cut -f 1,2 >$MOMENTS/NS_BS/NS_BS_sites2do 
angsd sites index $MOMENTS/NS_BS/NS_BS_sites2do

######	In the following lines, set minInd to 90% of each pop's sample size (that's 18 individuals in the North Sea and 81 in the Baltic Sea)

angsd -sites $MOMENTS/NS_BS/NS_BS_sites2do -b $MOMENTS/NS.bamlist -GL 1 -P 10 $TODO2 -minInd 18 -out $MOMENTS/NS_BS/SITES_STRICT/NS
angsd  -sites $MOMENTS/NS_BS/NS_BS_sites2do  -b $MOMENTS/BS.bamlist -GL 1 -P 10 $TODO2 -minInd 81 -out $MOMENTS/NS_BS/SITES_STRICT/BS

######	generating per-population SFS for the North Sea and Baltic Sea individuals
realSFS $MOMENTS/NS_BS/SITES_STRICT/NS.saf.idx >$MOMENTS/NS_BS/SITES_STRICT/NS.sfs
realSFS $MOMENTS/NS_BS/SITES_STRICT/BS.saf.idx >$MOMENTS/NS_BS/SITES_STRICT/BS.sfs

######	generating dadi-like posterior counts based on sfs priors

realSFS dadi $MOMENTS/NS_BS/SITES_STRICT/NS.saf.idx $MOMENTS/NS_BS/SITES_STRICT/BS.saf.idx -sfs $MOMENTS/NS_BS/SITES_STRICT/NS.sfs -sfs $MOMENTS/NS_BS/SITES_STRICT/BS.sfs -ref $REF -anc $REF >$MOMENTS/NS_BS/SITES_STRICT/NS_BS_dadi.out

realSFS dadi $MOMENTS/NS_BS/SITES_STRICT/NS.saf.idx -sfs $MOMENTS/NS_BS/SITES_STRICT/NS.sfs -ref $REF -anc $REF >$MOMENTS/NS_BS/SITES_STRICT/NS_dadi.out
realSFS dadi $MOMENTS/NS_BS/SITES_STRICT/BS.saf.idx -sfs $MOMENTS/NS_BS/SITES_STRICT/BS.sfs -ref $REF -anc $REF >$MOMENTS/NS_BS/SITES_STRICT/BS_dadi.out

######	converting to dadi-snp format understood by dadi an Moments:
######	(numbers after the input file name are numbers of individuals sampled per population)

cd $MOMENTS/NS_BS/SITES_STRICT
perl $SCRIPTS/realsfs2dadi.pl NS_BS_dadi.out 20 91 > NS_BS_dadi.data
perl $SCRIPTS/realsfs2dadi.pl NS_dadi.out 20 > NS_dadi.data
perl $SCRIPTS/realsfs2dadi.pl BS_dadi.out 91 > BS_dadi.data

######	Now we just change the scaffolds IDs to Linkage Group numbers 

cp $MOMENTS/NS_BS/SITES_RELAXED/NS_BS_dadi.out $MOMENTS/SITES_RELAXED/NS_BS_CHR_dadi.out
while read a b; do sed -i "s,$a,CHR$b,g" $MOMENTS/SITES_RELAXED/NS_BS_CHR_dadi.out  ; done < /home/DATA10/PAOLO/Turbot_2b/All_Trimmed/Variants/Fst/CHR_CONVERSION

cp $MOMENTS/NS_BS/SITES_STRICT/NS_BS_dadi.out $MOMENTS/SITES_STRICT/NS_BS_CHR_dadi.out
while read a b; do sed -i "s,$a,CHR$b,g" $MOMENTS/SITES_STRICT/NS_BS_CHR_dadi.out  ; done < /home/DATA10/PAOLO/Turbot_2b/All_Trimmed/Variants/Fst/CHR_CONVERSION
