load(file=paste("DenisovaMat.Rda",sep="")) load(file=paste("NeanderMat.Rda",sep="")) load(file=paste("AncestMat.Rda",sep="")) load(file=paste("PopulationMat.Rda",sep="")) library(ggplot2) library(doBy) cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7") 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) ################################################## LLL <- NeanderMat[,22] LLLAFR <- NeanderMat[AFR,22] LLLASN <- NeanderMat[ASN,22] LLLEUR <- NeanderMat[EUR,22] LLLAFRd <- NeanderMat[AFRd,22] LLLASNd <- NeanderMat[ASNd,22] LLLEURd <- NeanderMat[EURd,22] 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) 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 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) + labs(title = "Histogram IBD Lengths") + theme(legend.position="none")+ # theme(legend.background = element_rect(fill="gray90", size=0.3,linetype="solid",colour="black"))+ coord_cartesian(ylim = c(0,8300),xlim=c(0,60000)) + ylab("Counts") + xlab("IBD Length") dev.off() #histogramm length denisova 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) + labs(title = "Histogram Denisovan IBD Lengths") + theme(legend.position="none")+ # theme(legend.background = element_rect(fill="gray90", size=0.3,linetype="solid",colour="black"))+ coord_cartesian(ylim = c(0,220),xlim=c(0,60000)) + ylab("Counts") + xlab("IBD Length") dev.off() #Length densities 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.2) + labs(title = "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)) + theme(legend.position=c(.8, .5))+ theme(legend.background = element_rect(fill="gray90", size=0.3,linetype="solid",colour="black")) + coord_cartesian(ylim = c(0,5.2e-05),xlim=c(0,60000)) + ylab("Density") + xlab("IBD Length") dev.off() #difference of length densities 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.5) + scale_fill_manual(values=c(cbPalette[8],cbPalette[1])) + scale_colour_manual(values=c(cbPalette[4],cbPalette[1])) + geom_vline(data=df, xintercept=c(6400,24500,39500),linetype="dashed", size=0.5) + # geom_hline(data=df, aes(yintercept=c(0),colour="green"),linetype="solid", size=1.2) + labs(title = "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)) + theme(legend.position="none")+ # theme(legend.background = element_rect(fill="gray90", size=0.3,linetype="solid",colour="black"))+ 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 ################################################ 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.2) + geom_line(aes(y=Denisovan2_African,colour="Denisovan African"),size=1.2) + geom_line(aes(y=Denisovan3_Asian,colour="Denisovan Asian"),size=1.2) + geom_line(aes(y=Denisovan4_European,colour="Denisovan European"),size=1.2) + geom_hline(data=df, yintercept=c(0),linetype="solid", size=0.5) + geom_vline(data=df, aes(xintercept=c(8500,18000,24400,25500),colour=c("Denisovan African","Denisovan European","Denisovan Asian","Denisovan European")),linetype="dashed", size=0.8) + geom_vline(data=df, aes(xintercept=c(40000),colour=c("Denisovan European")),linetype="dashed", size=0.8) + scale_y_continuous(limits=c(0,0.000053)) + labs(title = "Denisovan Matching IBD Length Distributions per Population") + theme(legend.position=c(.8, .6))+ theme(legend.background = element_rect(fill="gray90", size=0.3,linetype="solid",colour="black"))+ 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 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) + labs(title = "Histogram Neandertal IBD Lengths") + theme(legend.position="none")+ # theme(legend.background = element_rect(fill="gray90", size=0.3,linetype="solid",colour="black"))+ coord_cartesian(ylim = c(0,350),xlim=c(0,60000)) + ylab("Counts") + xlab("IBD Length") dev.off() #Length densities 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.2) + geom_vline(data=cdf, aes(xintercept=c(43000,24200),color=IBD_Length), linetype="dashed", size=1.2) + labs(title = "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)) + theme(legend.position=c(.8, .5))+ theme(legend.background = element_rect(fill="gray90", size=0.3,linetype="solid",colour="black"))+ coord_cartesian(ylim = c(0,5.2e-05),xlim=c(0,60000)) + ylab("Density") + xlab("IBD Length") dev.off() #difference of length densities 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.5) + scale_fill_manual(values=c(cbPalette[8],cbPalette[1])) + scale_colour_manual(values=c(cbPalette[4],cbPalette[1])) + geom_vline(data=df, xintercept=c(11000,24200,43000),linetype="dashed", size=0.5) + # geom_hline(data=df, aes(yintercept=c(0),colour="green"),linetype="solid", size=1.2) + labs(title = "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)) + theme(legend.position="none")+ # theme(legend.background = element_rect(fill="gray90", size=0.3,linetype="solid",colour="black"))+ 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 ################################################ 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.2) + geom_line(aes(y=Neandertal2_African,colour="Neandertal African"),size=1.2) + geom_line(aes(y=Neandertal3_Asian,colour="Neandertal Asian"),size=1.2) + geom_line(aes(y=Neandertal4_European,colour="Neandertal European"),size=1.2) + 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.8) + scale_y_continuous(limits=c(0,0.00005)) + labs(title = "Neandertal Matching IBD Length Distributions per Population") + theme(legend.position=c(.8, .6))+ theme(legend.background = element_rect(fill="gray90", size=0.3,linetype="solid",colour="black"))+ 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 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) + labs(title = "Histogram Archaic IBD Lengths") + theme(legend.position="none")+ # theme(legend.background = element_rect(fill="gray90", size=0.3,linetype="solid",colour="black"))+ coord_cartesian(ylim = c(0,230),xlim=c(0,60000)) + ylab("Counts") + xlab("IBD Length") dev.off() #Length densities 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.2) + geom_vline(data=cdf, aes(xintercept=c(41000,24200),color=IBD_Length), linetype="dashed", size=1.2) + labs(title = "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)) + theme(legend.position=c(.8, .5))+ theme(legend.background = element_rect(fill="gray90", size=0.3,linetype="solid",colour="black"))+ coord_cartesian(ylim = c(0,5.2e-05),xlim=c(0,60000)) + ylab("Density") + xlab("IBD Length") dev.off() #difference of length densities 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.5) + scale_fill_manual(values=c(cbPalette[8],cbPalette[1])) + scale_colour_manual(values=c(cbPalette[4],cbPalette[1])) + geom_vline(data=df, xintercept=c(11200,24200,42000),linetype="dashed", size=0.5) + # geom_hline(data=df, aes(yintercept=c(0),colour="green"),linetype="solid", size=1.2) + labs(title = "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)) + theme(legend.position="none")+ # theme(legend.background = element_rect(fill="gray90", size=0.3,linetype="solid",colour="black"))+ 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 ################################################ 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.2) + geom_line(aes(y=Archaic2_African,colour="Archaic African"),size=1.2) + geom_line(aes(y=Archaic3_Asian,colour="Archaic Asian"),size=1.2) + geom_line(aes(y=Archaic4_European,colour="Archaic European"),size=1.2) + 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.8) + scale_y_continuous(limits=c(0,0.000038)) + labs(title = "Archaic Matching IBD Length Distributions per Population") + theme(legend.position=c(.8, .6))+ theme(legend.background = element_rect(fill="gray90", size=0.3,linetype="solid",colour="black"))+ coord_cartesian(ylim =c(0,0.000038),xlim=c(0,60000)) + ylab("Density") + xlab("IBD Length") dev.off() ################################################# ################################################ ################################################ ################################################ ################################################ ################################################ ############################################### ################################################ ################################################ ################################################ 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.2) + geom_line(aes(y=Denisovan2_African,colour="Denisovan 2African"),size=1.2) + geom_line(aes(y=Denisovan3_Asian,colour="Denisovan 3Asian"),size=1.2) + geom_line(aes(y=Denisovan4_European,colour="Denisovan 4European"),size=1.2) + 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.8) + scale_y_continuous(limits=c(0,0.00005)) + labs(title = "Denisovan Matching IBD Length Distributions per Exclusive Population") + theme(legend.position=c(.8, .6))+ theme(legend.background = element_rect(fill="gray90", size=0.3,linetype="solid",colour="black"))+ 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.2) + geom_line(aes(y=Neandertal2_African,colour="Neandertal 2African"),size=1.2) + geom_line(aes(y=Neandertal3_Asian,colour="Neandertal 3Asian"),size=1.2) + geom_line(aes(y=Neandertal4_European,colour="Neandertal 4European"),size=1.2) + 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.8) + scale_y_continuous(limits=c(0,0.00005)) + labs(title = "Neandertal Matching IBD Length Distributions per Exclusive Population") + theme(legend.position=c(.8, .6))+ theme(legend.background = element_rect(fill="gray90", size=0.3,linetype="solid",colour="black"))+ 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.2) + geom_line(aes(y=Archaic2_African,colour="Archaic 2African"),size=1.2) + geom_line(aes(y=Archaic3_Asian,colour="Archaic 3Asian"),size=1.2) + geom_line(aes(y=Archaic4_European,colour="Archaic 4European"),size=1.2) + 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 1All","Archaic 2African","Archaic 3Asian","Archaic 4European")), linetype="dashed", size=0.8) + scale_y_continuous(limits=c(0,0.000038)) + labs(title = "Archaic Matching IBD Length Distributions per Exclusive Population") + theme(legend.position=c(.8, .6))+ theme(legend.background = element_rect(fill="gray90", size=0.3,linetype="solid",colour="black"))+ 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.2) + geom_vline(data=cdf, aes(xintercept=c(10000,20000), colour=Population), linetype="dashed", size=1.2) + labs(title = "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)) + theme(legend.position=c(.8, .5))+ theme(legend.background = element_rect(fill="gray90", size=0.3,linetype="solid",colour="black")) + 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.2) + geom_vline(data=cdf, aes(xintercept=c(17000,21000), colour=Population), linetype="dashed", size=1.2) + labs(title = "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)) + theme(legend.position=c(.8, .5))+ theme(legend.background = element_rect(fill="gray90", size=0.3,linetype="solid",colour="black")) + 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.2) + geom_vline(data=cdf, aes(xintercept=c(14300,14700), colour=Population), linetype="dashed", size=1.2) + labs(title = "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)) + theme(legend.position=c(.8, .5))+ theme(legend.background = element_rect(fill="gray90", size=0.3,linetype="solid",colour="black")) + coord_cartesian(ylim = c(0,4.3e-05),xlim=c(0,60000)) + ylab("Density") + xlab("IBD Length") dev.off()