load(file=paste("DenisovaMat.Rda",sep=""))
load(file=paste("NeanderMat.Rda",sep=""))
load(file=paste("AncestMat.Rda",sep=""))
load(file=paste("PopulationMat.Rda",sep=""))
format(colSums(DenisovaMat),scientific=FALSE)
format(colSums(NeanderMat),scientific=FALSE)
format(colSums(AncestMat),scientific=FALSE)

source("funct.R")

library(ggplot2)
library(doBy)

cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

summary(DenisovaMat[,21])
summary(DenisovaMat[,22])
summary(DenisovaMat[,34])
# paper densiova/neanderthal ->  >= 0.3
#       population >= 0.2

shorT <- which(DenisovaMat[,5]==9)
shorTL <- DenisovaMat[shorT,22]
summary(shorTL)


de1 <- which(DenisovaMat[,12]>0.15)
de2 <- which(DenisovaMat[,35]>0.3)
de <- intersect(de1,de2)

DE1 <- which(DenisovaMat[,12]>0.3)
DE2 <- which(DenisovaMat[,35]>0.6)
DE <- intersect(DE1,DE2)

ne1 <- which(NeanderMat[,12]>0.15)
ne2 <- which(NeanderMat[,35]>0.3)
ne <- intersect(ne1,ne2)

NE1 <- which(NeanderMat[,12]>0.3)
NE2 <- which(NeanderMat[,35]>0.6)
NE <- intersect(NE1,NE2)



dene <- intersect(de,ne)

de <- setdiff(de,dene)
ne <- setdiff(ne,dene)

DENE <- intersect(DE,NE)

DE <- setdiff(DE,DENE)
NE <- setdiff(NE,DENE)




AFR <- which(PopulationMat[,6]>0.2) #147574

ASN <- which(PopulationMat[,8]>0.2) #7472

EUR <- which(PopulationMat[,9]>0.2) #10688


AFRd <- which(PopulationMat[,6]>0.5) #141262

ASNd <- which(PopulationMat[,8]>0.5) #5161

EURd <- which(PopulationMat[,9]>0.5) #5915


neL <- NeanderMat[ne,34]

deL <- DenisovaMat[de,34]


NEL <- NeanderMat[NE,34]

DEL <- DenisovaMat[DE,34]

deneL <- DenisovaMat[dene,34]


DENEL <- DenisovaMat[dene,34]



LLL <- DenisovaMat[,22]
LLLZ <- which(LLL>10)
LLL <- LLL[LLLZ]

LLLAFR <- DenisovaMat[AFR,22]
LLLASN <- DenisovaMat[ASN,22]
LLLEUR <- DenisovaMat[EUR,22]

LLLAFRd <- DenisovaMat[AFRd,22]
LLLASNd <- DenisovaMat[ASNd,22]
LLLEURd <- DenisovaMat[EURd,22]



DEL <- DenisovaMat[DE,34]

DLLL <- density(LLL,bw=1400,from=0,to=60000)
DLLL1 <- density(LLL,bw=1400,from=0,to=5000)
DLLL2 <- density(LLL,bw=1400,from=15000,to=60000)

DLLLAFR <- density(LLLAFR,bw=1400,from=0,to=60000)
DLLLAFR1 <- density(LLLAFR,bw=1400,from=0,to=5000)
DLLLAFR2 <- density(LLLAFR,bw=1400,from=15000,to=60000)

DLLLEUR <- density(LLLEUR,bw=1400,from=0,to=60000)
DLLLEUR1 <- density(LLLEUR,bw=1400,from=0,to=5000)
DLLLEUR2 <- density(LLLEUR,bw=1400,from=15000,to=60000)

DLLLASN <- density(LLLASN,bw=1400,from=0,to=60000)
DLLLASN1 <- density(LLLASN,bw=1400,from=0,to=5000)
DLLLASN2 <- density(LLLASN,bw=1400,from=15000,to=60000)


DLLLAFRd <- density(LLLAFRd,bw=1400,from=0,to=60000)
DLLLAFR1d <- density(LLLAFRd,bw=1400,from=0,to=5000)
DLLLAFR2d <- density(LLLAFRd,bw=1400,from=15000,to=60000)

DLLLEURd <- density(LLLEURd,bw=1400,from=0,to=60000)
DLLLEUR1d <- density(LLLEURd,bw=1400,from=0,to=5000)
DLLLEUR2d <- density(LLLEURd,bw=1400,from=15000,to=60000)

DLLLASNd <- density(LLLASNd,bw=1400,from=0,to=60000)
DLLLASN1d <- density(LLLASNd,bw=1400,from=0,to=5000)
DLLLASN2d <- density(LLLASNd,bw=1400,from=15000,to=60000)



DDEL <- density(DEL,bw=1400,from=0,to=60000)
DEASN <- intersect(ASN,DE)
DEASNL <- DenisovaMat[DEASN,34]
DDEASNL <- density(DEASNL,bw=1400,from=0,to=60000)
DEAFR <- intersect(AFR,DE)
DEAFRL <- DenisovaMat[DEAFR,34]
DDEAFRL <- density(DEAFRL,bw=1400,from=0,to=60000)
DEEUR <- intersect(EUR,DE)
DEEURL <- DenisovaMat[DEEUR,34]
DDEEURL <- density(DEEURL,bw=1400,from=0,to=60000)


DEASNd <- intersect(ASNd,DE) # 392
DEASNLd <- DenisovaMat[DEASNd,34]
DDEASNLd <- density(DEASNLd,bw=1400,from=0,to=60000)
DEAFRd <- intersect(AFRd,DE)
DEAFRLd <- DenisovaMat[DEAFRd,34]
DDEAFRLd <- density(DEAFRLd,bw=1400,from=0,to=60000)
DEEURd <- intersect(EURd,DE) # 206
DEEURLd <- DenisovaMat[DEEURd,34]
DDEEURLd <- density(DEEURLd,bw=1400,from=0,to=60000)



DdeL <- density(deL,bw=1400,from=0,to=60000)
DdeL1 <- density(deL,bw=1400,from=0,to=5000)
DdeL2 <- density(deL,bw=1400,from=15000,to=60000)
deASN <- intersect(ASN,de)
deASNL <- DenisovaMat[deASN,34]
DdeASNL <- density(deASNL,bw=1400,from=0,to=60000)
deAFR <- intersect(AFR,de)
deAFRL <- DenisovaMat[deAFR,34]
DdeAFRL <- density(deAFRL,bw=1400,from=0,to=60000)
deEUR <- intersect(EUR,de)
deEURL <- DenisovaMat[deEUR,34]
DdeEURL <- density(deEURL,bw=1400,from=0,to=60000)


deASNd <- intersect(ASNd,de) # 1303
deASNLd <- DenisovaMat[deASNd,34]
DdeASNLd <- density(deASNLd,bw=1400,from=0,to=60000)
deAFRd <- intersect(AFRd,de)
deAFRLd <- DenisovaMat[deAFRd,34]
DdeAFRLd <- density(deAFRLd,bw=1400,from=0,to=60000)
deEURd <- intersect(EURd,de) #886
deEURLd <- DenisovaMat[deEURd,34]
DdeEURLd <- density(deEURLd,bw=1400,from=0,to=60000)

##################################################


LLLAFR <- NeanderMat[AFR,22]
LLLASN <- NeanderMat[ASN,22]
LLLEUR <- NeanderMat[EUR,22]

LLLAFRd <- NeanderMat[AFRd,22]
LLLASNd <- NeanderMat[ASNd,22]
LLLEURd <- NeanderMat[EURd,22]



DLLLAFR <- density(LLLAFR,bw=1400,from=0,to=60000)
DLLLAFR1 <- density(LLLAFR,bw=1400,from=0,to=5000)
DLLLAFR2 <- density(LLLAFR,bw=1400,from=15000,to=60000)

