######################################################################################### # Calculates the number of Neandertal- and Denisovan-matching IBD segments # shared exclusively by a particular continental population for each chromosome and # the number of Neandertal- and Denisovan-matching IBD segments shared exclusively by # Africans that are observed in a particular subpopulation ######################################################################################### source("1000G_pop_functions.R") exPop <- c() obsPop <- c() exAFR <- c() notAFR <- c() denis <- c() neander <- c() for (chr in c(1:22, "X")) { wd <- paste0("chr", chr) setwd(wd) load(file = paste("DenisovaMatNew.Rda",sep = "")) load(file = paste("NeanderMatNew.Rda",sep = "")) load(file = paste("AncDenisovaMatNew.Rda",sep = "")) load(file = paste("AncNeanderMatNew.Rda",sep = "")) load(file = paste("PopulationMat.Rda",sep = "")) load(file = paste("SubPopMat.Rda",sep = "")) real1 <- which(NeanderMat[,"IBDlengthCorrected"] > 10) real2 <- which(NeanderMat[,"#SNVs"] > 7) real <- intersect(real1,real2) DenisovaMat <- DenisovaMat[real,] NeanderMat <- NeanderMat[real,] AncDenisovaMat <- AncDenisovaMat[real,] AncNeanderMat <- AncNeanderMat[real,] PopulationMat <- PopulationMat[real,] SubPopMat <- SubPopMat[real,] de1 <- which(DenisovaMat[,"AltR"] > 0.15) de2 <- which(DenisovaMat[,"ancientRatio"] > 0.3) de <- intersect(de1,de2) ne1 <- which(NeanderMat[,"AltR"] > 0.15) ne2 <- which(NeanderMat[,"ancientRatio"] > 0.3) ne <- intersect(ne1,ne2) #### ancestor correction #### nanca <- which(AncNeanderMat[,"ancestorRatio"] < 0.3) danca <- which(AncDenisovaMat[,"ancestorRatio"] < 0.3) dnanca <- intersect(nanca,danca) de <- intersect(de,danca) ne <- intersect(ne,nanca) ##### get populations #### # only this population: exclusively getoPops(PopulationMat) # not this population: not observed getNPops(PopulationMat) ############ intersections ####################### deoAFR <- intersect(oAFR2,de) deoEAS <- intersect(oEAS2,de) deoSAS <- intersect(oSAS2,de) deoEUR <- intersect(oEUR2,de) neoAFR <- intersect(oAFR2,ne) neoEAS <- intersect(oEAS2,ne) neoSAS <- intersect(oSAS2,ne) neoEUR <- intersect(oEUR2,ne) deNAFR <- intersect(NAFR,de) neNAFR <- intersect(NAFR,ne) #### get subpopulations #### # only this population: exclusively oACB <- which(SubPopMat[,"ACBr"] > 0.9999) oASW <- which(SubPopMat[,"ASWr"] > 0.9999) # at least once this population: observed getJSubPops(SubPopMat) NEASEUR <- intersect(NEAS,NEUR) NEASEURSAS <- intersect(NEASEUR, NSAS) afrd <- which((SubPopMat[,"ESNr"] + SubPopMat[,"GWDr"] + SubPopMat[,"LWKr"] + SubPopMat[,"MSLr"] + SubPopMat[,"YRIr"]) > 0.5) oafr <- intersect(afrd,NEASEURSAS) deoafr <- intersect(oafr,de) neoafr <- intersect(oafr,ne) jafr <- unique(c(JESN,JGWD,JLWK,JMSL,JYRI)) dejafr <- intersect(jafr,de) nejafr <- intersect(jafr,ne) # only AFR and at least once in subpopulation deoAFRJACB <- intersect(deoAFR,JACB) deoAFRJASW <- intersect(deoAFR,JASW) deoAFRJESN <- intersect(deoAFR,JESN) deoAFRJGWD <- intersect(deoAFR,JGWD) deoAFRJLWK <- intersect(deoAFR,JLWK) deoAFRJMSL <- intersect(deoAFR,JMSL) deoAFRJYRI <- intersect(deoAFR,JYRI) neoAFRJACB <- intersect(neoAFR,JACB) neoAFRJASW <- intersect(neoAFR,JASW) neoAFRJESN <- intersect(neoAFR,JESN) neoAFRJGWD <- intersect(neoAFR,JGWD) neoAFRJLWK <- intersect(neoAFR,JLWK) neoAFRJMSL <- intersect(neoAFR,JMSL) neoAFRJYRI <- intersect(neoAFR,JYRI) #### combine results #### exPop <- rbind(exPop, c(length(deoEAS),length(deoSAS),length(deoEUR),length(deoAFR), length(neoEAS),length(neoSAS),length(neoEUR),length(neoAFR))) exAFR <- rbind(exAFR, c(length(deoAFRJACB),length(deoAFRJASW),length(deoAFRJESN),length(deoAFRJGWD), length(deoAFRJLWK),length(deoAFRJMSL),length(deoAFRJYRI),length(neoAFRJACB), length(neoAFRJASW),length(neoAFRJESN),length(neoAFRJGWD),length(neoAFRJLWK), length(neoAFRJMSL),length(neoAFRJYRI))) notAFR <- rbind(notAFR, c(length(NAFR), length(deNAFR), length(neNAFR))) rest <- c(oAMR,oASW,oACB) denis <- c(denis, length(setdiff(de, rest))) neander <- c(neander, length(setdiff(ne, rest))) } colnames(exPop) <- c("deoEAS","deoSAS","deoEUR","deoAFR", "neoEAS","neoSAS","neoEUR","neoAFR") colnames(exAFR) <- c("dACB","dASW","dESN","dGWD","dLWK","dMSL","dYRI", "nACB","nASW","nESN","nGWD","nLWK","nMSL","nYRI") colnames(notAFR) <- c("NAFR", "deNAFR", "neNAFR") exPopFrac <- cbind(exPop[,1:4]/denis, exPop[,5:8]/neander) options(scipen = 999) exPop2 <- rbind(exPop[1:22,], colSums(exPop[1:22,]), exPop[23,]) denis2 <- c(denis[1:22], sum(denis[1:22]), denis[23]) neander2 <- c(neander[1:22], sum(neander[1:22]), neander[23]) exPopFrac2 <- cbind(exPop2[,1:4]/denis2, exPop2[,5:8]/neander2) denisex2 <- cbind(denis2, paste0(exPop2[,1], " (", round((exPopFrac2[,1]*100), digits = 1), "%)"), paste0(exPop2[,2], " (", round((exPopFrac2[,2]*100), digits = 1), "%)") , paste0(exPop2[,3], " (", round((exPopFrac2[,3]*100), digits = 1), "%)"), paste0(exPop2[,4], " (", round((exPopFrac2[,4]*100), digits = 1), "%)")) neanderex2 <- cbind(neander2, paste0(exPop2[,5], " (", round(exPopFrac2[,5]*100, digits = 1), "%)") , paste0(exPop2[,6], " (", round(exPopFrac2[,6]*100, digits = 1), "%)"), paste0(exPop2[,7], " (", round(exPopFrac2[,7]*100, digits = 1), "%)"), paste0(exPop2[,8], " (", round(exPopFrac2[,8]*100, digits = 1), "%)")) neanderdenisex2 <- cbind(neanderex2, denisex2) colnames(neanderdenisex2) <- c("Neandertal", "EAS", "SAS", "EUR", "AFR", "Denisova", "EAS", "SAS", "EUR", "AFR") write.csv(format(neanderdenisex2, big.mark = ","), file = "neanderdenisex.csv") write.csv(format(exAFR, big.mark = ","), file = "exAFR.csv") denisexAFR <- rbind(c(exPop2[23,4], colSums(exAFR[1:22,1:7])), c(exPop2[24,4], exAFR[23,1:7])) neanderexAFR <- rbind(c(exPop2[23,8], colSums(exAFR[1:22,8:14])), c(exPop2[24,8], exAFR[23,8:14])) neanderdenisexAFR <- rbind(neanderexAFR, denisexAFR) rownames(neanderdenisexAFR) <- c("Neandertal 1-22", "Neandertal X", "Denisova 1-22", "Denisova X") colnames(neanderdenisexAFR) <- c("AFR", "ACB", "ASW", "ESN", "GWD", "LWK", "MSL", "YRI") write.csv(format(neanderdenisexAFR, big.mark = ","), file = "neanderdenisexAFR.csv") #