


##################################
### 1 | Remove PCR duplication ###
##################################

DEDUP=/scratch/project_2000465/Whitefish/Example_files/trim2bRAD_2barcodes_dedup.pl
# "trim2bRAD_2barcodes_dedup.pl" is referred to https://github.com/z0on/2bRAD_denovo
cd /scratch/project_2000465/Whitefish/FASTQ
for i in *fq; do perl $DEDUP input=$i  sampleID=1 bc2=TGGT adaptor=AGA &>> Dedup.ou; done


################################
## 2 | Sequence reads mapping ##
################################

#!/bin/bash -l
#SBATCH --account=XXXXX
#SBATCH -o maps.ou
#SBATCH -e maps.err
#SBATCH -p small
#SBATCH -n 1
#SBATCH --cpus-per-task=4
#SBATCH -t 3-00:00:00
#SBATCH --mem-per-cpu=4000

BAM="/scratch/project_2000465/Whitefish/BAM/"
SAM="/scratch/project_2000465/Whitefish/SAM/"
tr0="/scratch/project_2000465/Whitefish/tr0"
GENOME=/scratch/project_2000465/Whitefish/GENOME/GCA_902810595.1_AWG_v2/GCA_902810595.1_AWG_v2_genomic.fna
SAMPLES=/scratch/project_2000465/Whitefish/Sample_ID_names.tab
module load biokit

while read a b
do bowtie2  --no-unal --score-min L,16,1 --local -L 16 -p 4 -x $GENOME -U $tr0/${a}.tr0  | samtools view -b  | samtools sort -o  $BAM/${b}.sorted.bam && picard AddOrReplaceReadGroups  INPUT=$BAM/${b}.sorted.bam OUTPUT= $BAM/${b}.bam RGID=group1 RGLB=lib1 RGPL=illumina RGPU=unit1 RGSM=$b && samtools index $BAM/${b}.bam
 done < $SAMPLES