DLLLEUR <- density(LLLEUR,bw=1400,from=0,to=60000)
DLLLEUR1 <- density(LLLEUR,bw=1400,from=0,to=5000)
DLLLEUR2 <- density(LLLEUR,bw=1400,from=15000,to=60000)

DLLLASN <- density(LLLASN,bw=1400,from=0,to=60000)
DLLLASN1 <- density(LLLASN,bw=1400,from=0,to=5000)
DLLLASN2 <- density(LLLASN,bw=1400,from=15000,to=60000)


DLLLAFRd <- density(LLLAFRd,bw=1400,from=0,to=60000)
DLLLAFR1d <- density(LLLAFRd,bw=1400,from=0,to=5000)
DLLLAFR2d <- density(LLLAFRd,bw=1400,from=15000,to=60000)

DLLLEURd <- density(LLLEURd,bw=1400,from=0,to=60000)
DLLLEUR1d <- density(LLLEURd,bw=1400,from=0,to=5000)
DLLLEUR2d <- density(LLLEURd,bw=1400,from=15000,to=60000)

DLLLASNd <- density(LLLASNd,bw=1400,from=0,to=60000)
DLLLASN1d <- density(LLLASNd,bw=1400,from=0,to=5000)
DLLLASN2d <- density(LLLASNd,bw=1400,from=15000,to=60000)

NEL <- NeanderMat[NE,34]


DNEL <- density(NEL,bw=1400,from=0,to=60000)
NEASN <- intersect(ASN,NE)
NEASNL <- NeanderMat[NEASN,34]
DNEASNL <- density(NEASNL,bw=1400,from=0,to=60000)
NEAFR <- intersect(AFR,NE)
NEAFRL <- NeanderMat[NEAFR,34]
DNEAFRL <- density(NEAFRL,bw=1400,from=0,to=60000)
NEEUR <- intersect(EUR,NE)
NEEURL <- NeanderMat[NEEUR,34]
DNEEURL <- density(NEEURL,bw=1400,from=0,to=60000)


NEASNd <- intersect(ASNd,NE) # 392
NEASNLd <- NeanderMat[NEASNd,34]
DNEASNLd <- density(NEASNLd,bw=1400,from=0,to=60000)
NEAFRd <- intersect(AFRd,NE)
NEAFRLd <- NeanderMat[NEAFRd,34]
DNEAFRLd <- density(NEAFRLd,bw=1400,from=0,to=60000)
NEEURd <- intersect(EURd,NE) # 206
NEEURLd <- NeanderMat[NEEURd,34]
DNEEURLd <- density(NEEURLd,bw=1400,from=0,to=60000)



DneL <- density(neL,bw=1400,from=0,to=60000)
DneL1 <- density(neL,bw=1400,from=0,to=5000)
DneL2 <- density(neL,bw=1400,from=15000,to=60000)
neASN <- intersect(ASN,ne)
neASNL <- NeanderMat[neASN,34]
DneASNL <- density(neASNL,bw=1400,from=0,to=60000)
neAFR <- intersect(AFR,ne)
neAFRL <- NeanderMat[neAFR,34]
DneAFRL <- density(neAFRL,bw=1400,from=0,to=60000)
neEUR <- intersect(EUR,ne)
neEURL <- NeanderMat[neEUR,34]
DneEURL <- density(neEURL,bw=1400,from=0,to=60000)


neASNd <- intersect(ASNd,ne) # 1303
neASNLd <- NeanderMat[neASNd,34]
DneASNLd <- density(neASNLd,bw=1400,from=0,to=60000)
neAFRd <- intersect(AFRd,ne)
neAFRLd <- NeanderMat[neAFRd,34]
DneAFRLd <- density(neAFRLd,bw=1400,from=0,to=60000)
neEURd <- intersect(EURd,ne) #886
neEURLd <- NeanderMat[neEURd,34]
DneEURLd <- density(neEURLd,bw=1400,from=0,to=60000)

################################################

DENEL <- DenisovaMat[DENE,34]


DDENEL <- density(DENEL,bw=1400,from=0,to=60000)
DENEASN <- intersect(ASN,DENE)
DENEASNL <- DenisovaMat[DENEASN,34]
DDENEASNL <- density(DENEASNL,bw=1400,from=0,to=60000)
DENEAFR <- intersect(AFR,DENE)
DENEAFRL <- DenisovaMat[DENEAFR,34]
DDENEAFRL <- density(DENEAFRL,bw=1400,from=0,to=60000)
DENEEUR <- intersect(EUR,DENE)
DENEEURL <- DenisovaMat[DENEEUR,34]
DDENEEURL <- density(DENEEURL,bw=1400,from=0,to=60000)


DENEASNd <- intersect(ASNd,DENE) # 392
DENEASNLd <- DenisovaMat[DENEASNd,34]
DDENEASNLd <- density(DENEASNLd,bw=1400,from=0,to=60000)
DENEAFRd <- intersect(AFRd,DENE)
DENEAFRLd <- DenisovaMat[DENEAFRd,34]
DDENEAFRLd <- density(DENEAFRLd,bw=1400,from=0,to=60000)
DENEEURd <- intersect(EURd,DENE) # 206
DENEEURLd <- DenisovaMat[DENEEURd,34]
DDENEEURLd <- density(DENEEURLd,bw=1400,from=0,to=60000)



DdeneL <- density(deneL,bw=1400,from=0,to=60000)
DdeneL1 <- density(deneL,bw=1400,from=0,to=5000)
DdeneL2 <- density(deneL,bw=1400,from=15000,to=60000)
deneASN <- intersect(ASN,dene)
deneASNL <- DenisovaMat[deneASN,34]
DdeneASNL <- density(deneASNL,bw=1400,from=0,to=60000)
deneAFR <- intersect(AFR,dene)
deneAFRL <- DenisovaMat[deneAFR,34]
DdeneAFRL <- density(deneAFRL,bw=1400,from=0,to=60000)
deneEUR <- intersect(EUR,dene)
deneEURL <- DenisovaMat[deneEUR,34]
DdeneEURL <- density(deneEURL,bw=1400,from=0,to=60000)


deneASNd <- intersect(ASNd,dene) # 1303
deneASNLd <- DenisovaMat[deneASNd,34]
DdeneASNLd <- density(deneASNLd,bw=1400,from=0,to=60000)
deneAFRd <- intersect(AFRd,dene)
deneAFRLd <- DenisovaMat[deneAFRd,34]
DdeneAFRLd <- density(deneAFRLd,bw=1400,from=0,to=60000)
deneEURd <- intersect(EURd,dene) #886
deneEURLd <- DenisovaMat[deneEURd,34]
DdeneEURLd <- density(deneEURLd,bw=1400,from=0,to=60000)

################################################
################################################
################################################


################################################
################################################
#DENISOVA LENGTH to all length
################################################

#Histogramm length
x11()
hist(LLL,breaks=500,ylim=c(0,5500),xlim=c(0,60000),plot=TRUE,col="blue",main="Histogram IBD Length",xlab="length in bases",ylab="counts")
abline(v=24000)

name="IBDlengthHistogram"
#pdf(paste(name,".pdf",sep=""),width=8)
#xfig(paste(name,".fig",sep=""),width=8)

df <- data.frame(counts=LLL)

ggplot(df, aes(x=counts)) +
  geom_histogram(aes(fill = ..count.. < 10000),colour="darkgreen",fill="coral",binwidth=1000) +
#  geom_histogram(aes(fill = cut(..count..,c(1000,2000,3000,4000,5000,6000))),binwidth=1000) +
    scale_y_continuous(breaks=c(0,1000,2000,3000,4000,5000,6000,7000,8000)) +
    scale_x_continuous(limits=c(0,60000)) +
    geom_vline(xintercept=c(24200),  linetype="dashed", size=0.5) +
    opts(title = expression("Histogram IBD Lengths")) +
    opts(legend.position="none")+
