install.packages("ggplot2", dep=TRUE)
install.packages("doBy", dep=TRUE)
#http://had.co.nz/ggplot2

library(ggplot2)
library(doBy)

################Denisova ASIAN#######################


de <- which(DenisovaMat[,12]>0.3)

ASN1D <- PopulationMat[de,8]

ASN2D <- PopulationMat[-de,8]


l1 <- length(ASN1D)
l2 <- length(ASN2D)

df <- data.frame(Matching_Genome = factor( c(rep("Not Denisovan",l2),rep("Denisovan",l1)) ), proportion = c(ASN2D,ASN1D))

cdf <- summaryBy(data=df, proportion ~ Matching_Genome, FUN=mean, na.rm=TRUE)

ggplot(df, aes(x=proportion, fill=Matching_Genome)) + geom_density(alpha=.3) +
    geom_vline(data=cdf, aes(xintercept=proportion.mean,  colour=Matching_Genome),
               linetype="dashed", size=1) +
#    opts(title = expression("Asian Proportions in Neandertal Hapoltypes")) +
    scale_y_continuous(breaks=c(0,2,5,7)) +
    opts(legend.position=c(.5, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(-0.5,10)) +
    ylab("Density") +
    xlab("Asian proportion of IBD segment occurrence")





name="denisovanAsian"
pdf(paste(name,".pdf",sep=""),width=8)
ggplot(df, aes(x=proportion, fill=Matching_Genome)) + geom_density(alpha=.3) +
    geom_vline(data=cdf, aes(xintercept=proportion.mean,  colour=Matching_Genome),
               linetype="dashed", size=1) +
    scale_y_continuous(breaks=c(0,2,5,7)) +
    opts(legend.position=c(.5, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(-0.5,10)) +
    ylab("Density") +
    xlab("Asian proportion of IBD segment occurrence")

dev.off()

xfig(paste(name,".fig",sep=""),width=8)
ggplot(df, aes(x=proportion, fill=Matching_Genome)) + geom_density(alpha=.3) +
    geom_vline(data=cdf, aes(xintercept=proportion.mean,  colour=Matching_Genome),
               linetype="dashed", size=1) +
    scale_y_continuous(breaks=c(0,2,5,7)) +
    opts(legend.position=c(.5, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(-0.5,10)) +
    ylab("Density") +
    xlab("Asian proportion of IBD segment occurrence")
dev.off()

################Denisova European#######################


EUR1D <- PopulationMat[de,9]

EUR2D <- PopulationMat[-de,9]

l1 <- length(EUR1D)
l2 <- length(EUR2D)

df <- data.frame(Matching_Genome = factor( c(rep("Not Denisovan",l2),rep("Denisovan",l1)) ), proportion = c(EUR2D,EUR1D))

cdf <- summaryBy(data=df, proportion ~ Matching_Genome, FUN=mean, na.rm=TRUE)

ggplot(df, aes(x=proportion, fill=Matching_Genome)) + geom_density(alpha=.3) +
    geom_vline(data=cdf, aes(xintercept=proportion.mean,  colour=Matching_Genome),
               linetype="dashed", size=1) +
    opts(legend.position=c(.5, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(-0.5,12)) +
    ylab("Density") +
    xlab("European proportion of IBD segment occurrence")



name="denisovanEuropean"
pdf(paste(name,".pdf",sep=""),width=8)
ggplot(df, aes(x=proportion, fill=Matching_Genome)) + geom_density(alpha=.3) +
    geom_vline(data=cdf, aes(xintercept=proportion.mean,  colour=Matching_Genome),
               linetype="dashed", size=1) +
    opts(legend.position=c(.5, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(-0.5,12)) +
    ylab("Density") +
    xlab("European proportion of IBD segment occurrence")
dev.off()

xfig(paste(name,".fig",sep=""),width=8)
ggplot(df, aes(x=proportion, fill=Matching_Genome)) + geom_density(alpha=.3) +
    geom_vline(data=cdf, aes(xintercept=proportion.mean,  colour=Matching_Genome),
               linetype="dashed", size=1) +
    opts(legend.position=c(.5, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(-0.5,12)) +
    ylab("Density") +
    xlab("European proportion of IBD segment occurrence")
dev.off()


################Denisova African#######################


AFR1D <- PopulationMat[de,6]

AFR2D <- PopulationMat[-de,6]

l1 <- length(AFR1D)
l2 <- length(AFR2D)

df <- data.frame(Matching_Genome = factor( c(rep("Not Denisovan",l2),rep("Denisovan",l1)) ), proportion = c(AFR2D,AFR1D))

cdf <- summaryBy(data=df, proportion ~ Matching_Genome, FUN=mean, na.rm=TRUE)

ggplot(df, aes(x=proportion, fill=Matching_Genome)) + geom_density(alpha=.3) +
    geom_vline(data=cdf, aes(xintercept=proportion.mean,  colour=Matching_Genome),
               linetype="dashed", size=1) +
    opts(legend.position=c(.5, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(-0.5,12)) +
    ylab("Density") +
    xlab("African proportion of IBD segment occurrence")



name="denisovanAfrican"
pdf(paste(name,".pdf",sep=""),width=8)
ggplot(df, aes(x=proportion, fill=Matching_Genome)) + geom_density(alpha=.3) +
    geom_vline(data=cdf, aes(xintercept=proportion.mean,  colour=Matching_Genome),
               linetype="dashed", size=1) +
    opts(legend.position=c(.5, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(-0.5,12)) +
    ylab("Density") +
    xlab("African proportion of IBD segment occurrence")
dev.off()

xfig(paste(name,".fig",sep=""),width=8)
ggplot(df, aes(x=proportion, fill=Matching_Genome)) + geom_density(alpha=.3) +
    geom_vline(data=cdf, aes(xintercept=proportion.mean,  colour=Matching_Genome),
               linetype="dashed", size=1) +
    opts(legend.position=c(.5, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(-0.5,12)) +
    ylab("Density") +
    xlab("African proportion of IBD segment occurrence")
dev.off()



################Neandertal ASIAN#######################


ne <- which(NeanderMat[,12]>0.3)

ASN1D <- PopulationMat[ne,8]

ASN2D <- PopulationMat[-ne,8]


l1 <- length(ASN1D)
l2 <- length(ASN2D)

df <- data.frame(Matching_Genome = factor( c(rep("Not Neandertal",l2),rep("Neandertal",l1)) ), proportion = c(ASN2D,ASN1D))

cdf <- summaryBy(data=df, proportion ~ Matching_Genome, FUN=mean, na.rm=TRUE)

ggplot(df, aes(x=proportion, fill=Matching_Genome)) + geom_density(alpha=.3) +
    geom_vline(data=cdf, aes(xintercept=proportion.mean,  colour=Matching_Genome),
               linetype="dashed", size=1) +
    opts(legend.position=c(.5, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(-0.5,10)) +
    ylab("Density") +
    xlab("Asian proportion of IBD segment occurrence")





name="neandertalAsian"
pdf(paste(name,".pdf",sep=""),width=8)
ggplot(df, aes(x=proportion, fill=Matching_Genome)) + geom_density(alpha=.3) +
    geom_vline(data=cdf, aes(xintercept=proportion.mean,  colour=Matching_Genome),
               linetype="dashed", size=1) +
    opts(legend.position=c(.5, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(-0.5,10)) +
    ylab("Density") +
    xlab("Asian proportion of IBD segment occurrence")
dev.off()

xfig(paste(name,".fig",sep=""),width=8)
ggplot(df, aes(x=proportion, fill=Matching_Genome)) + geom_density(alpha=.3) +
    geom_vline(data=cdf, aes(xintercept=proportion.mean,  colour=Matching_Genome),
               linetype="dashed", size=1) +
    opts(legend.position=c(.5, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(-0.5,10)) +
    ylab("Density") +
    xlab("Asian proportion of IBD segment occurrence")
dev.off()

################Neandertal European#######################


EUR1D <- PopulationMat[ne,9]

EUR2D <- PopulationMat[-ne,9]

l1 <- length(EUR1D)
l2 <- length(EUR2D)

df <- data.frame(Matching_Genome = factor(  c(rep("Not Neandertal",l2),rep("Neandertal",l1))  ), proportion = c(EUR2D,EUR1D))

cdf <- summaryBy(data=df, proportion ~ Matching_Genome, FUN=mean, na.rm=TRUE)

ggplot(df, aes(x=proportion, fill=Matching_Genome)) + geom_density(alpha=.3) +
    geom_vline(data=cdf, aes(xintercept=proportion.mean,  colour=Matching_Genome),
               linetype="dashed", size=1) +
    opts(legend.position=c(.5, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(-0.5,10)) +
    ylab("Density") +
    xlab("European proportion of IBD segment occurrence")



name="neandertalEuropean"
pdf(paste(name,".pdf",sep=""),width=8)
ggplot(df, aes(x=proportion, fill=Matching_Genome)) + geom_density(alpha=.3) +
    geom_vline(data=cdf, aes(xintercept=proportion.mean,  colour=Matching_Genome),
               linetype="dashed", size=1) +
    opts(legend.position=c(.5, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(-0.5,10)) +
    ylab("Density") +
    xlab("European proportion of IBD segment occurrence")
dev.off()

xfig(paste(name,".fig",sep=""),width=8)
ggplot(df, aes(x=proportion, fill=Matching_Genome)) + geom_density(alpha=.3) +
    geom_vline(data=cdf, aes(xintercept=proportion.mean,  colour=Matching_Genome),
               linetype="dashed", size=1) +
    opts(legend.position=c(.5, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(-0.5,10)) +
    ylab("Density") +
    xlab("European proportion of IBD segment occurrence")
dev.off()

################Neandertal African#######################


AFR1D <- PopulationMat[ne,6]

AFR2D <- PopulationMat[-ne,6]

l1 <- length(AFR1D)
l2 <- length(AFR2D)

df <- data.frame(Matching_Genome = factor(  c(rep("Not Neandertal",l2),rep("Neandertal",l1))  ), proportion = c(AFR2D,AFR1D))

cdf <- summaryBy(data=df, proportion ~ Matching_Genome, FUN=mean, na.rm=TRUE)

ggplot(df, aes(x=proportion, fill=Matching_Genome)) + geom_density(alpha=.3) +
    geom_vline(data=cdf, aes(xintercept=proportion.mean,  colour=Matching_Genome),
               linetype="dashed", size=1) +
    opts(legend.position=c(.5, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(-0.5,10)) +
    ylab("Density") +
    xlab("African proportion of IBD segment occurrence")



name="neandertalAfrican"
pdf(paste(name,".pdf",sep=""),width=8)
ggplot(df, aes(x=proportion, fill=Matching_Genome)) + geom_density(alpha=.3) +
    geom_vline(data=cdf, aes(xintercept=proportion.mean,  colour=Matching_Genome),
               linetype="dashed", size=1) +
    opts(legend.position=c(.5, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(-0.5,10)) +
    ylab("Density") +
    xlab("African proportion of IBD segment occurrence")
dev.off()

xfig(paste(name,".fig",sep=""),width=8)
ggplot(df, aes(x=proportion, fill=Matching_Genome)) + geom_density(alpha=.3) +
    geom_vline(data=cdf, aes(xintercept=proportion.mean,  colour=Matching_Genome),
               linetype="dashed", size=1) +
    opts(legend.position=c(.5, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(-0.5,10)) +
    ylab("Density") +
    xlab("African proportion of IBD segment occurrence")
dev.off()