/
PCA_vs_IBD.R
98 lines (89 loc) · 6.03 KB
/
PCA_vs_IBD.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
# Written 2013 by Peter Ralph and Graham Coop
#
# contact: petrel.harp@gmail.com
#
# To the extent possible under law, the author(s) have dedicated all copyright and related and neighboring rights to this software to the public domain worldwide. This software is distributed without any warranty.
#
# You should have received a copy of the CC0 Public Domain Dedication along with this software. If not, see <http://creativecommons.org/publicdomain/zero/1.0/>.
#
#
source(ibd-blocks-fns.R")
load("eda-data-fine.Rdata")
library("colorspace")
#####
pdf(file="Italy_PCA_vs_IBD.pdf")
plist <- c("France", "Italy", "Greece", "Turkey", "Cyprus")
pcols <- countrycols; pcols[plist] <- rainbow_hcl(length(plist),l=60,c=70)
layout(matrix(1:4,nrow=2))
par(mar=c(4,4,1,1))
#xlim=c(150,700), ylim=c(55,250),
with( subset(indivinfo,COUNTRY_SELF%in%plist), {
plot( United.Kingdom.blocks ~ PC1, pch=21,ylim=c(55,250),xlim=c(-0.07,.01), xlab="PC1", ylab="# blocks with UK", cex=1.5,
subset=(COUNTRY_SELF=="Italy"), bg=adjustcolor(pcols[COUNTRY_SELF],.75), col=NA, )
points( United.Kingdom.blocks ~ PC1, cex=1.5, subset=(COUNTRY_SELF!="Italy"), bg=adjustcolor(pcols[COUNTRY_SELF],.75),
pch=ifelse(COUNTRY_SELF%in%c("Switzerland","France"),22,23), col=adjustcolor("black",.25) )
} )
with( subset(indivinfo,COUNTRY_SELF%in%plist), {
plot( Swiss.French.blocks ~ PC1, pch=21,ylim=c(150,700),xlim=c(-0.07,.01), xlab="PC1", ylab="# blocks with CHf", cex=1.5,
subset=(COUNTRY_SELF=="Italy"), bg=adjustcolor(pcols[COUNTRY_SELF],.75), col=NA, main="" )
points( Swiss.French.blocks ~ PC1, cex=1.5, subset=(COUNTRY_SELF!="Italy"), bg=adjustcolor(pcols[COUNTRY_SELF],.75),
pch=ifelse(COUNTRY_SELF%in%c("Switzerland","France"),22,23), col=adjustcolor("black",.25) )
} )
with( subset(indivinfo,COUNTRY_SELF%in%plist), {
plot( United.Kingdom.blocks ~ PC2, pch=21,xlim=c(-0.03,0.05),ylim=c(55,250), xlab="PC2", ylab="# blocks with UK", cex=1.5,
subset=(COUNTRY_SELF=="Italy"), bg=adjustcolor(pcols[COUNTRY_SELF],.75), col=NA, main="" )
points( United.Kingdom.blocks ~ PC2, cex=1.5, subset=(COUNTRY_SELF!="Italy"), bg=adjustcolor(pcols[COUNTRY_SELF],.75),
pch=ifelse(COUNTRY_SELF%in%c("Switzerland","France"),22,23), col=adjustcolor("black",.25) )
} )
legend("topright", legend=countryabbrevs[plist], pch=c(22,21,23,23,23), pt.cex=1.5, bg="white", pt.bg=pcols[plist], col=ifelse(plist=="Italy",NA,"black") )
with( subset(indivinfo,COUNTRY_SELF%in%plist), {
plot( Swiss.French.blocks ~ PC2, pch=21,xlim=c(-0.03,0.05),ylim=c(150,700),, xlab="PC2", ylab="# blocks with CHf", cex=1.5,
subset=(COUNTRY_SELF=="Italy"), bg=adjustcolor(pcols[COUNTRY_SELF],.75), col=NA, main="" )
points( Swiss.French.blocks ~ PC2, cex=1.5, subset=(COUNTRY_SELF!="Italy"), bg=adjustcolor(pcols[COUNTRY_SELF],.75),
pch=ifelse(COUNTRY_SELF%in%c("Switzerland","France"),22,23), col=adjustcolor("black",.25) )
} )
dev.off()
####UK
pdf(file="UK_PCA_vs_IBD.pdf")
layout(matrix(1:4,nrow=2))
par(mar=c(4,4,1,1))
plist <- c("Germany","United Kingdom","Ireland")
pcols <- countrycols; pcols[plist] <- rainbow_hcl(length(plist),l=60,c=70)
with( subset(indivinfo,COUNTRY_SELF%in%plist), {
plot( Germany.blocks ~ PC1, pch=21, xlab="PC1", ylab="# blocks with Germany", bg=adjustcolor(pcols[COUNTRY_SELF],.75), cex=1.5, col=NA,
subset=(COUNTRY_SELF=="United Kingdom"), main="")
points( Germany.blocks ~ PC1, pch=ifelse(COUNTRY_SELF%in%c("Germany","Poland"),22,23), bg=adjustcolor(pcols[COUNTRY_SELF],.75),
col=adjustcolor("black",.25), cex=1.5, subset=(COUNTRY_SELF!="United Kingdom") )
} )
legend("topleft",pch=c(22,21,23),pt.bg=pcols[plist],legend=countryabbrevs[plist],pt.cex=1.5,col=ifelse(plist=="United Kingdom",NA,"black") )
with( subset(indivinfo,COUNTRY_SELF%in%plist), {
plot( Ireland.blocks ~ PC1, pch=21, xlab="PC1", ylab="# blocks with Ireland", bg=adjustcolor(pcols[COUNTRY_SELF],.75), cex=1.5, col=NA,
subset=(COUNTRY_SELF=="United Kingdom"), main="")
points( Ireland.blocks ~ PC1, pch=ifelse(COUNTRY_SELF%in%c("Germany","Poland"),22,23), bg=adjustcolor(pcols[COUNTRY_SELF],.75),
col=adjustcolor("black",.25), cex=1.5, subset=(COUNTRY_SELF!="United Kingdom") )
} )
with( subset(indivinfo,COUNTRY_SELF%in%plist), {
plot( Germany.blocks ~ PC2, pch=21, xlab="PC2", ylab="# blocks with Germany", bg=adjustcolor(pcols[COUNTRY_SELF],.75), cex=1.5, col=NA,
subset=(COUNTRY_SELF=="United Kingdom"), main="")
points( Germany.blocks ~ PC2, pch=ifelse(COUNTRY_SELF%in%c("Germany","Poland"),22,23), bg=adjustcolor(pcols[COUNTRY_SELF],.75),
col=adjustcolor("black",.25), cex=1.5, subset=(COUNTRY_SELF!="United Kingdom") )
} )
with( subset(indivinfo,COUNTRY_SELF%in%plist), {
plot( Ireland.blocks ~ PC2, pch=21, xlab="PC2", ylab="# blocks with Ireland", bg=adjustcolor(pcols[COUNTRY_SELF],.75), cex=1.5, col=NA,
subset=(COUNTRY_SELF=="United Kingdom"), main="")
points( Ireland.blocks ~ PC2, pch=ifelse(COUNTRY_SELF%in%c("Germany","Poland"),22,23), bg=adjustcolor(pcols[COUNTRY_SELF],.75),
col=adjustcolor("black",.25), cex=1.5, subset=(COUNTRY_SELF!="United Kingdom") )
} )
dev.off()
country.symbol <- rep(c(21,22,23,24,25),8)
names(country.symbol) <- levels(indivinfo$COUNTRY_SELF)
nsamples <- table(indivinfo$COUNTRY_SELF)
pdf("pca-map.pdf",width=6.5,height=6.5,pointsize=10)
par(mar=c(5,4,1,1)+.1)
pcols <- countrycols;
with(indivinfo,plot(PC1,PC2,bg=adjustcolor(pcols[COUNTRY_SELF],0.5),pch=country.symbol[COUNTRY_SELF],col=NA,xlim=c(-0.07,0.05)))
country.means <- do.call(rbind, with(indivinfo, tapply( 1:nrow(indivinfo), COUNTRY_SELF, function (k) { colMeans( indivinfo[k,c("PC1","PC2")] ) } ) ) )
text( country.means, labels=countryabbrevs )
points( country.means, pch=21, col=NA, bg=adjustcolor(countrycols,.25), cex=4 )
legend( "topright",pt.bg=countrycols,legend=countryabbrevs,pch=country.symbol,col=NA,cex=8/12,pt.cex=0.9)
dev.off()