#    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(0,8300),xlim=c(0,60000)) +
    ylab("Counts") +
    xlab("IBD Length")


dev.off()


#histogramm length denisova
x11()
hist(deL,breaks=500,ylim=c(0,300),xlim=c(0,60000),plot=TRUE,col="green")
abline(v=3100)
abline(v=24000)

name="IBDlengthHistogramDenisova"
#pdf(paste(name,".pdf",sep=""),width=8)
#xfig(paste(name,".fig",sep=""),width=8)


df <- data.frame(counts=deL)

ggplot(df, aes(x=counts)) +
  geom_histogram(aes(fill = ..count.. < 310),colour="darkgreen",fill="coral",binwidth=1000) +
#  geom_histogram(aes(fill = cut(..count..,c(1000,2000,3000,4000,5000,6000))),binwidth=1000) +
    scale_y_continuous(breaks=c(0,50,100,150,200)) +
    scale_x_continuous(limits=c(0,60000)) +
    geom_vline(xintercept=c(9000,24200),  linetype="dashed", size=0.5) +
    opts(title = expression("Histogram Denisovan IBD Lengths")) +
    opts(legend.position="none")+
#    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(0,220),xlim=c(0,60000)) +
    ylab("Counts") +
    xlab("IBD Length")

dev.off()


#Length densities
x11()
plot(DLLL$x,DLLL$y,type="l",ylim=c(0,0.0003))
abline(h=0)
points(DLLL$x,DdeL$y,type="l",col="red")
abline(v=500)
abline(v=24000)


name="IBDlengthDensityALLDenisova"
#pdf(paste(name,".pdf",sep=""),width=8)
#xfig(paste(name,".fig",sep=""),width=8)



l1 <- length(LLL)
l2 <- length(deL)

df <- data.frame(IBD_Length = factor( c(rep("zAll",l1),rep("Denisovan",l2)) ),
                   proportion = c(LLL,deL))

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

ggplot(df, aes(x=proportion, fill=IBD_Length)) + geom_density(alpha=.3,adjust=0.7) +
    geom_vline(data=cdf, aes(xintercept=c(9000,24200),  colour=IBD_Length),
               linetype="dashed", size=1) +
   opts(title = expression("IBD Length Distributions: All vs. Denisova")) +
   scale_y_continuous(breaks=c(0,1e-05,2e-05,3e-05,4e-05,5e-05)) +
#    scale_y_continuous(limits=c(0,4e-04)) +
    scale_x_continuous(limits=c(0,60000)) +
    opts(legend.position=c(.8, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(0,5.2e-05),xlim=c(0,60000)) +
    ylab("Density") +
    xlab("IBD Length")

dev.off()



#difference of length densities
x11()
plot(DLLL$x,(DdeL$y-DLLL$y),type="l",ylim=c(-1e-05,1e-05))
abline(h=0)
abline(v=6400)
abline(v=28000)
abline(v=23000)
abline(v=19000)



name="IBDlengthDensityDiffDenisova"
#pdf(paste(name,".pdf",sep=""),width=8)
#xfig(paste(name,".fig",sep=""),width=8)


                                        # The palette with grey:
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# The palette with black:
#cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")


y0=DdeL$y-DLLL$y

df <- data.frame(IBD_Length=DLLL$x,Density_Difference=y0)

ggplot(df,aes(x=IBD_Length,y=Density_Difference,fill="",colour="")) + geom_area(size=1) +
  scale_fill_manual(values=cbPalette[8]) +
  scale_colour_manual(values=cbPalette[4]) +
    geom_vline(data=df, aes(xintercept=c(6400,24500,39500)),linetype="dashed", size=0.5) +
    geom_hline(data=df, aes(yintercept=c(0),colour="green"),linetype="solid", size=1) +
    opts(title = expression("Difference IBD Length Distributions: Denisova - All")) +
    scale_y_continuous(breaks=c(-1e-05,-0.5e-05,0,0.5e-05,0.6e-05)) +
    scale_x_continuous(limits=c(0,60000)) +
     opts(legend.position="none")+
#    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(-1.2e-05,6e-06),xlim=c(0,60000)) +
    ylab("Density Difference") +
    xlab("IBD Length")

dev.off()

################################################
################################################
#DENISOVA LENGTH per Population
################################################

#Overall
plot(DLLL$x,DdeL$y,type="l",ylim=c(0,0.00003))
abline(h=0)

#ASN three peaks
points(DLLLASN$x,DdeASNL$y,type="l",col="red")

#AFR like overall
points(DLLLAFR$x,DdeAFRL$y,type="l",col="blue")

#EUR: two peaks which allmost coincide with ASN
points(DLLLEUR$x,DdeEURL$y,type="l",col="green")
abline(v=28000)
abline(v=23000)
abline(v=19000)


name="IBDlengthDensityDenisovaPop"
#pdf(paste(name,".pdf",sep=""),width=8)
#xfig(paste(name,".fig",sep=""),width=8)




df <- data.frame(Length=DdeL$x,Denisovan1_All=DdeL$y,Denisovan2_African=DdeAFRL$y,Denisovan3_Asian=DdeASNL$y,Denisovan4_European=DdeEURL$y)

ggplot(df, aes(x=Length)) +
    geom_line(aes(y=Denisovan1_All,colour="Denisovan All"),size=1) +
    geom_line(aes(y=Denisovan2_African,colour="Denisovan African"),size=1) +
    geom_line(aes(y=Denisovan3_Asian,colour="Denisovan Asian"),size=1) +
    geom_line(aes(y=Denisovan4_European,colour="Denisovan European"),size=1) +
    geom_hline(data=df, yintercept=c(0),linetype="solid", size=0.5) +
    geom_vline(data=df,  aes(xintercept=c(8500,18000,24400,25500,40000),colour=c("Denisovan African","Denisovan European","Denisovan Asian","Denisovan European","Denisovan European")),linetype="dashed", size=0.6) +
    scale_y_continuous(limits=c(0,0.000053)) +
    opts(title = expression("Denisovan Matching IBD Length Distributions per Population")) +
    opts(legend.position=c(.8, .6))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim =c(0,0.000053),xlim=c(0,60000)) +
    ylab("Density") +
    xlab("IBD Length")


dev.off()

################################################
################################################
################################################

################################################
################################################
################################################

################################################
################################################
################################################


################################################
################################################
################################################


################################################
################################################
#Neander LENGTH to all length
################################################

#Histogramm length
x11()
hist(LLL,breaks=500,ylim=c(0,5500),xlim=c(0,60000),plot=TRUE,col="blue")
abline(v=24000)

#histogramm length neander
x11()
hist(neL,breaks=500,ylim=c(0,350),xlim=c(0,60000),plot=TRUE,col="green")
abline(v=3100)
abline(v=24000)

name="IBDlengthHistogramNeander"
#pdf(paste(name,".pdf",sep=""),width=8)
#xfig(paste(name,".fig",sep=""),width=8)

df <- data.frame(counts=neL)

ggplot(df, aes(x=counts)) +
  geom_histogram(aes(fill = ..count.. < 360),colour="darkgreen",fill="coral",binwidth=1000) +
#  geom_histogram(aes(fill = cut(..count..,c(1000,2000,3000,4000,5000,6000))),binwidth=1000) +
    scale_y_continuous(breaks=c(0,50,100,150,200,250,300,350)) +
    scale_x_continuous(limits=c(0,60000)) +
    geom_vline(xintercept=c(18500,43000),  linetype="dashed", size=0.5) +
    opts(title = expression("Histogram Neandertal IBD Lengths")) +
    opts(legend.position="none")+
#    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(0,350),xlim=c(0,60000)) +
    ylab("Counts") +
    xlab("IBD Length")

dev.off()

#Length densities
x11()
plot(DLLL$x,DLLL$y,type="l",ylim=c(0,0.0003))
abline(h=0)
points(DLLL$x,DneL$y,type="l",col="red")
abline(v=500)
abline(v=24000)