rm -rf  $BAM/*.sorted.bam

#################################################
### 3 | Population genetics of all 61 samples ###
#################################################

GENOME=/scratch/project_2000465/Whitefish/GENOME/GCA_902810595.1_AWG_v2/GCA_902810595.1_AWG_v2_genomic.fna
BAM=/scratch/project_2000466/Whitefish/BAM
GENO=/scratch/project_2000466/Whitefish/GENO
SITES=/scratch/project_2000466/Whitefish/SITES
SAF=/scratch/project_2000466/Whitefish/Diversity/SAF
SFS=/scratch/project_2000466/Whitefish/Diversity/SFS
FST=/scratch/project_2000466/Whitefish/Diversity/FST
FILTERS="-uniqueOnly 1 -remove_bads 1 -minMapQ 20 -minQ 20 -C 50"
TODO="-dosaf 1 -GL 2"

module load bioconda
source activate angsd

#################################
## Calling genotype likelihood ##
#################################
# all 61 samples (pelagic + profundal + Profundal)
angsd -bam $BAM/bamlist_61IDs.txt \
-anc $GENOME -ref $GENOME \
-doMaf 2 -doGlf 2 -doMajorMinor 1 -doGeno 2 -doPost 1 -GL 1 -skipTriallelic 1 -doCounts -nThreads 16 \
$FILTERS \
-minInd 45 -minIndDepth 2 \
-SNP_pval 1e-6 \
-out $GENO/Geno_littoral_pelagic_profundal_filtered &>> $GENO/Geno_littoral_pelagic_profundal_filtered.log

#############
##-- PCA --##
#############
## PCA of 61 samples
cd /home/DATA10/BOHAO/Whitefish/PCANGSD/PCA
python $PCANGSD -beagle $GENO_WS/Geno_littoral_pelagic_profundal_filtered.beagle.gz -o Geno_littoral_pelagic_profundal_filtered_maf002 -threads 20 -sites_save -minMaf 0.02

##################
##-- NgsAdmix --##
##################
## NGSadmix of 61 samples
NGSadmix -likes $GENO_WS/Geno_littoral_pelagic_profundal_filtered.beagle.gz -K 2 -o Geno_littoral_pelagic_profundal_K2_NGSadmix -P 10 &>> Geno_littoral_pelagic_profundal_K2_NGSadmix.log
NGSadmix -likes $GENO_WS/Geno_littoral_pelagic_profundal_filtered.beagle.gz -K 3 -o Geno_littoral_pelagic_profundal_K3_NGSadmix -P 10 &>> Geno_littoral_pelagic_profundal_K3_NGSadmix.log
NGSadmix -likes $GENO_WS/Geno_littoral_pelagic_profundal_filtered.beagle.gz -K 4 -o Geno_littoral_pelagic_profundal_K4_NGSadmix -P 10 &>> Geno_littoral_pelagic_profundal_K4_NGSadmix.log


###############################################################
### 4 | Demographic history of littoral and pelgagic morphs ###
###############################################################

##  Calculate SFS ##
angsd -bam $BAM/bamlist_littoral-pelagic_pure_IDs.txt \
-anc $GENOME -ref $GENOME \
$FILTERS $TODO -minInd 11 \
-GL 1 \
-nThreads 20 \
-out $GENO/Geno_littoral_pelagic_demo &>> $GENO/Geno_littoral_pelagic_demo.log

## Remove sites with het > 0. (likely paralogs)
zcat $GENO/Geno_littoral_pelagic_demo.snpStat.gz | awk '($3+$4+$5+$6)>0' | awk '($12+$13+$14+$15)/($3+$4+$5+$6)<0.7' | cut -f 1,2 > $SITES/Geno_littoral_pelagic_demo_todo 
angsd sites index $SITES/Geno_littoral_pelagic_demo_todo

## SAF per-morph (set minInd to 90% of each pop's sample size)
angsd -bam $BAM/bamlist_littoral_pure_IDs.txt \
-anc $GENOME -ref $GENOME \
-sites $SITES/Geno_littoral_pelagic_demo_todo \
-doSaf 1 -GL 1 \
-minInd 5 \
-nThreads 20 \
-out $SAF/littoral-pelagic/littoral_pure_IDs_demo  &>> $SAF/littoral-pelagic/littoral_pure_IDs_demo.log

angsd -bam $BAM/bamlist_pelagic_pure_IDs.txt \
-anc $GENOME -ref $GENOME \
-sites $SITES/Geno_littoral_pelagic_demo_todo \
-doSaf 1 -GL 1 \
-minInd 6 \
-nThreads 20 \
-out $SAF/littoral-pelagic/pelagic_pure_IDs_demo  &>> $SAF/littoral-pelagic/pelagic_pure_IDs_demo.log

## SFS per-morph (set minInd to 90% of each pop's sample size)
realSFS $SAF/littoral-pelagic/littoral_pure_IDs_demo.saf.idx >$SFS/littoral-pelagic/littoral_pure_IDs_demo.sfs
realSFS $SAF/littoral-pelagic/pelagic_pure_IDs_demo.saf.idx >$SFS/littoral-pelagic/pelagic_pure_IDs_demo.sfs

## 2D-SFS
realSFS \
$SAF/littoral-pelagic/littoral_pure_IDs_demo.saf.idx \
$SAF/littoral-pelagic/pelagic_pure_IDs_demo.saf.idx \
> $SFS/littoral-pelagic/littoral-pelagic_pure_IDs_demo.sfs

# Bootstrap 2D-SFS
realSFS -bootstrap 100 -resample_chr 1 \
$SAF/littoral-pelagic/littoral_pure_IDs_demo.saf.idx \
$SAF/littoral-pelagic/pelagic_pure_IDs_demo.saf.idx \
> $SFS/littoral-pelagic/littoral-pelagic_pure_IDs_demo_bootstrapped.sfs



###############################
###--- Admixture mapping ---###
###############################

### 5 | Admixture mapping based on 41 admixed samples using GEMMA

angsd -bam $BAM_WS/Admix_41IDs.bamlist \
-anc $GENOME_WS \
-doMaf 1 -doGlf 2 -doMajorMinor 1 -doGeno 2 -doPost 1 -GL 1 -doCounts \
$FILTERS_GWAS -minIndDepth 2 -nThreads 20 -SNP_pval 1e-6 \
-minMaf 0.1 \
-minInd 30 \
-postCutoff 0.95 -doVcf 1 \
-sites $SITES_WS/Admixed_41IDs.sites \
-out $GENO_WS/Admixed_41IDs_VCF &>> $GENO_WS/Admixed_41IDs_VCF.log 

## imputation using BEAGLE program
beagle.java
cd /home/DATA10/BOHAO/Whitefish/GWAS/Imputation

java -jar beagle.18May20.d20.jar gt=$GENO_WS/Admixed_41IDs_VCF.vcf.gz \
impute=true out=Admixed_41IDs_VCF_imputed

## BED format
cd /scratch/project_2000466/Whitefish/GWAS/GEMMA
module load biokit plink

vcftools --vcf Admixed_41IDs_VCF_imputed.vcf --plink --out Admixed_41IDs
plink --file Admixed_41IDs --make-bed --out Admixed_41IDs

## Perform GEMMA
GEMMA=/projappl/project_2000466/PROGRAMS/GEMMA/gemma

## generate relatedness matrix
$GEMMA -bfile Admixed_41IDs -p Admixed_41IDs.recode.pheno.txt -gk -o Admixed_41IDs

## run GEMMA for 22 traits using univariate linear models
for a in `seq 1 22`
do $GEMMA -bfile Admixed_41IDs -k output/Admixed_41IDs.cXX.txt -lmm 4 -n $a -o Admixed_41IDs_MAF01_imputed_$a
done

## Apply Multiple Testing Corrections over p-values derived from GEMMA
## In R
nSNPs <- 10367 # the number of genome-wide SNPs in the analyses
nTraits <- 22 # 22 traits under admixture mapping
p_multi_test_correction <- p.adjust(p_wald_Trait, method = "fdr", n = nTraits * nSNPs) 


### 6 |  Admixture mapping using LDna-EMMAX 
## The analysis are performed with R-script "EMMAX-LDna_Scripts(Zenodo).R", deposited in Dryad as a separate file

# python $PCANGSD -beagle $GENO_WS/Geno_littoral_pelagic_profundal_filtered.beagle.gz -o Geno_littoral_pelagic_profundal_filtered_maf01 -threads 20 -sites_save -minMaf 0.1












