species location longitude latitude ancmean_z ancmean_auto
1 hybrid BCR -110.5530 43.00758 0.5240058 0.5057198
2 hybrid BIC -108.4012 45.16170 0.7421017 0.5023385
3 hybrid BLD -109.7161 43.52245 0.7123656 0.5104222
4 hybrid BNP -110.7212 44.93367 0.8183687 0.7400420
5 hybrid BTB -110.6820 43.63816 0.7240039 0.5661119
10 hybrid FRC -109.5776 43.72059 0.7196407 0.5175831
15 hybrid PIN -109.9762 43.74010 0.7060079 0.5621597
16 hybrid PSP -110.8493 42.74679 0.7100711 0.5327102
17 hybrid RDL -110.5465 44.36166 0.7260354 0.6058564
18 hybrid RNV -110.8847 43.59568 0.7210690 0.5687323
#read in the locations details table
poplocs<-read.table("../../../entropy/complete/mcmc/admixpoplocs.txt", header=T)
poplocs$species[4]<-"hybrid"
#subset JHL pops
hyblocs<-poplocs[poplocs$species == "hybrid",]
#get anc mean for each pop across all chromosomes
bcr.m<-mean(bcr3$mean)
bic.m<-mean(bic3$mean)
bld.m<-mean(bld3$mean)
btb.m<-mean(btb3$mean)
frc.m<-mean(frc3$mean)
pin.m<-mean(pin3$mean)
psp.m<-mean(psp3$mean)
rdl.m<-mean(rdl3$mean)
rnv.m<-mean(rnv3$mean)
## do z and autosome separately
bcr_z<-bcr3_mp_f[bcr3_mp_f$scaffold == 1631,]
bcr_auto<-bcr3_mp_f[!bcr3_mp_f$scaffold == 1631,]
dim(bcr_z)
dim(bcr_auto)
bcr_z_m<-mean(bcr_z$mean)
bcr_auto_m<-mean(bcr_auto$mean)
bic_z<-bic3_mp_f[bic3_mp_f$scaffold == 1631,]
bic_auto<-bic3_mp_f[!bic3_mp_f$scaffold == 1631,]
dim(bic_z)
dim(bic_auto)
bic_z_m<-mean(bic_z$mean)
bic_auto_m<-mean(bic_auto$mean)
bld_z<-bld3_mp_f[bld3_mp_f$scaffold == 1631,]
bld_auto<-bld3_mp_f[!bld3_mp_f$scaffold == 1631,]
dim(bld_z)
dim(bld_auto)
bld_z_m<-mean(bld_z$mean)
bld_auto_m<-mean(bld_auto$mean)
bnp_z<-bnp3_mp_f[bnp3_mp_f$scaffold == 1631,]
bnp_auto<-bnp3_mp_f[!bnp3_mp_f$scaffold == 1631,]
dim(bnp_z)
dim(bnp_auto)
bnp_z_m<-mean(bnp_z$mean)
bnp_auto_m<-mean(bnp_auto$mean)
btb_z<-btb3_mp_f[btb3_mp_f$scaffold == 1631,]
btb_auto<-btb3_mp_f[!btb3_mp_f$scaffold == 1631,]
dim(btb_z)
dim(btb_auto)
btb_z_m<-mean(btb_z$mean)
btb_auto_m<-mean(btb_auto$mean)
frc_z<-frc3_mp_f[frc3_mp_f$scaffold == 1631,]
frc_auto<-frc3_mp_f[!frc3_mp_f$scaffold == 1631,]
dim(frc_z)
dim(frc_auto)
frc_z_m<-mean(frc_z$mean)
frc_auto_m<-mean(frc_auto$mean)
pin_z<-pin3_mp_f[pin3_mp_f$scaffold == 1631,]
pin_auto<-pin3_mp_f[!pin3_mp_f$scaffold == 1631,]
dim(pin_z)
dim(pin_auto)
pin_z_m<-mean(pin_z$mean)
pin_auto_m<-mean(pin_auto$mean)
psp_z<-psp3_mp_f[psp3_mp_f$scaffold == 1631,]
psp_auto<-psp3_mp_f[!psp3_mp_f$scaffold == 1631,]
dim(psp_z)
dim(psp_auto)
psp_z_m<-mean(psp_z$mean)
psp_auto_m<-mean(psp_auto$mean)
rdl_z<-rdl3_mp_f[rdl3_mp_f$scaffold == 1631,]
rdl_auto<-rdl3_mp_f[!rdl3_mp_f$scaffold == 1631,]
dim(rdl_z)
dim(rdl_auto)
rdl_z_m<-mean(rdl_z$mean)
rdl_auto_m<-mean(rdl_auto$mean)
rnv_z<-rnv3_mp_f[rnv3_mp_f$scaffold == 1631,]
rnv_auto<-rnv3_mp_f[!rnv3_mp_f$scaffold == 1631,]
dim(rnv_z)
dim(rnv_auto)
rnv_z_m<-mean(rnv_z$mean)
rnv_auto_m<-mean(rnv_auto$mean)
#add this info to the locations table
ancmean_z<-c(bcr_z_m,bic_z_m,bld_z_m,bnp_z_m,btb_z_m,frc_z_m,pin_z_m,psp_z_m,rdl_z_m,rnv_z_m)
ancmean_auto<-c(bcr_auto_m,bic_auto_m,bld_auto_m,bnp_auto_m,btb_auto_m,frc_auto_m,pin_auto_m,psp_auto_m,rdl_auto_m,rnv_auto_m)
hyblocs<-cbind(hyblocs,ancmean_z,ancmean_auto)
ancmeanlg2<-1 - hyblocs$ancmeanlg
total<-hyblocs$ancmeanlg + hyblocs$ancmeanlg2
radius<-rep(0.03,9)
hyblocs<-cbind(hyblocs, ancmeanlg2, radius)
library(ggplot2)
library(ggmap)
library(maps)
library(mapdata)
library(scatterpie)
family<-as.factor(hyblocs$location)
usa <- map_data("usa")
state_dat<-map_data("state")
map_dat<-rbind(state_dat,usa)
#######################
library(gridBase)
library(grid)
#pdf("./mappcabarplot.pdf",width=10,height=10)
family<-as.factor(pops$V1)
usa <- map_data("usa")
state_dat<-map_data("state")
map_dat<-rbind(state_dat,usa)
### Sam is here for plotting ###
anc2<-1-hyblocs$ancmean_z
anc3<-1-hyblocs$ancmean_auto
hyblocs_z<-cbind(hyblocs, anc2)
hyblocs_auto<-cbind(hyblocs_auto,anc3)
hyblocs_z
colnames(hyblocs_z)<-c("hyblocs$location","longitude","latitude","L.idas","L.melissa")
hyblocs_auto<-data.frame(hyblocs_auto)
hyblocs_auto<-hyblocs_auto[,-c(1,2)]
hyblocs_auto<-cbind(hyblocs$location, hyblocs_auto)
colnames(hyblocs_auto)<-c("hyblocs$location","longitude","latitude","L.idas","L.melissa")
par(mfrow=c(1,2))
#z
p1<-ggplot() +
geom_polygon(data=map_dat, aes(x=long, y=lat, group=group), fill="white", color="black") +
ggtitle("C) Z chromosome") +
geom_scatterpie(aes(x=longitude, y=latitude, r=0.16), data=hyblocs_z, cols=c("L.idas","L.melissa"), color="black", alpha=0.9) +
scale_fill_manual(
labels = c("L. idas", "L. melissa"),
values = c("L.idas" = "gold",
"L.melissa" = "royalblue")) +
geom_text(data = hyblocs_z,
aes(longitude, latitude, label = hyblocs$location, fontface = 2),
nudge_x = 0.4,
size = 4,
color = "black") +
ylab("Latitude\n") + xlab("Longitude") +
coord_equal(xlim=c(-112,-107),ylim=c(41.5,46)) +
theme(
panel.background = element_rect(fill="lightsteelblue2"),
panel.grid.minor = element_line(colour="grey90", size=0.2),
panel.grid.major = element_line(colour="grey90", size=0.2),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.margin=unit(c(1,0.5,1,1),'cm'),
legend.position='none') +
theme(
axis.title.y = element_text(size = 20, vjust=1.3,hjust=0.5)) +
theme(axis.text.x = element_blank(),axis.text.y = element_blank(),axis.ticks=element_blank(),axis.title.x = element_text(size = 20,vjust=-1.3,hjust=0.5), legend.title = element_blank())
#auto
p2<-ggplot() +
geom_polygon(data=map_dat, aes(x=long, y=lat, group=group), fill="white", color="black") +
ggtitle("(D) Autosomes") +
geom_scatterpie(aes(x=longitude, y=latitude, r=0.16), data=hyblocs_auto, cols=c("L.idas","L.melissa"), color="black", alpha=0.9) +
scale_fill_manual(
labels = c("L. idas", "L. melissa"),
values = c("L.idas" = "gold",
"L.melissa" = "royalblue")) +
geom_text(data = hyblocs_auto,
aes(longitude, latitude, label = hyblocs$location, fontface = 2),
nudge_x = 0.4,
size = 4,
color = "black") +
ylab("Latitude\n") + xlab("Longitude") +
coord_equal(xlim=c(-112,-107),ylim=c(41.5,46)) +
theme(
panel.background = element_rect(fill="lightsteelblue2"),
panel.grid.minor = element_line(colour="grey90", size=0.2),
panel.grid.major = element_line(colour="grey90", size=0.2),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.margin=unit(c(1,0.5,1,1),'cm'),
legend.position=c(0.8,0.2), legend.text=element_text(size=6),legend.key.size = unit(0.5,"line")) +
theme(axis.title.y = element_blank(),axis.text.x = element_blank(),axis.text.y = element_blank(),axis.ticks=element_blank(),axis.title.x = element_text(size = 15,vjust=-1.3,hjust=0.5), legend.title = element_blank())
grid.arrange(p1, p2, ncol=2)
######## sam is here
layout(matrix(c(1,1,2,2,3,4), nrow=3, ncol=2, byrow = TRUE))
#bld
par(mar=c(6,7,5,2))
o_bld3_mp_f<-bld3_mp_f[order(as.integer(bld3_mp_f$linkage), decreasing=FALSE),]
o_bld3_mp_f$linkage<-factor(o_bld3_mp_f$linkage, levels=c(1:22, "Z"))
cols<-rep("lemonchiffon4","22")
cols<-c(cols,"darkred")
fills<-rep("lemonchiffon",22)
fills<-c(fills,"darksalmon")
boxplot(o_bld3_mp_f$mean ~ o_bld3_mp_f$linkage,col=fills,border=cols,xlab = "", ylab=expression(paste(italic("L.idas "),"ancestry")), xlab = "",cex.lab=2.5,cex.main=2,cex.axis=1.5,xaxt='n',yaxt='n',ylim=c(0,1))
axis(side=2,at=c(0,0.5,1),cex.axis=1.5,labels=c(0,0.5,1))
axis(side=1, at=mticks,cex.axis=1.3, labels=c("1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","16","17","18","19","20","21","22","Z"))
title(main="(A) Bald Mountain, WY",cex.main = 2, adj=0)
#pin
par(mar=c(6,7,5,2))
#boxplot
o_pin3_mp_f<-pin3_mp_f[order(as.integer(pin3_mp_f$linkage), decreasing=FALSE),]
o_pin3_mp_f$linkage<-factor(o_pin3_mp_f$linkage, levels=c(1:22, "Z"))
cols<-rep("lemonchiffon4","22")
cols<-c(cols,"darkred")
fills<-rep("lemonchiffon",22)
fills<-c(fills,"darksalmon")
boxplot(o_pin3_mp_f$mean ~ o_pin3_mp_f$linkage,col=fills,border=cols,xlab = "", ylab=expression(paste(italic("L.idas "),"ancestry")), xlab = "Linkage group",cex.lab=2.5,cex.main=2,cex.axis=1.5,xaxt='n',yaxt='n',ylim=c(0,1))
axis(side=2,at=c(0,0.5,1),cex.axis=1.5,labels=c(0,0.5,1))
axis(side=1, at=mticks,cex.axis=1.3, labels=c("1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","16","17","18","19","20","21","22","Z"))
title(main="(B) Pinnacle, MT",cex.main = 2, adj=0)
#ggplot
#z
p1<-ggplot() +
geom_polygon(data=map_dat, aes(x=long, y=lat, group=group), fill="white", color="black") +
geom_scatterpie(aes(x=longitude, y=latitude, r=0.16), data=hyblocs_z, cols=c("L.idas","L.melissa"), color="black", alpha=0.9) +
scale_fill_manual(
labels = c("L. idas", "L. melissa"),
values = c("L.idas" = "gold",
"L.melissa" = "royalblue")) +
geom_text(data = hyblocs_z,
aes(longitude, latitude, label = hyblocs$location, fontface = 2),
nudge_x = 0.4,
size = 3,
color = "black") +
ylab("Latitude\n") + xlab("Longitude") +
coord_equal(xlim=c(-112,-107),ylim=c(42,45.8)) +
theme(
panel.background = element_rect(fill="lightsteelblue2"),
panel.grid.minor = element_line(colour="grey90", size=0.2),
panel.grid.major = element_line(colour="grey90", size=0.2),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.margin=unit(c(0.1,0.1,0.5,1),'cm'),
legend.position='none') +
theme(axis.text.x = element_blank(),axis.text.y = element_blank(),axis.ticks=element_blank(),
axis.title.y = element_text(size = 20, vjust=1.3,hjust=0.5)) +
theme(axis.title.x = element_text(size = 20,vjust=-1.3,hjust=0.5), legend.title = element_blank())
title(main="(C) Z chromosome",cex.main = 2, adj=0)
#auto
p2<-ggplot() +
geom_polygon(data=map_dat, aes(x=long, y=lat, group=group), fill="white", color="black") +
geom_scatterpie(aes(x=longitude, y=latitude, r=0.16), data=hyblocs_auto, cols=c("L.idas","L.melissa"), color="black", alpha=0.9) +
scale_fill_manual(
labels = c("L. idas", "L. melissa"),
values = c("L.idas" = "gold",
"L.melissa" = "royalblue")) +
geom_text(data = hyblocs_auto,
aes(longitude, latitude, label = hyblocs$location, fontface = 2),
nudge_x = 0.4,
size = 3,
color = "black") +
xlab("Longitude") +
coord_equal(xlim=c(-112,-107),ylim=c(42,45.8)) +
theme(
panel.background = element_rect(fill="lightsteelblue2"),
panel.grid.minor = element_line(colour="grey90", size=0.2),
panel.grid.major = element_line(colour="grey90", size=0.2),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.margin=unit(c(0.1,0.1,0.5,1),'cm'),
legend.position=c(0.8,0.2), legend.text=element_text(size=18,face= "italic"),legend.key.size = unit(0.8,"line"),legend.spacing.x = unit(0.3,'cm')) +
theme(axis.text.x = element_blank(),axis.text.y = element_blank(),axis.ticks=element_blank(),
axis.title.y = element_blank()) +
theme(axis.title.x = element_text(size = 20,vjust=-1.3,hjust=0.5), legend.title = element_blank())
print(p2, vp=vp.bottomleft)
vp.bottomright <- viewport(height=unit(.35, "npc"), width=unit(0.4, "npc"),
just=c("right","bottom"),
y=0, x=0.5)
print(p1, vp=vp.bottomright)
vp.bottomleft <- viewport(height=unit(.35, "npc"), width=unit(0.4, "npc"),
just=c("left","bottom"),
y=0, x=0.5)
print(p2, vp=vp.bottomleft)