name="IBDlengthDensityALLNeander"
#pdf(paste(name,".pdf",sep=""),width=8)
xfig(paste(name,".fig",sep=""),width=8)


l1 <- length(LLL)
l2 <- length(neL)

df <- data.frame(IBD_Length = factor( c(rep("zAll",l1),rep("Neandertal",l2)) ),
                   proportion = c(LLL,neL))

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

ggplot(df, aes(x=proportion, fill=IBD_Length)) + geom_density(alpha=.3,adjust=0.6) +
    geom_vline(data=cdf, aes(xintercept=c(18000,24200),color=IBD_Length),
               linetype="dashed", size=1) +
    geom_vline(data=cdf, aes(xintercept=c(43000,24200),color=IBD_Length),
               linetype="dashed", size=1) +
   opts(title = expression("IBD Length Distributions: All vs. Neandertal")) +
    scale_y_continuous(breaks=c(0,1e-05,2e-05,3e-05,4e-05,5e-05)) +
#    scale_y_continuous(limits=c(0,4e-04)) +
    scale_x_continuous(limits=c(0,60000)) +
    opts(legend.position=c(.8, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(0,5.2e-05),xlim=c(0,60000)) +
    ylab("Density") +
    xlab("IBD Length")

dev.off()


#difference of length densities
x11()
plot(DLLL$x,(DneL$y-DLLL$y),type="l",ylim=c(-1e-05,1e-05))
abline(h=0)
abline(v=6800)
abline(v=26700)
abline(v=21500)

name="IBDlengthDensityDiffNeander"
#pdf(paste(name,".pdf",sep=""),width=8)
#xfig(paste(name,".fig",sep=""),width=8)


                                        # The palette with grey:
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# The palette with black:
#cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")


y0=DneL$y-DLLL$y

df <- data.frame(IBD_Length=DLLL$x,Density_Difference=y0)

ggplot(df,aes(x=IBD_Length,y=Density_Difference,fill="",colour="")) + geom_area(size=1) +
  scale_fill_manual(values=cbPalette[8]) +
  scale_colour_manual(values=cbPalette[4]) +
    geom_vline(data=df, aes(xintercept=c(11000,24200,43000)),linetype="dashed", size=0.5) +
    geom_hline(data=df, aes(yintercept=c(0),colour="green"),linetype="solid", size=1) +
    opts(title = expression("Difference IBD Length Distributions: Neandertal - All")) +
    scale_y_continuous(breaks=c(-0.8e-05,-0.5e-05,0,0.4e-05)) +
    scale_x_continuous(limits=c(0,60000)) +
     opts(legend.position="none")+
#    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(-0.8e-05,0.4e-05),xlim=c(0,60000)) +
    ylab("Density Difference") +
    xlab("IBD Length")

dev.off()



################################################
################################################
#Neander LENGTH per Population
################################################

#Overall
plot(DLLL$x,DneL$y,type="l",ylim=c(0,0.00003))
abline(h=0)

#ASN two peaks
points(DLLLASN$x,DneASNL$y,type="l",col="red")

#AFR like overall
points(DLLLAFR$x,DneAFRL$y,type="l",col="blue")

#EUR: two peaks which allmost coincide with ASN
points(DLLLEUR$x,DneEURL$y,type="l",col="green")
abline(v=26700)
abline(v=21500)


name="IBDlengthDensityNeanderPop"
#pdf(paste(name,".pdf",sep=""),width=8)
#xfig(paste(name,".fig",sep=""),width=8)


df <- data.frame(Length=DneL$x,Neandertal1_All=DneL$y,Neandertal2_African=DneAFRL$y,Neandertal3_Asian=DneASNL$y,Neandertal4_European=DneEURL$y)

ggplot(df, aes(x=Length)) +
    geom_line(aes(y=Neandertal1_All,colour="Neandertal All"),size=1) +
    geom_line(aes(y=Neandertal2_African,colour="Neandertal African"),size=1) +
    geom_line(aes(y=Neandertal3_Asian,colour="Neandertal Asian"),size=1) +
    geom_line(aes(y=Neandertal4_European,colour="Neandertal European"),size=1) +
    geom_hline(data=df, yintercept=c(0),linetype="solid", size=0.5) +
    geom_vline(data=df, aes(xintercept=c(17000,24000,25800,42000),colour=c("Neandertal African","Neandertal European","Neandertal Asian","Neandertal Asian")), linetype="dashed", size=0.6) +
    scale_y_continuous(limits=c(0,0.00005)) +
    opts(title = expression("Neandertal Matching IBD Length Distributions per Population")) +
    opts(legend.position=c(.8, .6))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim =c(0,0.00005),xlim=c(0,60000)) +
    ylab("Density") +
    xlab("IBD Length")

dev.off()

################################################
################################################
################################################

################################################
################################################
################################################

################################################
################################################
################################################


################################################
################################################
#Archaic LENGTH to all length
################################################

#Histogramm length
x11()
hist(LLL,breaks=500,ylim=c(0,5500),xlim=c(0,60000),plot=TRUE,col="blue")
abline(v=24000)

#histogramm length archaic
x11()
hist(neL,breaks=500,ylim=c(0,350),xlim=c(0,60000),plot=TRUE,col="green")
abline(v=3100)
abline(v=24000)

name="IBDlengthHistogramArchaic"
#pdf(paste(name,".pdf",sep=""),width=8)
#xfig(paste(name,".fig",sep=""),width=8)

df <- data.frame(counts=deneL)

ggplot(df, aes(x=counts)) +
  geom_histogram(aes(fill = ..count.. < 360),colour="darkgreen",fill="coral",binwidth=1000) +
#  geom_histogram(aes(fill = cut(..count..,c(1000,2000,3000,4000,5000,6000))),binwidth=1000) +
    scale_y_continuous(breaks=c(0,50,100,150,200)) +
    scale_x_continuous(limits=c(0,60000)) +
    geom_vline(xintercept=c(11200,24200,41000),  linetype="dashed", size=0.5) +
    opts(title = expression("Histogram Archaic IBD Lengths")) +
    opts(legend.position="none")+
#    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(0,230),xlim=c(0,60000)) +
    ylab("Counts") +
    xlab("IBD Length")

dev.off()

#Length densities
x11()
plot(DLLL$x,DLLL$y,type="l",ylim=c(0,0.0003))
abline(h=0)
points(DLLL$x,DdeneL$y,type="l",col="red")
abline(v=500)
abline(v=24000)

name="IBDlengthDensityALLArchaic"
#pdf(paste(name,".pdf",sep=""),width=8)
xfig(paste(name,".fig",sep=""),width=8)


l1 <- length(LLL)
l2 <- length(deneL)

df <- data.frame(IBD_Length = factor( c(rep("zAll",l1),rep("Archaic",l2)) ),
                   proportion = c(LLL,deneL))

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

ggplot(df, aes(x=proportion, fill=IBD_Length)) + geom_density(alpha=.3,adjust=0.6) +
    geom_vline(data=cdf, aes(xintercept=c(11200,24200),color=IBD_Length),
               linetype="dashed", size=1) +
    geom_vline(data=cdf, aes(xintercept=c(41000,24200),color=IBD_Length),
               linetype="dashed", size=1) +
   opts(title = expression("IBD Length Distributions: All vs. Archaic")) +
    scale_y_continuous(breaks=c(0,1e-05,2e-05,3e-05,4e-05,5e-05)) +
