# This assumes 97 chromosomes! It will probably work for fewer, but not much more

ANGSD=$1
REF=$2
BAMS=$3
PROJ=$4
SAMPS=$5
POPS=$6

# convert angsd maf file to tab separated chromosome\tposition per line
python3 /home/asharo/scripts/prepPositionFilesForPileupCaller.py $ANGSD $PROJ

# create pileup file from bams
samtools mpileup -R -B -q20 -Q20 -f $REF -b $BAMS -l $PROJ.tsv > $PROJ.txt

# change chromsome names to numbers
cat $PROJ.tsv | cut -f 1 | uniq > $PROJ.allChrmsInOrder.txt
python3 /home/asharo/scripts/renameChrmsForPileupCaller.py $PROJ

# run pileupCaller  CAN ADD "--skipTransitions" here for ancient samples
pileupCaller --randomHaploid --skipTransitions --sampleNameFile $SAMPS --samplePopName temp -f $PROJ.numberedChr.snp -e $PROJ.pileupCaller < $PROJ.numberedChr.txt

# add reference as a population, remove chr 5
python3 /home/asharo/scripts/RetainChr5InversionAddRef.py $PROJ

# optionally prune variants

# first, create a conversion file to ped plink format

echo indivname: $PROJ.pileupCaller.ref.pops.ind > par.EIGENSTRAT.PED
echo snpname: $PROJ.pileupCaller.ref.snp >> par.EIGENSTRAT.PED
echo genotypename: $PROJ.pileupCaller.ref.geno >> par.EIGENSTRAT.PED
echo outputformat: PED >> par.EIGENSTRAT.PED
echo genotypeoutname: $PROJ.ped >> par.EIGENSTRAT.PED
echo snpoutname: $PROJ.map >> par.EIGENSTRAT.PED
echo indivoutname: $PROJ.pedind >> par.EIGENSTRAT.PED

# create indiv file
awk '{FS=OFS="\t"}FNR==NR{a[NR]=$1;next}{$3=a[FNR]}1' $POPS $PROJ.pileupCaller.ref.ind > $PROJ.pileupCaller.ref.pops.ind

# convert to plink format

/redser4/personal/jonas/src/AdmixTools/bin/convertf -p ./par.EIGENSTRAT.PED

# run plink
plink --file $PROJ --indep-pairwise 50 5 0.5 --allow-extra-chr --out plink --chr-set 95

# extract just the unpruned vars
plink --file $PROJ --extract plink.prune.in --allow-extra-chr --chr-set 95 --recode --out plink.pruned

# create a qpDstat file

cat $POPS | sort | uniq > $PROJ.qpDstat

# create f3 parfile
echo indivname: $PROJ.pileupCaller.ref.pops.ind > $PROJ.f3.parfile
#echo snpname: $PROJ.pileupCaller.snp >> $PROJ.parfile
echo snpname: plink.pruned.map >> $PROJ.f3.parfile
# echo genotypename: $PROJ.pileupCaller.geno >> $PROJ.parfile
echo genotypename: plink.pruned.ped >> $PROJ.f3.parfile
# echo evecoutname: $PROJ.evecout >> $PROJ.parfile
# echo evaloutname: $PROJ.evalout >> $PROJ.parfile
echo numchrom: 32 >> $PROJ.f3.parfile
echo outgroupmode: YES >> $PROJ.f3.parfile
#echo outliersigmathresh: 8 >> $PROJ.parfile
echo popfilename: $PROJ.qp3pop >> $PROJ.f3.parfile

# create d-stat parfile
echo indivname: $PROJ.pileupCaller.ref.pops.ind > $PROJ.d.parfile
#echo snpname: $PROJ.pileupCaller.snp >> $PROJ.parfile
echo snpname: plink.pruned.map >> $PROJ.d.parfile
# echo genotypename: $PROJ.pileupCaller.geno >> $PROJ.parfile
echo genotypename: plink.pruned.ped >> $PROJ.d.parfile
# echo evecoutname: $PROJ.evecout >> $PROJ.parfile
# echo evaloutname: $PROJ.evalout >> $PROJ.parfile
echo poplistname: $PROJ.qpDstat >> $PROJ.d.parfile
echo numchrom: 32 >> $PROJ.d.parfile
#echo outgroupmode: YES >> $PROJ.f3.parfile
#echo outliersigmathresh: 8 >> $PROJ.parfile
#echo popfilename: $PROJ.qp3Pop >> $PROJ.parfile

# create smartpca parfile
echo indivname: $PROJ.pileupCaller.ref.pops.ind > $PROJ.smart.parfile
echo snpname: plink.pruned.map >> $PROJ.smart.parfile
echo genotypename: plink.pruned.ped >> $PROJ.smart.parfile
echo numchrom: 32 >> $PROJ.smart.parfile
echo evecoutname: $PROJ.evecout >> $PROJ.smart.parfile
echo evaloutname: $PROJ.evalout >> $PROJ.smart.parfile
echo outliersigmathresh: 12 >> $PROJ.smart.parfile

# run smartpca
export LD_LIBRARY_PATH=
smartpca -p $PROJ.smart.parfile > $PROJ.smart.logfile

# create qp3 file
# python3 /home/asharo/scripts/makeAllPopFiles.py $POPS $PROJ.qp3pop

# run f3
# /redser4/personal/jonas/src/AdmixTools/bin/qp3Pop -p $PROJ.f3.parfile > $PROJ.f3.logfile

# run d-stat
# /redser4/personal/jonas/src/AdmixTools/bin/qpDstat -p $PROJ.d.parfile > $PROJ.d.logfile

