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("/seppdata/sepp/linkage/release/scripts/funct.R") library(ggplot2) library(doBy) library(fabia) cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7") summary(DenisovaMat[,21]) # paper densiova/neanderthal -> >= 0.3 # population >= 0.2 de <- which(DenisovaMat[,12]>0.3) #11561 DE <- which(DenisovaMat[,12]>0.6) #2562 ne <- which(NeanderMat[,12]>0.3) NE <- which(NeanderMat[,12]>0.6) 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,21] deL <- DenisovaMat[de,21] NEL <- NeanderMat[NE,21] DEL <- DenisovaMat[DE,21] LLL <- DenisovaMat[,21] LLLAFR <- DenisovaMat[AFR,21] LLLASN <- DenisovaMat[ASN,21] LLLEUR <- DenisovaMat[EUR,21] LLLAFRd <- DenisovaMat[AFRd,21] LLLASNd <- DenisovaMat[ASNd,21] LLLEURd <- DenisovaMat[EURd,21] DEL <- DenisovaMat[DE,21] DLLL <- density(LLL,bw=700,from=0,to=80000) DLLL1 <- density(LLL,bw=700,from=0,to=5000) DLLL2 <- density(LLL,bw=700,from=15000,to=80000) DLLLAFR <- density(LLLAFR,bw=700,from=0,to=80000) DLLLAFR1 <- density(LLLAFR,bw=700,from=0,to=5000) DLLLAFR2 <- density(LLLAFR,bw=700,from=15000,to=80000) DLLLEUR <- density(LLLEUR,bw=700,from=0,to=80000) DLLLEUR1 <- density(LLLEUR,bw=700,from=0,to=5000) DLLLEUR2 <- density(LLLEUR,bw=700,from=15000,to=80000) DLLLASN <- density(LLLASN,bw=700,from=0,to=80000) DLLLASN1 <- density(LLLASN,bw=700,from=0,to=5000) DLLLASN2 <- density(LLLASN,bw=700,from=15000,to=80000) DLLLAFRd <- density(LLLAFRd,bw=700,from=0,to=80000) DLLLAFR1d <- density(LLLAFRd,bw=700,from=0,to=5000) DLLLAFR2d <- density(LLLAFRd,bw=700,from=15000,to=80000) DLLLEURd <- density(LLLEURd,bw=700,from=0,to=80000) DLLLEUR1d <- density(LLLEURd,bw=700,from=0,to=5000) DLLLEUR2d <- density(LLLEURd,bw=700,from=15000,to=80000) DLLLASNd <- density(LLLASNd,bw=700,from=0,to=80000) DLLLASN1d <- density(LLLASNd,bw=700,from=0,to=5000) DLLLASN2d <- density(LLLASNd,bw=700,from=15000,to=80000) DDEL <- density(DEL,bw=700,from=0,to=80000) DEASN <- intersect(ASN,DE) DEASNL <- DenisovaMat[DEASN,21] DDEASNL <- density(DEASNL,bw=700,from=0,to=80000) DEAFR <- intersect(AFR,DE) DEAFRL <- DenisovaMat[DEAFR,21] DDEAFRL <- density(DEAFRL,bw=700,from=0,to=80000) DEEUR <- intersect(EUR,DE) DEEURL <- DenisovaMat[DEEUR,21] DDEEURL <- density(DEEURL,bw=700,from=0,to=80000) DEASNd <- intersect(ASNd,DE) # 392 DEASNLd <- DenisovaMat[DEASNd,21] DDEASNLd <- density(DEASNLd,bw=700,from=0,to=80000) DEAFRd <- intersect(AFRd,DE) DEAFRLd <- DenisovaMat[DEAFRd,21] DDEAFRLd <- density(DEAFRLd,bw=700,from=0,to=80000) DEEURd <- intersect(EURd,DE) # 206 DEEURLd <- DenisovaMat[DEEURd,21] DDEEURLd <- density(DEEURLd,bw=700,from=0,to=80000) DdeL <- density(deL,bw=700,from=0,to=80000) DdeL1 <- density(deL,bw=700,from=0,to=5000) DdeL2 <- density(deL,bw=700,from=15000,to=80000) deASN <- intersect(ASN,de) deASNL <- DenisovaMat[deASN,21] DdeASNL <- density(deASNL,bw=700,from=0,to=80000) deAFR <- intersect(AFR,de) deAFRL <- DenisovaMat[deAFR,21] DdeAFRL <- density(deAFRL,bw=700,from=0,to=80000) deEUR <- intersect(EUR,de) deEURL <- DenisovaMat[deEUR,21] DdeEURL <- density(deEURL,bw=700,from=0,to=80000) deASNd <- intersect(ASNd,de) # 1303 deASNLd <- DenisovaMat[deASNd,21] DdeASNLd <- density(deASNLd,bw=700,from=0,to=80000) deAFRd <- intersect(AFRd,de) deAFRLd <- DenisovaMat[deAFRd,21] DdeAFRLd <- density(deAFRLd,bw=700,from=0,to=80000) deEURd <- intersect(EURd,de) #886 deEURLd <- DenisovaMat[deEURd,21] DdeEURLd <- density(deEURLd,bw=700,from=0,to=80000) ################################################## LLL <- NeanderMat[,21] LLLAFR <- NeanderMat[AFR,21] LLLASN <- NeanderMat[ASN,21] LLLEUR <- NeanderMat[EUR,21] LLLAFRd <- NeanderMat[AFRd,21] LLLASNd <- NeanderMat[ASNd,21] LLLEURd <- NeanderMat[EURd,21] DLLL <- density(LLL,bw=700,from=0,to=80000) DLLL1 <- density(LLL,bw=700,from=0,to=5000) DLLL2 <- density(LLL,bw=700,from=15000,to=80000) DLLLAFR <- density(LLLAFR,bw=700,from=0,to=80000) DLLLAFR1 <- density(LLLAFR,bw=700,from=0,to=5000) DLLLAFR2 <- density(LLLAFR,bw=700,from=15000,to=80000) DLLLEUR <- density(LLLEUR,bw=700,from=0,to=80000) DLLLEUR1 <- density(LLLEUR,bw=700,from=0,to=5000) DLLLEUR2 <- density(LLLEUR,bw=700,from=15000,to=80000) DLLLASN <- density(LLLASN,bw=700,from=0,to=80000) DLLLASN1 <- density(LLLASN,bw=700,from=0,to=5000) DLLLASN2 <- density(LLLASN,bw=700,from=15000,to=80000) DLLLAFRd <- density(LLLAFRd,bw=700,from=0,to=80000) DLLLAFR1d <- density(LLLAFRd,bw=700,from=0,to=5000) DLLLAFR2d <- density(LLLAFRd,bw=700,from=15000,to=80000) DLLLEURd <- density(LLLEURd,bw=700,from=0,to=80000) DLLLEUR1d <- density(LLLEURd,bw=700,from=0,to=5000) DLLLEUR2d <- density(LLLEURd,bw=700,from=15000,to=80000) DLLLASNd <- density(LLLASNd,bw=700,from=0,to=80000) DLLLASN1d <- density(LLLASNd,bw=700,from=0,to=5000) DLLLASN2d <- density(LLLASNd,bw=700,from=15000,to=80000) NEL <- NeanderMat[NE,21] DNEL <- density(NEL,bw=700,from=0,to=80000) NEASN <- intersect(ASN,NE) NEASNL <- NeanderMat[NEASN,21] DNEASNL <- density(NEASNL,bw=700,from=0,to=80000) NEAFR <- intersect(AFR,NE) NEAFRL <- NeanderMat[NEAFR,21] DNEAFRL <- density(NEAFRL,bw=700,from=0,to=80000) NEEUR <- intersect(EUR,NE) NEEURL <- NeanderMat[NEEUR,21] DNEEURL <- density(NEEURL,bw=700,from=0,to=80000) NEASNd <- intersect(ASNd,NE) # 392 NEASNLd <- NeanderMat[NEASNd,21] DNEASNLd <- density(NEASNLd,bw=700,from=0,to=80000) NEAFRd <- intersect(AFRd,NE) NEAFRLd <- NeanderMat[NEAFRd,21] DNEAFRLd <- density(NEAFRLd,bw=700,from=0,to=80000) NEEURd <- intersect(EURd,NE) # 206 NEEURLd <- NeanderMat[NEEURd,21] DNEEURLd <- density(NEEURLd,bw=700,from=0,to=80000) DneL <- density(neL,bw=700,from=0,to=80000) DneL1 <- density(neL,bw=700,from=0,to=5000) DneL2 <- density(neL,bw=700,from=15000,to=80000) neASN <- intersect(ASN,ne) neASNL <- NeanderMat[neASN,21] DneASNL <- density(neASNL,bw=700,from=0,to=80000) neAFR <- intersect(AFR,ne) neAFRL <- NeanderMat[neAFR,21] DneAFRL <- density(neAFRL,bw=700,from=0,to=80000) neEUR <- intersect(EUR,ne) neEURL <- NeanderMat[neEUR,21] DneEURL <- density(neEURL,bw=700,from=0,to=80000) neASNd <- intersect(ASNd,ne) # 1303 neASNLd <- NeanderMat[neASNd,21] DneASNLd <- density(neASNLd,bw=700,from=0,to=80000) neAFRd <- intersect(AFRd,ne) neAFRLd <- NeanderMat[neAFRd,21] DneAFRLd <- density(neAFRLd,bw=700,from=0,to=80000) neEURd <- intersect(EURd,ne) #886 neEURLd <- NeanderMat[neEURd,21] DneEURLd <- density(neEURLd,bw=700,from=0,to=80000) ################################################ ################################################ ################################################ ################################################ ################################################ #DENISOVA LENGTH to all length ################################################ #Histogramm length x11() hist(LLL,breaks=500,ylim=c(0,5500),xlim=c(0,80000),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.. < 6000),colour="darkgreen",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)) + scale_x_continuous(limits=c(0,80000)) + geom_vline(xintercept=c(24000), 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,5500),xlim=c(0,80000)) + ylab("Counts") + xlab("IBD Length") dev.off() #histogramm length denisova x11() hist(deL,breaks=500,ylim=c(0,300),xlim=c(0,80000),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",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)) + scale_x_continuous(limits=c(0,80000)) + geom_vline(xintercept=c(3100,24000), 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,310),xlim=c(0,80000)) + 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.3) + geom_vline(data=cdf, aes(xintercept=c(6400,26000), colour=IBD_Length), linetype="dashed", size=1) + geom_vline(data=cdf, aes(xintercept=c(39500,26000), 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,80000)) + opts(legend.position=c(.8, .5))+ opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+ coord_cartesian(ylim = c(0,5.5e-05),xlim=c(0,80000)) + 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,80000)) + 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,80000)) + 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) l1 <- length(DLLL) df <- data.frame(Length=DLLL$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=cdf, aes(xintercept=c(26000,28000,43000,57500),colour=c("Denisovan African","Denisovan Asian","Denisovan European","Denisovan Asian")),linetype="dashed", size=0.6) + scale_y_continuous(limits=c(0,0.00005)) + opts(title = expression("IBD Length Distributions")) + 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,80000)) + 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,80000),plot=TRUE,col="blue") abline(v=24000) #histogramm length neander x11() hist(neL,breaks=500,ylim=c(0,350),xlim=c(0,80000),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",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)) + scale_x_continuous(limits=c(0,80000)) + geom_vline(xintercept=c(3100,24000), 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,360),xlim=c(0,80000)) + 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.3) + geom_vline(data=cdf, aes(xintercept=c(7000,26700),color=IBD_Length), linetype="dashed", size=1) + geom_vline(data=cdf, aes(xintercept=c(44000,26700),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,80000)) + opts(legend.position=c(.8, .5))+ opts(legend.background = theme_rect(fill="gray90", size=0.3,linetype="solid"))+ coord_cartesian(ylim = c(0,5.5e-05),xlim=c(0,80000)) + 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(7000,23200,44000)),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,80000)) + 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,80000)) + 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) l1 <- length(DLLL) df <- data.frame(Length=DLLL$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=cdf, aes(xintercept=c(24800,27000,43500),colour=c("Neandertal African","Neandertal Asian","Neandertal European")), linetype="dashed", size=0.6) + scale_y_continuous(limits=c(0,0.00005)) + opts(title = expression("IBD Length Distributions")) + 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,80000)) + 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,80000),plot=TRUE,col="red") x11() hist(LLLAFR,breaks=500,ylim=c(0,5000),xlim=c(0,80000),plot=TRUE,col="blue") x11() hist(LLLEUR,breaks=500,ylim=c(0,300),xlim=c(0,80000),plot=TRUE,col="green") x11() hist(deASNL,breaks=200,ylim=c(0,40),xlim=c(0,80000),plot=TRUE,col="red") x11() hist(deEURL,breaks=200,ylim=c(0,30),xlim=c(0,80000),plot=TRUE,col="green") x11() hist(deAFRL,breaks=500,ylim=c(0,300),xlim=c(0,80000),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,80000),plot=TRUE,col="blue") hist(DEEURL[which(DEEURL>1000)],breaks=100,xlim=c(0,80000),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,80000),plot=TRUE,col="blue") hist(deASNL[which(deASNL>1000)],breaks=100,xlim=c(0,80000),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,21] plot(density(DEASNL,bw=700,from=0,to=80000)) DDEASNL <- density(DEASNL,bw=700,from=0,to=80000) plot(DLLL$x,(DDEASNL$y-DLLL$y),type="l") abline(h=0) DEAFR <- intersect(AFR,DE) DEAFRL <- DenisovaMat[DEAFR,21] plot(density(DEAFRL,bw=700,from=0,to=80000)) DDEAFRL <- density(DEAFRL,bw=700,from=0,to=80000) plot(DLLL$x,(DDEAFRL$y-DLLL$y),type="l") abline(h=0) DEEUR <- intersect(EUR,DE) DEEURL <- DenisovaMat[DEEUR,21] plot(density(DEEURL,bw=700,from=0,to=80000)) DDEEURL <- density(DEEURL,bw=700,from=0,to=80000) plot(DLLL$x,(DDEEURL$y-DLLL$y),type="l") abline(h=0) DDEEURL1 <- density(DEEURL,bw=700,from=0,to=5000) plot(DLLL1$x,(DDEEURL1$y-DLLL1$y),type="l") abline(h=0) DDEEURL2 <- density(DEEURL,bw=700,from=15000,to=80000) plot(DLLL2$x,(DDEEURL2$y-DLLL2$y),type="l") abline(h=0) neL <- NeanderMat[ne,21] plot(density(neL)) hist(neL,breaks=10000,xlim=c(0,80000),plot=TRUE) hist(neL,breaks=10000,xlim=c(0,1000),plot=TRUE) plot(density(neL,bw=700,from=20000,to=30000)) neL <- NeanderMat[ne,21] deL <- DenisovaMat[de,21] plot(density(deL)) plot(density(deL,bw=700,from=0,to=80000)) plot(density(deL,bw=700,from=20000,to=30000)) > hist(deL,breaks=10000,xlim=c(0,80000),plot=TRUE) > hist(deL,breaks=10000,xlim=c(0,1000),plot=TRUE) neL <- NeanderMat[ne,21] plot(density(neL)) hist(neL,breaks=10000,xlim=c(0,80000),plot=TRUE) hist(neL,breaks=10000,xlim=c(0,1000),plot=TRUE) plot(density(neL,bw=700,from=20000,to=30000)) deASN <- intersect(ASN,de) deASNL <- DenisovaMat[deASN,21] plot(density(deASNL)) neASN <- intersect(ASN,ne) neASNL <- NeanderMat[neASN,21] plot(density(neASNL,bw=700,from=0,to=80000)) deAFR <- intersect(AFR,de) deAFRL <- DenisovaMat[deAFR,21] plot(density(deAFRL,bw=700,from=0,to=80000)) neAFR <- intersect(AFR,ne) neAFRL <- NeanderMat[neAFR,21] plot(density(neAFRL,bw=700,from=0,to=80000)) deEUR <- intersect(EUR,de) deEURL <- DenisovaMat[deEUR,21] plot(density(deEURL,bw=700,from=0,to=80000)) neEUR <- intersect(EUR,ne) neEURL <- NeanderMat[neEUR,21] plot(density(neEURL,bw=700,from=0,to=80000)) 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=700,from=0,to=80000)) points(density(DEAFRL,bw=700,from=0,to=80000),type="l",col="blue") points(density(DEASNL,bw=700,from=0,to=80000),type="l",col="red") points(density(DEEURL,bw=700,from=0,to=80000),type="l",col="green") plot(density(DEL,bw=700,from=0,to=5000)) points(density(DEAFRL,bw=700,from=0,to=5000),type="l",col="blue") points(density(DEASNL,bw=700,from=0,to=5000),type="l",col="red") points(density(DEEURL,bw=700,from=0,to=5000),type="l",col="green") plot(density(DEL,bw=700,from=15000,to=35000)) points(density(DEAFRL,bw=700,from=15000,to=35000),type="l",col="blue") points(density(DEASNL,bw=700,from=15000,to=35000),type="l",col="red") points(density(DEEURL,bw=700,from=15000,to=35000),type="l",col="green") plot(density(deL,bw=700,from=0,to=80000)) points(density(deAFRL,bw=700,from=0,to=80000),type="l",col="blue") points(density(deASNL,bw=700,from=0,to=80000),type="l",col="red") points(density(deEURL,bw=700,from=0,to=80000),type="l",col="green") plot(density(deL,bw=700,from=0,to=5000)) points(density(deAFRL,bw=700,from=0,to=5000),type="l",col="blue") points(density(deASNL,bw=700,from=0,to=5000),type="l",col="red") points(density(deEURL,bw=700,from=0,to=5000),type="l",col="green") plot(density(deL,bw=700,from=15000,to=35000)) points(density(deAFRL,bw=700,from=15000,to=35000),type="l",col="blue") points(density(deASNL,bw=700,from=15000,to=35000),type="l",col="red") points(density(deEURL,bw=700,from=15000,to=35000),type="l",col="green") plot(density(DEL,bw=700,from=0,to=80000)) points(density(DEAFRL,bw=700,from=0,to=80000),type="l",col="blue") points(density(DEASNL,bw=700,from=0,to=80000),type="l",col="red") points(density(DEEURL,bw=700,from=0,to=80000),type="l",col="green") hist(DEL,breaks=100,xlim=c(0,80000),plot=TRUE) hist(DEAFRL,breaks=100,xlim=c(0,80000),plot=TRUE,col="blue") hist(DEASNL,breaks=100,xlim=c(0,80000),plot=TRUE,col="red") hist(DEEURL,breaks=100,xlim=c(0,80000),plot=TRUE,col="green") hist(deL,breaks=1000,xlim=c(0,80000),plot=TRUE) hist(deAFRL,breaks=1000,xlim=c(0,80000),plot=TRUE,col="blue") hist(deASNL,breaks=400,xlim=c(0,80000),plot=TRUE,col="red") hist(deEURL,breaks=400,xlim=c(0,80000),plot=TRUE,col="green") ######################################################################### ######################################################################### ######################################################################### summary(DenisovaMat[,21]) summary(DenisovaMat[DE,21]) summary(DenisovaMat[NE,21]) summary(DenisovaMat[de,21]) summary(DenisovaMat[ne,21]) > summary(DenisovaMat[,21]) Min. 1st Qu. Median Mean 3rd Qu. Max. 56 20120 25380 26970 30550 21250000 > summary(DenisovaMat[DE,21]) Min. 1st Qu. Median Mean 3rd Qu. Max. 252 17800 24650 24540 30080 158100 > summary(DenisovaMat[NE,21]) Min. 1st Qu. Median Mean 3rd Qu. Max. 509 19390 25600 26170 31370 201400 > summary(DenisovaMat[de,21]) Min. 1st Qu. Median Mean 3rd Qu. Max. 252 18750 25560 26460 31740 373300 > summary(DenisovaMat[ne,21]) Min. 1st Qu. Median Mean 3rd Qu. Max. 509 19760 25980 26960 32150 535500