#    scale_y_continuous(limits=c(0,4e-04)) +
    scale_x_continuous(limits=c(0,60000)) +
    opts(legend.position=c(.8, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(0,5.2e-05),xlim=c(0,60000)) +
    ylab("Density") +
    xlab("IBD Length")

dev.off()


#difference of length densities
x11()
plot(DLLL$x,(DdeneL$y-DLLL$y),type="l",ylim=c(-1e-05,1e-05))
abline(h=0)
abline(v=6800)
abline(v=26700)
abline(v=21500)

name="IBDlengthDensityDiffArchaic"
#pdf(paste(name,".pdf",sep=""),width=8)
#xfig(paste(name,".fig",sep=""),width=8)


                                        # The palette with grey:
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# The palette with black:
#cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")


y0=DdeneL$y-DLLL$y

df <- data.frame(IBD_Length=DLLL$x,Density_Difference=y0)

ggplot(df,aes(x=IBD_Length,y=Density_Difference,fill="",colour="")) + geom_area(size=1) +
  scale_fill_manual(values=cbPalette[8]) +
  scale_colour_manual(values=cbPalette[4]) +
    geom_vline(data=df, aes(xintercept=c(11200,24200,42000)),linetype="dashed", size=0.5) +
    geom_hline(data=df, aes(yintercept=c(0),colour="green"),linetype="solid", size=1) +
    opts(title = expression("Difference IBD Length Distributions: Archaic - All")) +
    scale_y_continuous(breaks=c(-0.8e-05,-0.5e-05,0,0.4e-05)) +
    scale_x_continuous(limits=c(0,60000)) +
     opts(legend.position="none")+
#    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim = c(-0.8e-05,0.4e-05),xlim=c(0,60000)) +
    ylab("Density Difference") +
    xlab("IBD Length")

dev.off()



################################################
################################################
#Archaic LENGTH per Population
################################################

#Overall
plot(DLLL$x,DdeneL$y,type="l",ylim=c(0,0.00003))
abline(h=0)

#ASN two peaks
points(DLLLASN$x,DdeneASNL$y,type="l",col="red")

#AFR like overall
points(DLLLAFR$x,DdeneAFRL$y,type="l",col="blue")

#EUR: two peaks which allmost coincide with ASN
points(DLLLEUR$x,DdeneEURL$y,type="l",col="green")
abline(v=26700)
abline(v=21500)


name="IBDlengthDensityArchaicPop"
#pdf(paste(name,".pdf",sep=""),width=8)
#xfig(paste(name,".fig",sep=""),width=8)


df <- data.frame(Length=DdeneL$x,Archaic1_All=DdeneL$y,Archaic2_African=DdeneAFRL$y,Archaic3_Asian=DdeneASNL$y,Archaic4_European=DdeneEURL$y)

ggplot(df, aes(x=Length)) +
    geom_line(aes(y=Archaic1_All,colour="Archaic All"),size=1) +
    geom_line(aes(y=Archaic2_African,colour="Archaic African"),size=1) +
    geom_line(aes(y=Archaic3_Asian,colour="Archaic Asian"),size=1) +
    geom_line(aes(y=Archaic4_European,colour="Archaic European"),size=1) +
    geom_hline(data=df, yintercept=c(0),linetype="solid", size=0.5) +
    geom_vline(data=df, aes(xintercept=c(12000,22000,24400,42000),colour=c("Archaic European","Archaic African","Archaic Asian","Archaic Asian")), linetype="dashed", size=0.6) +
    scale_y_continuous(limits=c(0,0.000038)) +
    opts(title = expression("Archaic Matching IBD Length Distributions per Population")) +
    opts(legend.position=c(.8, .6))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+
    coord_cartesian(ylim =c(0,0.000038),xlim=c(0,60000)) +
    ylab("Density") +
    xlab("IBD Length")

dev.off()







################################################
################################################
################################################

################################################
################################################
################################################

################################################
################################################
################################################

# Where do  Denisova  contribute to
# length ofAsians, Africans, Europeans

plot(DLLL$x,(DdeL$y-DLLL$y),type="l",ylim=c(-2.8e-05,2.5e-05))
abline(h=0)

points(DLLLASN$x,(DdeASNL$y-DLLLASN$y),type="l",col="red")
abline(h=0)

points(DLLLAFR$x,(DdeAFRL$y-DLLLAFR$y),type="l",col="blue")
abline(h=0)

points(DLLLEUR$x,(DdeEURL$y-DLLLEUR$y),type="l",col="green")
abline(h=0)



# Where do Asians, Africans, Europeans contribute to
# length of Denisova
plot(DLLL$x,(DdeL$y-DLLL$y),type="l")
abline(h=0)


points(DLLLASN$x,(DdeASNL$y-DdeL$y),type="l",col="red")
abline(h=0)

points(DLLLAFR$x,(DdeAFRL$y-DdeL$y),type="l",col="blue")
abline(h=0)

points(DLLLEUR$x,(DdeEURL$y-DdeL$y),type="l",col="green")
abline(h=0)





plot(DLLL$x,DdeL$y,type="l",ylim=c(0,0.00004))
abline(h=0)


points(DLLLASN$x,DdeASNL$y,type="l",col="red")
abline(h=0)

points(DLLLAFR$x,DdeAFRL$y,type="l",col="blue")
abline(h=0)

points(DLLLEUR$x,DdeEURL$y,type="l",col="green")
abline(h=0)



plot(DLLL$x,DLLL$y,type="l")
abline(h=0)


points(DLLLASN$x,DLLLASN$y,type="l",col="red")

points(DLLLAFR$x,DLLLAFR$y,type="l",col="blue")

points(DLLLEUR$x,DLLLEUR$y,type="l",col="green")



plot(DLLL$x,DLLL$y,type="l",ylim=c(-2.8e-05,0.00012))
abline(h=0)


points(DLLLASN$x,(DLLLASN$y-DLLL$y),type="l",col="red")

points(DLLLAFR$x,(DLLLAFR$y-DLLL$y),type="l",col="blue")

points(DLLLEUR$x,(DLLLEUR$y-DLLL$y),type="l",col="green")

abline(v=500)


x11()
hist(LLLASN,breaks=500,ylim=c(0,200),xlim=c(0,60000),plot=TRUE,col="red")

x11()
hist(LLLAFR,breaks=500,ylim=c(0,5000),xlim=c(0,60000),plot=TRUE,col="blue")

x11()
hist(LLLEUR,breaks=500,ylim=c(0,300),xlim=c(0,60000),plot=TRUE,col="green")


x11()
hist(deASNL,breaks=200,ylim=c(0,40),xlim=c(0,60000),plot=TRUE,col="red")

x11()
hist(deEURL,breaks=200,ylim=c(0,30),xlim=c(0,60000),plot=TRUE,col="green")

x11()
hist(deAFRL,breaks=500,ylim=c(0,300),xlim=c(0,60000),plot=TRUE,col="blue")

###################################################
###################################################


plot(DLLL$x,DLLL$y,type="l")
points(DdeL$x,DdeL$y,type="l",col="red")

plot(DLLL1$x,DLLL1$y,type="l")
points(DdeL1$x,DdeL1$y,type="l",col="red")

plot(DLLL2$x,DLLL2$y,type="l")
points(DdeL2$x,DdeL2$y,type="l",col="red")

################################################

> plot(DLLL$x,(DdeL$y-DLLL$y),type="l")
> abline(h=0)
> points(DLLL$x,(DDEL$y-DLLL$y),type="l",col="red")
################################################


plot(DLLL$x,(DDEL$y-DLLL$y),type="l")
abline(h=0)


points(DLLLASNd$x,(DDEASNLd$y-DLLLASNd$y),type="l",col="red")
abline(h=0)

points(DLLLAFRd$x,(DDEAFRLd$y-DLLLAFRd$y),type="l",col="blue")
abline(h=0)

points(DLLLEURd$x,(DDEEURLd$y-DLLLEURd$y),type="l",col="green")
abline(h=0)


hist(DEEURLd,breaks=100,ylim=c(0,30),xlim=c(0,60000),plot=TRUE,col="blue")

hist(DEEURL[which(DEEURL>1000)],breaks=100,xlim=c(0,60000),plot=TRUE,col="blue")



plot(DLLL$x,DDEL$y,type="l")
abline(h=0)


points(DLLLASNd$x,DLLLASNd$y,type="l",col="red")
abline(h=0)

points(DLLLAFRd$x,DLLLAFRd$y,type="l",col="blue")
abline(h=0)

points(DLLLEURd$x,DLLLEURd$y,type="l",col="green")
abline(h=0)

################################################
################################################


plot(DLLL$x,(DdeL$y-DLLL$y),type="l")
abline(h=0)

plot(DLLL1$x,(DdeL1$y-DLLL1$y),type="l")
abline(h=0)

plot(DLLL2$x,(DdeL2$y-DLLL2$y),type="l")
abline(h=0)



points(DLLLASNd$x,(DdeASNLd$y-DLLLASNd$y),type="l",col="red")
abline(h=0)

points(DLLLAFRd$x,(DdeAFRLd$y-DLLLAFRd$y),type="l",col="blue")
abline(h=0)

points(DLLLEURd$x,(DdeEURLd$y-DLLLEURd$y),type="l",col="green")
abline(h=0)



plot(DLLL$x,DdeL$y,type="l")
abline(h=0)


points(DLLLASNd$x,DdeASNLd$y,type="l",col="red")
abline(h=0)

points(DLLLAFRd$x,DdeAFRLd$y,type="l",col="blue")
abline(h=0)

points(DLLLEURd$x,DdeEURLd$y,type="l",col="green")
abline(h=0)


hist(deASNLd,breaks=100,ylim=c(0,60),xlim=c(0,60000),plot=TRUE,col="blue")

hist(deASNL[which(deASNL>1000)],breaks=100,xlim=c(0,60000),plot=TRUE,col="blue")



plot(DLLL$x,DdeL$y,type="l")
abline(h=0)


points(DLLLASNd$x,DLLLASNd$y,type="l",col="red")
abline(h=0)

points(DLLLAFRd$x,DLLLAFRd$y,type="l",col="blue")
abline(h=0)

points(DLLLEURd$x,DLLLEURd$y,type="l",col="green")
abline(h=0)

################################################
plot(DLLL$x,(DDEL$y-DLLL$y),type="l")
abline(h=0)


points(DLLLASN$x,(DDEASNL$y-DLLLASN$y),type="l",col="red")
abline(h=0)

points(DLLLAFR$x,(DDEAFRL$y-DLLLAFR$y),type="l",col="blue")
abline(h=0)

points(DLLLEUR$x,(DDEEURL$y-DLLLEUR$y),type="l",col="green")
abline(h=0)



plot(DLLL$x,DDEL$y,type="l")
abline(h=0)


points(DLLLASN$x,DLLLASN$y,type="l",col="red")
abline(h=0)

points(DLLLAFR$x,DLLLAFR$y,type="l",col="blue")
abline(h=0)

points(DLLLEUR$x,DLLLEUR$y,type="l",col="green")
abline(h=0)


####################################################

plot(DLLL1$x,(DdeL1$y-DLLL1$y),type="l")
abline(h=0)

plot(DLLL$x,(DdeL$y-DLLL$y),type="l")
abline(h=0)


points(DLLLASN$x,(DdeASNL$y-DLLLASN$y),type="l",col="red")
abline(h=0)

points(DLLLAFR$x,(DdeAFRL$y-DLLLAFR$y),type="l",col="blue")
abline(h=0)

points(DLLLEUR$x,(DdeEURL$y-DLLLEUR$y),type="l",col="green")
abline(h=0)




plot(DLLL$x,(DdeL$y-DLLL$y),type="l")
abline(h=0)


points(DLLLASN$x,(DdeASNL$y-DdeL$y),type="l",col="red")
abline(h=0)

points(DLLLAFR$x,(DdeAFRL$y-DdeL$y),type="l",col="blue")
abline(h=0)

points(DLLLEUR$x,(DdeEURL$y-DdeL$y),type="l",col="green")
abline(h=0)







plot(DLLL$x,DdeL$y,type="l")
abline(h=0)


points(DLLLASN$x,DdeASNL$y,type="l",col="red")
abline(h=0)

points(DLLLAFR$x,DdeAFRL$y,type="l",col="blue")
abline(h=0)

points(DLLLEUR$x,DdeEURL$y,type="l",col="green")
abline(h=0)



plot(DLLL$x,DLLL$y,type="l")
abline(h=0)


points(DLLLASN$x,DLLLASN$y,type="l",col="red")
abline(h=0)

points(DLLLAFR$x,DLLLAFR$y,type="l",col="blue")
abline(h=0)

points(DLLLEUR$x,DLLLEUR$y,type="l",col="green")
abline(h=0)

##################################################




DEASN <- intersect(ASN,DE)

DEASNL <- DenisovaMat[DEASN,34]

plot(density(DEASNL,bw=1400,from=0,to=60000))


DDEASNL <- density(DEASNL,bw=1400,from=0,to=60000)

plot(DLLL$x,(DDEASNL$y-DLLL$y),type="l")
abline(h=0)




DEAFR <- intersect(AFR,DE)

DEAFRL <- DenisovaMat[DEAFR,34]

plot(density(DEAFRL,bw=1400,from=0,to=60000))

DDEAFRL <- density(DEAFRL,bw=1400,from=0,to=60000)

plot(DLLL$x,(DDEAFRL$y-DLLL$y),type="l")
abline(h=0)



DEEUR <- intersect(EUR,DE)

DEEURL <- DenisovaMat[DEEUR,34]

plot(density(DEEURL,bw=1400,from=0,to=60000))

DDEEURL <- density(DEEURL,bw=1400,from=0,to=60000)

plot(DLLL$x,(DDEEURL$y-DLLL$y),type="l")
abline(h=0)


DDEEURL1 <- density(DEEURL,bw=1400,from=0,to=5000)
plot(DLLL1$x,(DDEEURL1$y-DLLL1$y),type="l")
abline(h=0)


DDEEURL2 <- density(DEEURL,bw=1400,from=15000,to=60000)
plot(DLLL2$x,(DDEEURL2$y-DLLL2$y),type="l")
abline(h=0)


neL <- NeanderMat[ne,34]

plot(density(neL))

hist(neL,breaks=10000,xlim=c(0,60000),plot=TRUE)
hist(neL,breaks=10000,xlim=c(0,1000),plot=TRUE)

plot(density(neL,bw=1400,from=20000,to=30000))


neL <- NeanderMat[ne,34]

deL <- DenisovaMat[de,34]

plot(density(deL))

plot(density(deL,bw=1400,from=0,to=60000))
plot(density(deL,bw=1400,from=20000,to=30000))

> hist(deL,breaks=10000,xlim=c(0,60000),plot=TRUE)
> hist(deL,breaks=10000,xlim=c(0,1000),plot=TRUE)

neL <- NeanderMat[ne,34]

plot(density(neL))

hist(neL,breaks=10000,xlim=c(0,60000),plot=TRUE)
hist(neL,breaks=10000,xlim=c(0,1000),plot=TRUE)

plot(density(neL,bw=1400,from=20000,to=30000))





deASN <- intersect(ASN,de)

deASNL <- DenisovaMat[deASN,34]

plot(density(deASNL))

neASN <- intersect(ASN,ne)

neASNL <-  NeanderMat[neASN,34]

plot(density(neASNL,bw=1400,from=0,to=60000))



deAFR <- intersect(AFR,de)

deAFRL <- DenisovaMat[deAFR,34]

plot(density(deAFRL,bw=1400,from=0,to=60000))

neAFR <- intersect(AFR,ne)

neAFRL <-  NeanderMat[neAFR,34]

plot(density(neAFRL,bw=1400,from=0,to=60000))



deEUR <- intersect(EUR,de)

deEURL <- DenisovaMat[deEUR,34]

plot(density(deEURL,bw=1400,from=0,to=60000))

neEUR <- intersect(EUR,ne)

neEURL <-  NeanderMat[neEUR,34]

plot(density(neEURL,bw=1400,from=0,to=60000))




deANSAFR <-  intersect(deASN,AFR1)




aeASN <- intersect(ASN1,ae)
aeANSAFR <-  intersect(aeASN,AFR1)

au <- which(AncestMat[,5]>30)

auANSAFR <- intersect(au,aeANSAFR)



plot(density(DEL,bw=1400,from=0,to=60000))
points(density(DEAFRL,bw=1400,from=0,to=60000),type="l",col="blue")
points(density(DEASNL,bw=1400,from=0,to=60000),type="l",col="red")
points(density(DEEURL,bw=1400,from=0,to=60000),type="l",col="green")

plot(density(DEL,bw=1400,from=0,to=5000))
points(density(DEAFRL,bw=1400,from=0,to=5000),type="l",col="blue")
points(density(DEASNL,bw=1400,from=0,to=5000),type="l",col="red")
points(density(DEEURL,bw=1400,from=0,to=5000),type="l",col="green")

plot(density(DEL,bw=1400,from=15000,to=35000))
points(density(DEAFRL,bw=1400,from=15000,to=35000),type="l",col="blue")
points(density(DEASNL,bw=1400,from=15000,to=35000),type="l",col="red")
points(density(DEEURL,bw=1400,from=15000,to=35000),type="l",col="green")



plot(density(deL,bw=1400,from=0,to=60000))
points(density(deAFRL,bw=1400,from=0,to=60000),type="l",col="blue")
points(density(deASNL,bw=1400,from=0,to=60000),type="l",col="red")
points(density(deEURL,bw=1400,from=0,to=60000),type="l",col="green")

plot(density(deL,bw=1400,from=0,to=5000))
points(density(deAFRL,bw=1400,from=0,to=5000),type="l",col="blue")
points(density(deASNL,bw=1400,from=0,to=5000),type="l",col="red")
points(density(deEURL,bw=1400,from=0,to=5000),type="l",col="green")

plot(density(deL,bw=1400,from=15000,to=35000))
points(density(deAFRL,bw=1400,from=15000,to=35000),type="l",col="blue")
points(density(deASNL,bw=1400,from=15000,to=35000),type="l",col="red")
points(density(deEURL,bw=1400,from=15000,to=35000),type="l",col="green")




plot(density(DEL,bw=1400,from=0,to=60000))
points(density(DEAFRL,bw=1400,from=0,to=60000),type="l",col="blue")
points(density(DEASNL,bw=1400,from=0,to=60000),type="l",col="red")
points(density(DEEURL,bw=1400,from=0,to=60000),type="l",col="green")

hist(DEL,breaks=100,xlim=c(0,60000),plot=TRUE)
hist(DEAFRL,breaks=100,xlim=c(0,60000),plot=TRUE,col="blue")
hist(DEASNL,breaks=100,xlim=c(0,60000),plot=TRUE,col="red")
hist(DEEURL,breaks=100,xlim=c(0,60000),plot=TRUE,col="green")



hist(deL,breaks=1000,xlim=c(0,60000),plot=TRUE)
hist(deAFRL,breaks=1000,xlim=c(0,60000),plot=TRUE,col="blue")
hist(deASNL,breaks=400,xlim=c(0,60000),plot=TRUE,col="red")
hist(deEURL,breaks=400,xlim=c(0,60000),plot=TRUE,col="green")




#########################################################################
#########################################################################
#########################################################################

summary(DenisovaMat[,22])
summary(DenisovaMat[DE,22])
summary(DenisovaMat[NE,22])
summary(DenisovaMat[de,22])
summary(DenisovaMat[ne,22])

#########################################################################
#########################################################################
#########################################################################

#########################################################################
#########################################################################
#########################################################################


oAFR <- which(PopulationMat[,6]>0.9999)

oAMR <- which(PopulationMat[,7]>0.9999)

oASN <- which(PopulationMat[,8]>0.9999)

oEUR <- which(PopulationMat[,9]>0.9999)


NAFR <- which(PopulationMat[,6]<0.0001)

NAMR <- which(PopulationMat[,7]<0.0001)

NASN <- which(PopulationMat[,8]<0.0001)

NEUR <- which(PopulationMat[,9]<0.0001)

JAFR <- which(PopulationMat[,6]>0.0001)

JAMR <- which(PopulationMat[,7]>0.0001)

JASN <- which(PopulationMat[,8]>0.0001)

JEUR <- which(PopulationMat[,9]>0.0001)



deoASN <- intersect(oASN,de)
deoASNL <- DenisovaMat[deoASN,34]
DdeoASNL <- density(deoASNL,bw=3000,from=0,to=60000)
deoAFR <- intersect(oAFR,de)
deoAFRL <- DenisovaMat[deoAFR,34]
DdeoAFRL <- density(deoAFRL,bw=3000,from=0,to=60000)
deoEUR <- intersect(oEUR,de)
deoEURL <- DenisovaMat[deoEUR,34]
DdeoEURL <- density(deoEURL,bw=3000,from=0,to=60000)

neoASN <- intersect(oASN,ne)
neoASNL <- NeanderMat[neoASN,34]
DneoASNL <- density(neoASNL,bw=3000,from=0,to=60000)
neoAFR <- intersect(oAFR,ne)
neoAFRL <- NeanderMat[neoAFR,34]
DneoAFRL <- density(neoAFRL,bw=3000,from=0,to=60000)
neoEUR <- intersect(oEUR,ne)
neoEURL <- NeanderMat[neoEUR,34]
DneoEURL <- density(neoEURL,bw=3000,from=0,to=60000)

deneoASN <- intersect(oASN,dene)
deneoASNL <- DenisovaMat[deneoASN,34]
DdeneoASNL <- density(deneoASNL,bw=4000,from=0,to=60000)
deneoAFR <- intersect(oAFR,dene)
deneoAFRL <- DenisovaMat[deneoAFR,34]
DdeneoAFRL <- density(deneoAFRL,bw=4000,from=0,to=60000)
deneoEUR <- intersect(oEUR,dene)
deneoEURL <- DenisovaMat[deneoEUR,34]
DdeneoEURL <- density(deneoEURL,bw=4000,from=0,to=60000)



deNAFR <- intersect(NAFR,de)
deNAFRL <- DenisovaMat[deNAFR,34]
DdeNAFRL <- density(deNAFRL,bw=1400,from=0,to=60000)

neNAFR <- intersect(NAFR,ne)
neNAFRL <- NeanderMat[neNAFR,34]
DneNAFRL <- density(neNAFRL,bw=1400,from=0,to=60000)

deneNAFR <- intersect(NAFR,dene)
deneNAFRL <- DenisovaMat[deneNAFR,34]
DdeneNAFRL <- density(deneNAFRL,bw=1400,from=0,to=60000)


DdeL1 <- density(deL,bw=3000,from=0,to=60000)
DdeL2 <- density(deL,bw=4000,from=0,to=60000)

DneL1 <- density(neL,bw=3000,from=0,to=60000)
DneL2 <- density(neL,bw=4000,from=0,to=60000)

DdeneL1 <- density(deneL,bw=3000,from=0,to=60000)
DdeneL2 <- density(deneL,bw=4000,from=0,to=60000)

name="IBDlengthDensityDenisovaPopOnly"
#pdf(paste(name,".pdf",sep=""),width=8)
xfig(paste(name,".fig",sep=""),width=8)




df <- data.frame(Length=DdeL1$x,Denisovan1_All=DdeL1$y,Denisovan2_African=DdeoAFRL$y,Denisovan3_Asian=DdeoASNL$y,Denisovan4_European=DdeoEURL$y)

ggplot(df, aes(x=Length)) +
    geom_line(aes(y=Denisovan1_All,colour="Denisovan 1All"),size=1) +
    geom_line(aes(y=Denisovan2_African,colour="Denisovan 2African"),size=1) +
    geom_line(aes(y=Denisovan3_Asian,colour="Denisovan 3Asian"),size=1) +
    geom_line(aes(y=Denisovan4_European,colour="Denisovan 4European"),size=1) +
    geom_hline(data=df, yintercept=c(0),linetype="solid", size=0.5) +
    geom_vline(data=df,  aes(xintercept=c(8500,9700,14800,27000),colour=c("Denisovan 1All","Denisovan 2African","Denisovan 4European","Denisovan 3Asian")),linetype="dashed", size=0.6) +
    scale_y_continuous(limits=c(0,0.00005)) +
    opts(title = expression("Denisovan Matching IBD Length Distributions per Exclusive Population")) +
    opts(legend.position=c(.8, .6))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid")) +
    coord_cartesian(ylim =c(0,0.00005),xlim=c(0,60000)) +
    ylab("Density") +
    xlab("IBD Length")


dev.off()


name="IBDlengthDensityNeanderPopOnly"
#pdf(paste(name,".pdf",sep=""),width=8)
xfig(paste(name,".fig",sep=""),width=8)




df <- data.frame(Length=DneL1$x,Neandertal1_All=DneL1$y,Neandertal2_African=DneoAFRL$y,Neandertal3_Asian=DneoASNL$y,Neandertal4_European=DneoEURL$y)

ggplot(df, aes(x=Length)) +
    geom_line(aes(y=Neandertal1_All,colour="Neandertal 1All"),size=1) +
    geom_line(aes(y=Neandertal2_African,colour="Neandertal 2African"),size=1) +
    geom_line(aes(y=Neandertal3_Asian,colour="Neandertal 3Asian"),size=1) +
    geom_line(aes(y=Neandertal4_European,colour="Neandertal 4European"),size=1) +
    geom_hline(data=df, yintercept=c(0),linetype="solid", size=0.5) +
    geom_vline(data=df,  aes(xintercept=c(19500,15500,22800,23800),colour=c("Neandertal 1All","Neandertal 2African","Neandertal 4European","Neandertal 3Asian")),linetype="dashed", size=0.6) +
    scale_y_continuous(limits=c(0,0.00005)) +
    opts(title = expression("Neandertal Matching IBD Length Distributions per Exclusive Population")) +
    opts(legend.position=c(.8, .6))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid")) +
    coord_cartesian(ylim =c(0,0.00005),xlim=c(0,60000)) +
    ylab("Density") +
    xlab("IBD Length")

dev.off()



name="IBDlengthDensityArchaicPopOnly"
#pdf(paste(name,".pdf",sep=""),width=8)
#xfig(paste(name,".fig",sep=""),width=8)



df <- data.frame(Length=DdeneL2$x,Archaic1_All=DdeneL2$y,Archaic2_African=DdeneoAFRL$y,Archaic3_Asian=DdeneoASNL$y,Archaic4_European=DdeneoEURL$y)

ggplot(df, aes(x=Length)) +
    geom_line(aes(y=Archaic1_All,colour="Archaic 1All"),size=1) +
    geom_line(aes(y=Archaic2_African,colour="Archaic 2African"),size=1) +
    geom_line(aes(y=Archaic3_Asian,colour="Archaic 3Asian"),size=1) +
    geom_line(aes(y=Archaic4_European,colour="Archaic 4European"),size=1) +
    geom_hline(data=df, yintercept=c(0),linetype="solid", size=0.5) +
    geom_vline(data=df, aes(xintercept=c(13000,21000,19500,13600),colour=c("Archaic 1All","Archaic 2African","Archaic 3Asian","Archaic 4European")), linetype="dashed", size=0.6) +
    scale_y_continuous(limits=c(0,0.000038)) +
    opts(title = expression("Archaic Matching IBD Length Distributions per Exclusive Population")) +
    opts(legend.position=c(.8, .6))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid")) +
    coord_cartesian(ylim =c(0,0.000038),xlim=c(0,60000)) +
    ylab("Density") +
    xlab("IBD Length")

dev.off()



name="IBDlengthDensityDenisovaAfrica-NoAfrica"
#pdf(paste(name,".pdf",sep=""),width=8)
#xfig(paste(name,".fig",sep=""),width=8)



l1 <- length(deoAFRL)
l2 <- length(deNAFRL)

df <- data.frame(Population = factor( c(rep("Africans",l1),rep("No Africans",l2)) ),
                   proportion = c(deoAFRL,deNAFRL))

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

ggplot(df, aes(x=proportion, fill=Population)) + geom_density(alpha=.3,adjust=0.7) +
    geom_vline(data=cdf, aes(xintercept=c(10000,28000),  colour=Population),
               linetype="dashed", size=1) +
    geom_vline(data=cdf, aes(xintercept=c(10000,20000),  colour=Population),
               linetype="dashed", size=1) +
   opts(title = expression("Denisovan IBD Length Distributions: Africans vs. No Africans")) +
   scale_y_continuous(breaks=c(0,1e-05,2e-05,3e-05,4e-05,5e-05)) +
#    scale_y_continuous(limits=c(0,4e-04)) +
    scale_x_continuous(limits=c(0,60000)) +
    opts(legend.position=c(.8, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid")) +
    coord_cartesian(ylim = c(0,5.4e-05),xlim=c(0,60000)) +
    ylab("Density") +
    xlab("IBD Length")

dev.off()



name="IBDlengthDensityNeandertalAfrica-NoAfrica"
#pdf(paste(name,".pdf",sep=""),width=8)
#xfig(paste(name,".fig",sep=""),width=8)



l1 <- length(neoAFRL)
l2 <- length(neNAFRL)

df <- data.frame(Population = factor( c(rep("Africans",l1),rep("No Africans",l2)) ),
                   proportion = c(neoAFRL,neNAFRL))

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

ggplot(df, aes(x=proportion, fill=Population)) + geom_density(alpha=.3,adjust=0.7) +
    geom_vline(data=cdf, aes(xintercept=c(17000,25800),  colour=Population),
               linetype="dashed", size=1) +
    geom_vline(data=cdf, aes(xintercept=c(17000,21000),  colour=Population),
               linetype="dashed", size=1) +
   opts(title = expression("Neandertal IBD Length Distributions: Africans vs. No Africans")) +
   scale_y_continuous(breaks=c(0,1e-05,2e-05,3e-05,4e-05,5e-05)) +
#    scale_y_continuous(limits=c(0,4e-04)) +
    scale_x_continuous(limits=c(0,60000)) +
    opts(legend.position=c(.8, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid")) +
    coord_cartesian(ylim = c(0,5.4e-05),xlim=c(0,60000)) +
    ylab("Density") +
    xlab("IBD Length")

dev.off()



name="IBDlengthDensityArchaicAfrica-NoAfrica"
#pdf(paste(name,".pdf",sep=""),width=8)
xfig(paste(name,".fig",sep=""),width=8)



l1 <- length(deneoAFRL)
l2 <- length(deneNAFRL)

df <- data.frame(Population = factor( c(rep("Africans",l1),rep("No Africans",l2)) ),
                   proportion = c(deneoAFRL,deneNAFRL))

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

ggplot(df, aes(x=proportion, fill=Population)) + geom_density(alpha=.3,adjust=0.7) +
    geom_vline(data=cdf, aes(xintercept=c(22000,24000),  colour=Population),
               linetype="dashed", size=1) +
    geom_vline(data=cdf, aes(xintercept=c(14300,14700),  colour=Population),
               linetype="dashed", size=1) +
   opts(title = expression("Archaic IBD Length Distributions: Africans vs. No Africans")) +
   scale_y_continuous(breaks=c(0,1e-05,2e-05,3e-05,4e-05)) +
#    scale_y_continuous(limits=c(0,4e-04)) +
    scale_x_continuous(limits=c(0,60000)) +
    opts(legend.position=c(.8, .5))+
    opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid")) +
    coord_cartesian(ylim = c(0,4.3e-05),xlim=c(0,60000)) +
    ylab("Density") +
    xlab("IBD Length")

dev.off()