#!/bin/bash
bcftools merge -l combine.txt -0 -Oz -o immigrant_merged.vcf.gz
bcftools index immigrant_merged.vcf.gz
vcftools --gzvcf immigrant_merged.vcf.gz --max-missing 0.8 --mac 3 --minQ 20 --recode --recode-INFO-all --maf 0.05 --hwe 0.001 --min-alleles 2 --max
-alleles 2 --out immigrant_M2_min_maf_hwe

cd part1_Ho_analysis
mv ../immigrant_M2_min_maf_hwe.recode.vcf ./

plink --vcf immigrant_M2_min_maf_hwe.recode.vcf --make-bed --out immigrant_out1 --allow-extra-chr

wc -l immigrant_out1.bim > snp_number1.txt

sed -e "s/NC_064243.1/1/g" -e "s/NC_064244.1/2/g" -e "s/NC_064245.1/3/g" -e "s/NC_064246.1/4/g" -e "s/NC_064247.1/5/g" -e "s/NC_064248.1/6/g" -e "s/
NC_064249.1/7/g" -e "s/NC_064250.1/8/g" -e "s/NC_064251.1/9/g" -e "s/NC_064252.1/10/g" -e "s/NC_064253.1/11/g" -e "s/NC_064254.1/12/g" -e "s/NC_0642
55.1/13/g" -e "s/NC_064256.1/14/g" -e "s/NC_064257.1/15/g" -e "s/NC_064258.1/16/g" -e "s/NC_064259.1/17/g" -e "s/NC_064260.1/18/g" -e "s/NC_064261.1
/19/g" -e "s/NC_064262.1/20/g" -e "s/NC_064263.1/21/g" -e "s/NC_064264.1/22/g" -e "s/NC_064265.1/23/g" -e "s/NC_064266.1/24/g" -e "s/NC_064267.1/25/
g" -e "s/NC_064268.1/26/g" -e "s/NC_064269.1/27/g" -e "s/NC_064270.1/28/g" -e "s/NC_064271.1/29/g" -e "s/NC_064272.1/30/g" -e "s/NC_064273.1/31/g" -
e "s/NC_064274.1/32/g" -e "s/NC_064275.1/33/g" -e "s/NC_064276.1/34/g" -e "s/NC_064277.1/35/g" -e "s/NC_064278.1/36/g" -e "s/NC_064279.1/37/g" -e "s
/NC_064280.1/38/g" -i immigrant_out1.bim

awk '{print $1 ,"\t", $2="snp"NR, "\t", $3 , "\t", $4 , "\t", $5 ,"\t", $6}' immigrant_out1.bim > immigrant_out1.txt

rm -rf immigrant_out1.bim | cp immigrant_out1.txt immigrant_out1.bim

rm -rf immigrant_out1.txt

awk '{ if($1 >=1 && $1 <= 38) print $2}' immigrant_out1.bim > snp_1_38.txt

plink --bfile immigrant_out1 --extract snp_1_38.txt --make-bed --out immigrant_out2 --allow-extra-chr --chr-set 38

rm -rf immigrant_out1*

wc -l immigrant_out2.bim > snp_number2.txt

plink --bfile immigrant_out2 --recode --out immigrant_out2 --allow-extra-chr --chr-set 38

plink --file immigrant_out2 --hardy --out immigrant_het_test --allow-extra-chr --chr-set 38

awk '{print $7 "\t" $8}' immigrant_het_test.hwe | sed '1d' > count_het.txt

awk '{print $1}' count_het.txt | datamash mean 1 sstdev 1 > immigrant_Ho.txt

awk '{print $2}' count_het.txt | datamash mean 1 sstdev 1 > immigrant_He.txt

plink --bfile immigrant_out2 --export vcf --out immigrant_1_38 --chr-set 38

vcftools --vcf immigrant_1_38.vcf --SNPdensity 100000 --out immigrant_SNPdensity

awk '{print $4}' immigrant_SNPdensity.snpden | sed '1d' > immigrant_SNPdensity.txt

awk '{print $1}' immigrant_SNPdensity.txt | datamash mean 1 sstdev 1 > immigrant_snpdensity.txt

vcftools --vcf immigrant_1_38.vcf --window-pi 100000 --out immigrant_pi

awk '{print $5}' immigrant_pi.windowed.pi | sed '1d' > A.txt

sort -n -r A.txt > B.txt

awk '{print $1}' A.txt | datamash mean 1 sstdev 1 > immigrant_mean_pi.txt

vcftools --vcf immigrant_1_38.vcf --TajimaD 10000  --out immigrant_TajimaD

rm -rf *log