# slide 16 #

# make the .par file:
echo "genotypename: dutch_1kG.ped" >> dutch_1kG.par
echo "snpname: dutch_1kG.map" >> dutch_1kG.par
echo "indivname: dutch_1kG.fam" >> dutch_1kG.par
echo "evecoutname: dutch_1kG.evec" >> dutch_1kG.par 
echo "evaloutname: dutch_1kG.eval" >> dutch_1kG.par
echo "numoutevec: 10" >> dutch_1kG.par 
echo "numoutlieriter: 0" >> dutch_1kG.par
echo "poplistname: poplist_1kG.txt" >> dutch_1kG.par
echo "snpweightoutname: dutch_1kG.snpweight" >> dutch_1kG.par

# make the poplistname file:
echo "4\n5\n6\n7\n8\n9\n10\n11\n12\n13" > poplist_1kG.txt

	# slide 17 #

# run EIGENSTRAT:
smartpca -p dutch_1kG.par > dutch_1kG.log

	# slide 18 #

# replace all colons with a space
sed 's/:/ /g' dutch_1kG.evec > dutch_1kG.R.evec

# run the R script
R CMD BATCH outliers.R

	# slide 22 #

# extract all 1000 Genomes individuals:
awk '$6>3{print $1,$2}' dutch_1kG.fam > 1kG.ids 

# make a file with all 1000 Genomes individuals and Dutch individuals with a non-European ancestry
cat outliers.txt 1kG.ids > remove_outliers.ids

# make plink file excluding these individuals
plink --bfile dutch_1kG --remove remove_outliers.ids --make-bed --out dutch 
plink --bfile dutch --recode --out dutch

	# slide 23 #

# make the .par file:
echo "genotypename: dutch.ped" >> dutch.par
echo "snpname: dutch.map" >> dutch.par
echo "indivname: dutch.fam" >> dutch.par
echo "evecoutname: dutch.evec" >> dutch.par 
echo "evaloutname: dutch.eval" >> dutch.par
echo "numoutevec: 10" >> dutch.par 
echo "numoutlieriter: 0" >> dutch.par
echo "poplistname: poplist_NL.txt" >> dutch.par
echo "snpweightoutname: dutch.snpweight" >> dutch.par

# make the poplistname file:
echo "3" > poplist_NL.txt

	# slide 24 #

# run EIGENSTRAT:
smartpca -p dutch.par > dutch.log

	# slide 25 #

# replace all colons with a space
sed 's/:/ /g' dutch.evec > dutch.R.evec

# run R script to make plot
R CMD BATCH plot_NL.R

	# slide 47 #

# compute ROHs for the dutch dataset
plink --bfile dutch --homozyg --homozyg-snp 65 --homozyg-window-het 0 --homozyg-density 200 --homozyg-gap 500 --homozyg-group --out dutch.ROHs