# first we read in two files: the EIGENSTRAT output file, and a file with three columns: family ID, person ID, province (of current living address)

PCs_NL <- read.delim(file="dutch.R.evec",skip=1,head=FALSE,sep="")
province <- read.delim(file="NL_provinces.txt",skip=1,head=FALSE,sep="")

# then we merge the two matrices:
PCs_NL <- merge(PCs_NL,province,by.x=c("V1","V2"),by.y=c("V1","V2"))

	## plot ##

# choose colors for plot
cls=as.character(PCs_NL[,14])
cls[cls=="1"]="#FF82AB"	# Noord-Holland	(525)
cls[cls=="2"]="#BF3EFF"	# Zuid-Holland 	(96)
cls[cls=="3"]="#87CEEB"	# Zeeland 	(589)
cls[cls=="4"]="#FFB90F"	# Utrecht	(76)
cls[cls=="5"]="#4876FF"	# Noord-Brabant (563)
cls[cls=="6"]="#000080"	# Limburg 	(490)
cls[cls=="7"]="#FF7F24"	# Gelderland	(53)
cls[cls=="8"]="#8B0000"	# Overijssel	(556)
cls[cls=="9"]="#FF0000"	# Flevoland	(255)
cls[cls=="10"]="#54FF9F"	# Friesland	(575)
cls[cls=="11"]="#32CD32"	# Drenthe	(448)
cls[cls=="12"]="#6E8B3D"	# Groningen	(89)
 
# plot
pdf(file="NL_PC1_PC2.pdf",width=12.27,height=8.27)
par(xpd=NA, mfrow=c(1,1), oma=c(1,1,1,20))
plot(PCs_NL[,4],PCs_NL[,3],col=cls,type="p",pch=20,cex=2.5,bty="l",xlab="PC 2",ylab="PC 1")

legend(0.26,0.15,bty="n",
legend=c("Friesland","Groningen","Drenthe","Noord-Holland","Zuid-Holland","Utrecht","Flevoland",
"Gelderland","Overijssel","Zeeland","Noord-Brabant","Limburg"),
col=c("#54FF9F","#6E8B3D","#32CD32","#FF82AB","#BF3EFF","#FFB90F","#FF0000","#FF7F24","#8B0000",
"#87CEEB","#4876FF","#000080"), cex=1.4, lty = c(-1), pch = c(20),lwd=c(-1))

dev.off()