require(ggplot2) require(DESeq2) require(compcodeR) ### Examples ##simulation II #run.experiments.tax.glo(experiment = 'sim2', nIter = 3, do.ROC = TRUE, do.plot = TRUE) ##real data #run.experiments.tax.glo(experiment = 'real', file = '/home/myfile.csv', condA = '13days', condB = '15days') run.experiments.tax.glo <- function(experiment= 'none', file = NULL, nIter = 10, do.ROC = FALSE, condA = NULL, condB = NULL, do.plot = TRUE, pval = 0.05, returnFLAG = FALSE) { # experiment: The type of experiment to perform # file: The csv file provided in the supplemental materials # nIter: Number of iterations for the simulations # do.ROC: Should a ROC curve be plotted # condX: What conditions should be compared for the real data # do.plot: Produce boxplots for the simulations # pval: pvalue threshold # returnFLAG: return the results as a list of lists # returns a list for each iteration with the following entries: #"T.Mat": taxon-specific scaling: nOrg * comparison label; DEF, TP DEF, opposite direction DEF, NDE, FP DEF (are NDE), missed DEF #"G.Mat": global scaling: nOrg * comparison label; DEF, TP DEF, opposite direction DEF, NDE, FP DEF (are NDE), missed DEF #"res.tax.padj" taxon-specific scaling adjusted p-values (BH) #"res.glo.padj" global scaling adjusted p-values (BH)_ #"res.tax.log2FoldChange" taxon-specific scaling log2 fold changes #"res.glo.log2FoldChange" global scaling log2 fold changes #"UP.reg" features generated as up-regulated DEF #"DOWN.reg" features generated as down-regulated DEF set.seed(42); CHOICES = c('sim1','sim2','sim3','sim4','real') Res <- list() if (experiment %in% CHOICES) { } else { cat('please select one of the following experiments to run:',CHOICES,'\n') return() } if (experiment %in% CHOICES[1:3]) { for (i in 1:nIter) { Res[[i]] = run.sim.tax.glo(experiment, pval) } } if (experiment == CHOICES[4]) { nCut.glo = vector(mode = 'numeric', length = nIter) nCut.tax = vector(mode = 'numeric', length = nIter) for (i in 1:nIter) { Res[[i]] = run.sim.tax.glo4(experiment, pval) } for (i in 1:nIter) { nCut.tax[i] = all.cuts(X.Stacked = Res[[i]]$X.stacked, Y.Stacked = Res[[i]]$Y.stacked, X.Super = Res[[i]]$X.super.tax, Y.Super = Res[[i]]$Y.super.tax, UP.Vec = Res[[i]]$UP.reg, Down.Vec = Res[[1]]$DOWN.reg) nCut.glo[i] = all.cuts(X.Stacked = Res[[i]]$X.stacked, Y.Stacked = Res[[i]]$Y.stacked, X.Super = Res[[i]]$X.super.glo, Y.Super = Res[[i]]$Y.super.glo, UP.Vec = Res[[i]]$UP.reg, Down.Vec = Res[[1]]$DOWN.reg) } cat('average cuts with diagonal global: ',mean(nCut.glo),'\n') cat('average cuts with diagonal taxon-specific: ',mean(nCut.tax),'\n') } if (experiment == CHOICES[5]) { if (is.null(file)) { cat('please provide the file for analysis.\n') return(-1) } if (is.null(condA) | is.null(condB)) { cat('no conditions selected, using default condition 13days vs 27days\n') cat('availible conditions are:\n') cat(c("13days", "15days", "16days", "27days", "29days", "30days"),'\n') condA = '13days' condB = '27days' } ret = real.data(file, condA, condB, pval) return() } if (do.plot) { if (experiment %in% CHOICES[1:3]) { df = generate.plot.data(Result.List = Res) df$method = factor(df$method, c('glo','tax')) x11() g <- ggplot(df, aes(x = factor(method), y = x)) g = g + geom_boxplot(aes(fill = factor(Org))) + facet_grid(type~., scales = 'free_y') + xlab('normalization method') + ylab('#features') print(g) } if (experiment == CHOICES[4]) { X = plot.prepare.sim4(Results = Res, dir = 'UP') for (i in 1:nIter) { x11() print(X$arrow.plot[[i]]) x11() print(X$some.plot[[i]]) } } } if(do.ROC & experiment %in% CHOICES[1:3]) { generate.roc(Res,nIter) } if (returnFLAG) { return(Res) } } run.sim.tax.glo <- function(experiment, p.thresh) { require(compcodeR) require(DESeq2) nOrg = 5 #Simulation: number of organism to simulate nSamples = 6 #Simulation: number of samples per condition CountVec = c(1e7,5e6,1e6,5e5,1e5) #Simulation: base library size Sampling.Rate.Mat = matrix(1,ncol = nOrg, nrow = nSamples*2) #Simulation: sampling rates for each organism and sample cond.vec <- c(rep('A',nSamples),rep('B',nSamples)) #conditions for each sample type.vec <- rep('1',2*nSamples) #type of read UP.reg = c() #features labeled as up-regulated by the data generation process DOWN.reg = c() #features labeled as downr-egulated by the data generation process if (experiment == 'sim2') { for (i in 1:nOrg) { Sampling.Rate.Mat[,i] = sample(seq(0.5,2,0.1), size = nSamples*2, replace = TRUE) } } if (experiment == 'sim3') { for (i in 1:nOrg) { if (i == 1) { Sampling.Rate.Mat[,i] = c(sample(seq(1.5,2,0.05), size = nSamples, replace = TRUE), sample(seq(0.5,0.667,0.05), size = nSamples, replace = TRUE)) } else if (i == 2) { Sampling.Rate.Mat[,i] = c(sample(seq(0.5,0.667,0.05), size = nSamples, replace = TRUE), sample(seq(1.5,2,0.05), size = nSamples, replace = TRUE)) } else { Sampling.Rate.Mat[,i] = sample(seq(0.5,2,0.1), size = nSamples*2, replace = TRUE) } } } resA = generate.org.mat(nSamples = nSamples, nOrg = nOrg, CountVec = CountVec, Sampling.Rate.Mat = Sampling.Rate.Mat) for (i in 1:nOrg) { UP.reg = c(UP.reg,resA[[i]]@variable.annotations$upregulation) DOWN.reg = c(DOWN.reg,resA[[i]]@variable.annotations$downregulation) if (i == 1) { G.Mat <- resA[[i]]@count.matrix T.Mat.norm <- DESeq2.norm.mat(resA[[i]]@count.matrix, cond.vec,type.vec) } else { G.Mat = rbind(G.Mat,resA[[i]]@count.matrix) T.Mat.norm = rbind(T.Mat.norm,DESeq2.norm.mat(resA[[i]]@count.matrix, cond.vec,type.vec)) } } rownames(G.Mat) <- as.character(1:dim(G.Mat)[1]) rownames(T.Mat.norm) <- as.character(1:dim(T.Mat.norm)[1]) G.Mat.norm = DESeq2.norm.mat(G.Mat, cond.vec,type.vec) res.tax <- DESeq2.result(Xmat = T.Mat.norm, cond = cond.vec, type = type.vec) res.glo <- DESeq2.result(Xmat = G.Mat.norm, cond = cond.vec, type = type.vec) res.tax$padj[is.na(res.tax$padj)] = 1 res.glo$padj[is.na(res.glo$padj)] = 1 T.Mat = de.counts.with.direction.per.org(Compcoder.List = resA, DESeq.result = res.tax, p.thresh = p.thresh) G.Mat = de.counts.with.direction.per.org(Compcoder.List = resA, DESeq.result = res.glo, p.thresh = p.thresh) return(list(T.Mat = T.Mat, G.Mat = G.Mat, res.tax.padj = res.tax$padj, res.glo.padj = res.glo$padj, res.tax.log2FoldChange = res.tax$log2FoldChange, res.glo.log2FoldChange = res.glo$log2FoldChange, UP.reg = UP.reg, DOWN.reg = DOWN.reg)) } run.sim.tax.glo4 <- function(experiment, p.thresh) { require(compcodeR) require(DESeq2) nOrg = 5 #Simulation: number of organism to simulate nSamples = 6 #Simulation: number of samples per condition CountVec = c(1e7,5e6,1e6,5e5,1e5) #Simulation: base library size Sampling.Rate.Mat = matrix(1,ncol = nOrg, nrow = nSamples*2) #Simulation: sampling rates for each organism and sample cond.vec <- c(rep('A',nSamples),rep('B',nSamples)) #conditions for each sample type.vec <- rep('1',2*nSamples) #type of read if (experiment == 'sim4') { nOrg = 2 #Simulation: number of organism to simulate CountVec = c(1e7,1e7) #Simulation: base library size Sampling.Rate.Mat = matrix(1,ncol = nOrg, nrow = nSamples*2) #Simulation: sampling rates for each organism and sample for (i in 1:nOrg) { if (i == 1) { Sampling.Rate.Mat[,i] = c(sample(seq(1.5,2,0.05), size = nSamples, replace = TRUE), sample(seq(0.5,0.667,0.05), size = nSamples, replace = TRUE)) } else if (i == 2) { Sampling.Rate.Mat[,i] = c(sample(seq(0.5,0.667,0.05), size = nSamples, replace = TRUE), sample(seq(1.5,2,0.05), size = nSamples, replace = TRUE)) } } } resA = generate.org.mat(nSamples = nSamples, nOrg = nOrg, CountVec = CountVec, Sampling.Rate.Mat = Sampling.Rate.Mat) UP.reg = c() DOWN.reg = c() for (i in 1:nOrg) { UP.reg = c(UP.reg,resA[[i]]@variable.annotations$upregulation) DOWN.reg = c(DOWN.reg,resA[[i]]@variable.annotations$downregulation) if (i == 1) { G.Mat <- resA[[i]]@count.matrix T.Mat.norm <- DESeq2.norm.mat(resA[[i]]@count.matrix, cond.vec,type.vec) } else { G.Mat = G.Mat + resA[[i]]@count.matrix val = DESeq2.norm.mat(resA[[i]]@count.matrix, cond.vec,type.vec) T.Mat.norm2 = T.Mat.norm + val T.Mat.norm = rbind(T.Mat.norm,val) } } rownames(G.Mat) <- as.character(1:dim(G.Mat)[1]) rownames(T.Mat.norm) <- as.character(1:dim(T.Mat.norm)[1]) G.Mat.norm = DESeq2.norm.mat(G.Mat, cond.vec,type.vec) res.tax <- DESeq2.result(Xmat = T.Mat.norm, cond = cond.vec, type = type.vec) res.tax2 <- DESeq2.result(Xmat = T.Mat.norm2, cond = cond.vec, type = type.vec) res.glo <- DESeq2.result(Xmat = G.Mat.norm, cond = cond.vec, type = type.vec) res.tax$padj[is.na(res.tax$padj)] = 1 res.glo$padj[is.na(res.glo$padj)] = 1 X.stacked = log(rowMeans(T.Mat.norm[,1:6])) Y.stacked = log(rowMeans(T.Mat.norm[,7:12])) X.super.tax = log(rowMeans(T.Mat.norm2[,1:6])) Y.super.tax = log(rowMeans(T.Mat.norm2[,7:12])) X.super.glo = log(rowMeans(G.Mat.norm[,1:6])) Y.super.glo = log(rowMeans(G.Mat.norm[,7:12])) T.Mat = de.counts.with.direction.per.org(Compcoder.List = resA, DESeq.result = res.tax, p.thresh = p.thresh) G.Mat = de.counts.with.direction.per.org(Compcoder.List = resA, DESeq.result = res.glo, p.thresh = p.thresh) return(list(T.Mat = T.Mat, G.Mat = G.Mat, X.stacked = X.stacked, Y.stacked = Y.stacked, X.super.tax = X.super.tax, Y.super.tax = Y.super.tax, X.super.glo = X.super.glo, Y.super.glo = Y.super.glo, UP.reg = UP.reg, DOWN.reg = DOWN.reg)) } real.data <- function(file, condA = '13days', condB = '27days', pval = 0.05) { X = read.csv(file) if (sum(dim(X) == c(210210,5)) != 2) { cat('incorrect csv provided.','\n') return(-1) } Start.Vec = seq(1,210210,3185) End.Vec = seq(3185,210210,3185) cond.names <- c("13days", "15days", "16days", "27days", "29days", "30days") spec.names <- c("BACCAC", "BACOVA", "BACUNI", "BDI", "BT", "BVU", "CLOSCI", "CLOSPI", "COLAER", "DORLON", "RUMOBE") type.vec = rep('1',8) #1 org, 6 cond k = 1 List <- list() for (i in 1:11) { List[[i]] <- list() for (j in 1:6) { List[[i]][[j]] <- X[Start.Vec[k]:End.Vec[k],2:5] rownames(List[[i]][[j]]) <- Pfam.annot() k = k + 1 } names(List[[i]]) <- cond.names } names(List) <- spec.names Result.Transcriptome.List <- list() Result.Metatranscriptome.List <- list() Normalized.Transcriptome <-list() Global.Mat <- matrix(0,3185,8) Tax.Mat.norm <- Global.Mat for (i in 1:11) { cond.vec <- c(rep(condA,4),rep(condB,4)) XMat <- cbind(as.matrix(List[[i]][[condA]]),as.matrix(List[[i]][[condB]])) ZMat = XMat != 0 if (sum(rowSums(ZMat) == 8) > 0) { YMat <- DESeq2.norm.mat(XMat, cond.vec,type.vec) Global.Mat = Global.Mat + XMat Tax.Mat.norm = Tax.Mat.norm + YMat Normalized.Transcriptome[[spec.names[i]]] <- YMat res <- DESeq2.result(Xmat = YMat, cond = cond.vec, type = type.vec) Result.Transcriptome.List[[spec.names[i]]] <- res } } Global.Mat.norm <- DESeq2.norm.mat(Global.Mat, cond.vec,type.vec) res.glo <- DESeq2.result(Xmat = Global.Mat.norm, cond = cond.vec, type = type.vec) res.tax <- DESeq2.result(Xmat = Tax.Mat.norm, cond = cond.vec, type = type.vec) res.glo$padj[is.na(res.glo$padj)] = 1 res.tax$padj[is.na(res.tax$padj)] = 1 res.glo$log2FoldChange[is.na(res.glo$log2FoldChange)] = 0 res.tax$log2FoldChange[is.na(res.tax$log2FoldChange)] = 0 T.Mat = sapply(names(Result.Transcriptome.List), function(x) Result.Transcriptome.List[[x]]$padj < pval) T.Mat[is.na(T.Mat)] = FALSE T.Vec = rowSums(T.Mat) shared.up = sum(res.glo$padj < pval & res.tax$padj < pval & res.glo$log2FoldChange > 0 & res.tax$log2FoldChange > 0) shared.down = sum(res.glo$padj < pval & res.tax$padj < pval & res.glo$log2FoldChange < 0 & res.tax$log2FoldChange < 0) diff.dir.up = sum(res.glo$padj < pval & res.tax$padj < pval & res.glo$log2FoldChange > 0 & res.tax$log2FoldChange < 0) diff.dir.down = sum(res.glo$padj < pval & res.tax$padj < pval & res.glo$log2FoldChange < 0 & res.tax$log2FoldChange > 0) extra.tax = sum(res.glo$padj > pval & res.tax$padj < pval) extra.glo = sum(res.glo$padj < pval & res.tax$padj > pval) shared.sum = shared.up + shared.down diff.dir.sum = diff.dir.up + diff.dir.down Tax.Vec = c(extra.tax,diff.dir.sum,shared.sum) Glo.Vec = c(extra.glo,diff.dir.sum,shared.sum) x11() par(mfrow = c(1, 3)) plot(table(T.Vec), main = 'transcriptome', xlab = 'shared DEF between transcriptomes', ylab = 'count') plot(table(T.Vec[res.glo$padj < pval]), main = 'glo', xlab = 'DEF shared with single transcriptomes', ylab = 'count') plot(table(T.Vec[res.tax$padj < pval]), main = 'tax', xlab = 'DEF shared with single transcriptomes', ylab = 'count') x11() df <- data.frame(x = c(Tax.Vec,Glo.Vec), y = rep(c('tax','glo'), each = 3), col = rep(factor(c('unique','opposite direction','shared'), levels = c('unique','opposite direction','shared')),2)) p <- ggplot(data = df, aes(x = y, y = x)) p <- p + geom_bar(stat = 'identity', aes(fill = col)) + ggtitle(paste0('DEF for ',condA,' vs ',condB)) + xlab('method') + ylab('count') print(p) return(p + geom_abline(intercept = 0, slope = 1, alpha = 0.25, linetype = 'dashed') + facet_wrap(~type) + theme(legend.position="bottom")) } Pfam.annot <- function() { X = c(4 ,5 ,6 ,11 ,12 ,13 ,23 ,27 ,34 ,35 ,37 ,41 ,44 ,56 ,69 ,72 ,75 ,76 ,78 ,79 ,80 ,81 ,82 ,83 ,85 ,89 ,91 ,92 ,106 , 107 ,109 ,111 ,112 ,113 ,117 ,118 ,119 ,120 ,121 ,122 ,126 ,128 ,132 ,133 ,135 ,137 ,140 ,144 ,145 ,146 ,148 ,149 ,150 ,152 ,154 ,155 ,156 ,158 , 160 ,162 ,163 ,164 ,165 ,166 ,171 ,175 ,176 ,177 ,180 ,181 ,183 ,185 ,186 ,188 ,190 ,196 ,198 ,199 ,202 ,203 ,204 ,205 ,206 ,207 ,208 ,209 ,210 , 213 ,215 ,216 ,218 ,221 ,224 ,226 ,230 ,231 ,232 ,237 ,238 ,239 ,245 ,246 ,248 ,251 ,252 ,253 ,254 ,255 ,258 ,265 ,266 ,268 ,269 ,270 ,271 ,275 , 276 ,278 ,281 ,282 ,285 ,288 ,289 ,290 ,291 ,293 ,294 ,295 ,297 ,298 ,300 ,301 ,302 ,303 ,306 ,308 ,309 ,310 ,311 ,312 ,313 ,317 ,318 ,326 ,327 , 328 ,329 ,330 ,331 ,333 ,334 ,338 ,342 ,343 ,344 ,346 ,347 ,348 ,349 ,350 ,355 ,356 ,358 ,359 ,361 ,364 ,365 ,366 ,367 ,370 ,374 ,375 ,376 ,378 , 380 ,381 ,383 ,384 ,388 ,389 ,390 ,391 ,392 ,393 ,398 ,401 ,403 ,408 ,410 ,411 ,413 ,416 ,420 ,425 ,430 ,436 ,437 ,438 ,440 ,441 ,444 ,448 ,453 , 455 ,456 ,459 ,462 ,464 ,465 ,466 ,467 ,468 ,471 ,472 ,474 ,475 ,476 ,478 ,479 ,480 ,481 ,482 ,483 ,484 ,485 ,486 ,488 ,490 ,491 ,496 ,497 ,498 , 499 ,501 ,507 ,512 ,515 ,520 ,521 ,528 ,529 ,532 ,533 ,534 ,535 ,542 ,543 ,545 ,549 ,550 ,551 ,557 ,561 ,562 ,563 ,565 ,570 ,571 ,572 ,573 ,574 , 575 ,578 ,579 ,580 ,581 ,582 ,583 ,584 ,586 ,587 ,588 ,589 ,590 ,591 ,593 ,595 ,596 ,623 ,625 ,633 ,636 ,639 ,648 ,654 ,657 ,662 ,664 ,665 ,668 , 670 ,672 ,673 ,675 ,676 ,677 ,679 ,682 ,684 ,686 ,687 ,689 ,690 ,691 ,692 ,694 ,696 ,697 ,698 ,701 ,703 ,704 ,707 ,708 ,709 ,710 ,712 ,717 ,722 , 723 ,724 ,728 ,730 ,731 ,733 ,742 ,745 ,746 ,749 ,750 ,753 ,754 ,756 ,759 ,763 ,764 ,766 ,768 ,772 ,781 ,793 ,795 ,800 ,801 ,805 ,809 ,814 ,815 , 817 ,825 ,828 ,829 ,830 ,831 ,834 ,842 ,847 ,849 ,852 ,854 ,857 ,860 ,861 ,871 ,872 ,873 ,874 ,877 ,881 ,882 ,884 ,885 ,886 ,889 ,890 ,891 ,892 , 893 ,899 ,902 ,903 ,905 ,908 ,909 ,912 ,916 ,919 ,920 ,923 ,924 ,925 ,926 ,929 ,930 ,932 ,933 ,936 ,939 ,941 ,950 ,953 ,958 ,959 ,962 ,970 ,977 , 982 ,984 ,986 ,988 ,989 ,990 ,994 ,999 ,1000 ,1008 ,1011 ,1012 ,1016 ,1018 ,1022 ,1025 ,1026 ,1029 ,1032 ,1035 ,1037 ,1039 ,1040 ,1041 ,1042 ,1043 ,1047 ,1048 ,1050 , 1053 ,1055 ,1058 ,1061 ,1063 ,1066 ,1070 ,1071 ,1074 ,1075 ,1076 ,1077 ,1078 ,1081 ,1084 ,1087 ,1095 ,1098 ,1103 ,1106 ,1112 ,1113 ,1116 ,1118 ,1119 ,1120 ,1121 ,1128 ,1131 , 1132 ,1134 ,1136 ,1138 ,1139 ,1144 ,1145 ,1149 ,1155 ,1156 ,1165 ,1168 ,1170 ,1171 ,1174 ,1175 ,1176 ,1177 ,1180 ,1182 ,1183 ,1187 ,1189 ,1192 ,1193 ,1195 ,1196 ,1197 ,1202 , 1204 ,1205 ,1206 ,1207 ,1208 ,1209 ,1210 ,1212 ,1219 ,1220 ,1223 ,1225 ,1226 ,1227 ,1228 ,1230 ,1232 ,1235 ,1238 ,1242 ,1243 ,1244 ,1245 ,1248 ,1250 ,1252 ,1253 ,1255 ,1256 , 1258 ,1259 ,1261 ,1262 ,1263 ,1264 ,1266 ,1268 ,1272 ,1276 ,1281 ,1288 ,1289 ,1292 ,1293 ,1297 ,1300 ,1301 ,1314 ,1315 ,1316 ,1321 ,1325 ,1326 ,1327 ,1330 ,1336 ,1337 ,1343 , 1344 ,1345 ,1346 ,1348 ,1351 ,1364 ,1367 ,1368 ,1370 ,1371 ,1379 ,1380 ,1381 ,1384 ,1385 ,1386 ,1396 ,1402 ,1406 ,1408 ,1409 ,1411 ,1416 ,1418 ,1420 ,1424 ,1425 ,1427 ,1430 , 1431 ,1432 ,1433 ,1434 ,1435 ,1436 ,1443 ,1450 ,1451 ,1455 ,1458 ,1464 ,1467 ,1470 ,1471 ,1472 ,1473 ,1475 ,1476 ,1478 ,1479 ,1487 ,1488 ,1493 ,1496 ,1497 ,1501 ,1502 ,1503 , 1507 ,1509 ,1510 ,1512 ,1513 ,1515 ,1520 ,1522 ,1523 ,1527 ,1531 ,1541 ,1544 ,1545 ,1546 ,1547 ,1548 ,1551 ,1553 ,1554 ,1555 ,1556 ,1557 ,1558 ,1565 ,1566 ,1568 ,1569 ,1571 , 1575 ,1578 ,1580 ,1588 ,1590 ,1592 ,1594 ,1595 ,1596 ,1597 ,1609 ,1610 ,1612 ,1613 ,1614 ,1618 ,1619 ,1624 ,1625 ,1627 ,1628 ,1632 ,1634 ,1636 ,1637 ,1638 ,1640 ,1641 ,1642 , 1643 ,1645 ,1648 ,1649 ,1653 ,1654 ,1656 ,1658 ,1661 ,1663 ,1668 ,1676 ,1678 ,1680 ,1687 ,1693 ,1694 ,1695 ,1699 ,1702 ,1709 ,1713 ,1715 ,1717 ,1719 ,1725 ,1726 ,1728 ,1729 , 1734 ,1738 ,1740 ,1741 ,1743 ,1746 ,1750 ,1751 ,1755 ,1757 ,1758 ,1761 ,1765 ,1769 ,1773 ,1782 ,1783 ,1784 ,1790 ,1791 ,1795 ,1797 ,1799 ,1807 ,1808 ,1809 ,1810 ,1812 ,1813 , 1814 ,1817 ,1820 ,1823 ,1832 ,1833 ,1835 ,1841 ,1842 ,1844 ,1850 ,1855 ,1862 ,1863 ,1865 ,1867 ,1869 ,1872 ,1878 ,1880 ,1882 ,1883 ,1888 ,1890 ,1891 ,1894 ,1895 ,1902 ,1906 , 1909 ,1910 ,1914 ,1915 ,1923 ,1924 ,1925 ,1926 ,1928 ,1930 ,1934 ,1935 ,1936 ,1937 ,1938 ,1943 ,1944 ,1948 ,1955 ,1957 ,1958 ,1960 ,1963 ,1964 ,1965 ,1966 ,1967 ,1969 ,1975 , 1978 ,1979 ,1980 ,1985 ,1987 ,1988 ,1990 ,1991 ,1992 ,1996 ,2001 ,2012 ,2016 ,2018 ,2021 ,2026 ,2028 ,2033 ,2037 ,2055 ,2056 ,2065 ,2073 ,2074 ,2075 ,2080 ,2082 ,2086 ,2091 , 2092 ,2096 ,2110 ,2113 ,2127 ,2129 ,2130 ,2132 ,2133 ,2142 ,2146 ,2151 ,2152 ,2153 ,2156 ,2163 ,2169 ,2178 ,2190 ,2195 ,2219 ,2222 ,2224 ,2225 ,2229 ,2230 ,2237 ,2253 ,2254 , 2255 ,2256 ,2261 ,2272 ,2275 ,2277 ,2278 ,2283 ,2286 ,2287 ,2288 ,2302 ,2308 ,2310 ,2311 ,2317 ,2321 ,2322 ,2325 ,2329 ,2335 ,2347 ,2348 ,2350 ,2353 ,2355 ,2357 ,2358 ,2361 , 2367 ,2368 ,2371 ,2374 ,2378 ,2381 ,2384 ,2386 ,2388 ,2390 ,2397 ,2401 ,2403 ,2405 ,2416 ,2417 ,2421 ,2423 ,2424 ,2435 ,2436 ,2441 ,2445 ,2446 ,2447 ,2449 ,2452 ,2457 ,2463 , 2464 ,2469 ,2470 ,2472 ,2481 ,2482 ,2485 ,2486 ,2491 ,2492 ,2493 ,2494 ,2498 ,2502 ,2503 ,2504 ,2508 ,2511 ,2514 ,2515 ,2517 ,2518 ,2525 ,2527 ,2530 ,2534 ,2535 ,2537 ,2540 , 2541 ,2542 ,2544 ,2545 ,2547 ,2548 ,2550 ,2553 ,2554 ,2556 ,2557 ,2558 ,2559 ,2562 ,2563 ,2565 ,2566 ,2568 ,2569 ,2570 ,2571 ,2572 ,2574 ,2575 ,2576 ,2577 ,2578 ,2579 ,2580 , 2581 ,2583 ,2585 ,2586 ,2588 ,2589 ,2590 ,2591 ,2592 ,2595 ,2596 ,2597 ,2601 ,2602 ,2603 ,2604 ,2606 ,2607 ,2608 ,2609 ,2610 ,2614 ,2615 ,2616 ,2617 ,2618 ,2620 ,2622 ,2625 , 2627 ,2629 ,2630 ,2631 ,2632 ,2633 ,2635 ,2637 ,2638 ,2639 ,2643 ,2645 ,2646 ,2650 ,2652 ,2653 ,2654 ,2655 ,2657 ,2659 ,2660 ,2661 ,2662 ,2664 ,2666 ,2669 ,2670 ,2673 ,2674 , 2675 ,2677 ,2678 ,2683 ,2684 ,2686 ,2687 ,2689 ,2690 ,2698 ,2699 ,2702 ,2705 ,2706 ,2709 ,2719 ,2729 ,2730 ,2733 ,2734 ,2738 ,2739 ,2742 ,2743 ,2744 ,2746 ,2748 ,2749 ,2754 , 2767 ,2768 ,2769 ,2770 ,2771 ,2772 ,2773 ,2774 ,2775 ,2776 ,2777 ,2779 ,2780 ,2781 ,2782 ,2784 ,2785 ,2786 ,2787 ,2796 ,2800 ,2801 ,2806 ,2810 ,2811 ,2812 ,2817 ,2823 ,2824 , 2826 ,2829 ,2833 ,2836 ,2837 ,2838 ,2843 ,2844 ,2852 ,2861 ,2863 ,2866 ,2867 ,2870 ,2872 ,2873 ,2874 ,2875 ,2878 ,2879 ,2880 ,2881 ,2882 ,2884 ,2885 ,2887 ,2894 ,2896 ,2897 , 2899 ,2901 ,2903 ,2906 ,2910 ,2911 ,2912 ,2913 ,2915 ,2922 ,2926 ,2927 ,2929 ,2934 ,2938 ,2952 ,2954 ,2963 ,2965 ,2971 ,2976 ,2978 ,3006 ,3008 ,3009 ,3023 ,3030 ,3031 ,3033 , 3050 ,3051 ,3060 ,3061 ,3062 ,3063 ,3065 ,3073 ,3092 ,3099 ,3102 ,3109 ,3116 ,3118 ,3119 ,3120 ,3123 ,3129 ,3143 ,3144 ,3147 ,3150 ,3160 ,3167 ,3169 ,3173 ,3176 ,3180 ,3186 , 3193 ,3200 ,3205 ,3235 ,3237 ,3253 ,3255 ,3264 ,3266 ,3275 ,3279 ,3288 ,3309 ,3313 ,3315 ,3319 ,3321 ,3323 ,3328 ,3330 ,3331 ,3349 ,3352 ,3354 ,3358 ,3372 ,3374 ,3389 ,3400 , 3412 ,3413 ,3414 ,3415 ,3417 ,3418 ,3419 ,3422 ,3432 ,3435 ,3437 ,3442 ,3444 ,3446 ,3447 ,3448 ,3449 ,3450 ,3453 ,3454 ,3457 ,3458 ,3459 ,3460 ,3461 ,3462 ,3466 ,3471 ,3473 , 3477 ,3479 ,3480 ,3481 ,3483 ,3484 ,3485 ,3534 ,3544 ,3547 ,3548 ,3551 ,3553 ,3572 ,3575 ,3576 ,3577 ,3588 ,3590 ,3591 ,3592 ,3595 ,3598 ,3599 ,3600 ,3601 ,3605 ,3606 ,3609 , 3610 ,3611 ,3613 ,3616 ,3618 ,3625 ,3629 ,3630 ,3631 ,3632 ,3633 ,3636 ,3641 ,3648 ,3649 ,3652 ,3663 ,3668 ,3672 ,3690 ,3692 ,3702 ,3703 ,3704 ,3706 ,3711 ,3717 ,3719 ,3720 , 3721 ,3724 ,3725 ,3726 ,3727 ,3729 ,3733 ,3734 ,3738 ,3739 ,3740 ,3747 ,3749 ,3755 ,3764 ,3772 ,3773 ,3775 ,3776 ,3780 ,3781 ,3786 ,3787 ,3788 ,3793 ,3796 ,3799 ,3802 ,3806 , 3807 ,3808 ,3812 ,3814 ,3816 ,3819 ,3830 ,3837 ,3838 ,3840 ,3841 ,3845 ,3852 ,3853 ,3861 ,3862 ,3869 ,3880 ,3883 ,3894 ,3900 ,3928 ,3929 ,3932 ,3938 ,3946 ,3947 ,3948 ,3949 , 3950 ,3951 ,3952 ,3956 ,3960 ,3965 ,3976 ,3977 ,3979 ,3989 ,3992 ,3993 ,4002 ,4011 ,4014 ,4015 ,4018 ,4020 ,4023 ,4024 ,4025 ,4026 ,4028 ,4041 ,4055 ,4060 ,4069 ,4070 ,4071 , 4072 ,4073 ,4074 ,4079 ,4085 ,4087 ,4093 ,4095 ,4101 ,4127 ,4131 ,4138 ,4143 ,4166 ,4167 ,4170 ,4172 ,4173 ,4187 ,4198 ,4199 ,4203 ,4204 ,4205 ,4221 ,4224 ,4226 ,4227 ,4229 , 4230 ,4232 ,4233 ,4235 ,4237 ,4239 ,4245 ,4246 ,4253 ,4255 ,4261 ,4263 ,4265 ,4266 ,4277 ,4290 ,4294 ,4295 ,4296 ,4297 ,4298 ,4301 ,4304 ,4306 ,4309 ,4313 ,4321 ,4324 ,4326 , 4327 ,4336 ,4342 ,4343 ,4346 ,4355 ,4357 ,4371 ,4378 ,4383 ,4389 ,4390 ,4392 ,4397 ,4402 ,4404 ,4405 ,4413 ,4422 ,4432 ,4434 ,4439 ,4445 ,4448 ,4450 ,4452 ,4459 ,4461 ,4463 , 4464 ,4466 ,4468 ,4471 ,4472 ,4474 ,4480 ,4488 ,4519 ,4531 ,4536 ,4539 ,4542 ,4545 ,4551 ,4552 ,4556 ,4560 ,4561 ,4563 ,4565 ,4577 ,4586 ,4587 ,4607 ,4608 ,4613 ,4616 ,4630 , 4647 ,4657 ,4685 ,4715 ,4723 ,4754 ,4760 ,4773 ,4794 ,4816 ,4851 ,4854 ,4860 ,4865 ,4879 ,4883 ,4892 ,4898 ,4909 ,4916 ,4932 ,4952 ,4960 ,4961 ,4962 ,4963 ,4965 ,4970 ,4973 , 4977 ,4983 ,4984 ,4986 ,4991 ,4993 ,4994 ,4997 ,4998 ,5000 ,5015 ,5016 ,5017 ,5025 ,5031 ,5036 ,5043 ,5050 ,5063 ,5065 ,5066 ,5067 ,5089 ,5096 ,5099 ,5103 ,5105 ,5107 ,5117 , 5119 ,5124 ,5133 ,5135 ,5139 ,5140 ,5147 ,5164 ,5167 ,5168 ,5170 ,5173 ,5175 ,5188 ,5190 ,5191 ,5192 ,5193 ,5194 ,5195 ,5198 ,5201 ,5221 ,5222 ,5223 ,5239 ,5256 ,5257 ,5258 , 5272 ,5336 ,5343 ,5345 ,5362 ,5389 ,5402 ,5406 ,5426 ,5430 ,5437 ,5448 ,5488 ,5491 ,5496 ,5521 ,5523 ,5524 ,5525 ,5534 ,5544 ,5552 ,5569 ,5572 ,5576 ,5580 ,5582 ,5592 ,5593 , 5598 ,5635 ,5636 ,5649 ,5656 ,5670 ,5673 ,5681 ,5683 ,5684 ,5685 ,5686 ,5691 ,5693 ,5697 ,5698 ,5704 ,5708 ,5709 ,5711 ,5713 ,5717 ,5731 ,5738 ,5746 ,5762 ,5766 ,5816 ,5838 , 5857 ,5866 ,5869 ,5895 ,5896 ,5913 ,5922 ,5935 ,5958 ,5960 ,5970 ,5971 ,5973 ,5977 ,5979 ,5991 ,6018 ,6026 ,6032 ,6050 ,6054 ,6056 ,6071 ,6074 ,6082 ,6094 ,6100 ,6107 ,6114 , 6115 ,6123 ,6125 ,6130 ,6133 ,6134 ,6135 ,6144 ,6149 ,6152 ,6153 ,6160 ,6161 ,6165 ,6166 ,6177 ,6180 ,6182 ,6199 ,6202 ,6207 ,6210 ,6224 ,6245 ,6250 ,6257 ,6271 ,6280 ,6283 , 6293 ,6296 ,6321 ,6353 ,6379 ,6397 ,6406 ,6414 ,6415 ,6418 ,6421 ,6439 ,6445 ,6452 ,6458 ,6463 ,6470 ,6480 ,6496 ,6508 ,6537 ,6541 ,6560 ,6574 ,6580 ,6603 ,6605 ,6612 ,6628 , 6675 ,6686 ,6689 ,6695 ,6713 ,6723 ,6725 ,6733 ,6738 ,6739 ,6750 ,6772 ,6782 ,6792 ,6803 ,6808 ,6810 ,6824 ,6826 ,6827 ,6831 ,6832 ,6835 ,6838 ,6841 ,6854 ,6857 ,6866 ,6874 , 6889 ,6898 ,6902 ,6908 ,6912 ,6925 ,6930 ,6935 ,6940 ,6941 ,6949 ,6953 ,6961 ,6962 ,6964 ,6965 ,6968 ,6969 ,6970 ,6971 ,7005 ,7009 ,7022 ,7030 ,7075 ,7081 ,7083 ,7085 ,7090 , 7136 ,7155 ,7168 ,7171 ,7179 ,7228 ,7242 ,7244 ,7247 ,7261 ,7275 ,7282 ,7295 ,7302 ,7313 ,7314 ,7332 ,7352 ,7355 ,7364 ,7366 ,7371 ,7377 ,7385 ,7396 ,7432 ,7441 ,7450 ,7451 , 7454 ,7456 ,7463 ,7470 ,7475 ,7477 ,7478 ,7479 ,7486 ,7488 ,7494 ,7495 ,7497 ,7498 ,7499 ,7501 ,7503 ,7505 ,7508 ,7510 ,7516 ,7517 ,7520 ,7521 ,7523 ,7532 ,7537 ,7538 ,7549 , 7554 ,7555 ,7556 ,7561 ,7581 ,7582 ,7584 ,7602 ,7603 ,7610 ,7632 ,7638 ,7650 ,7653 ,7659 ,7660 ,7661 ,7662 ,7664 ,7669 ,7670 ,7676 ,7680 ,7681 ,7683 ,7685 ,7687 ,7690 ,7691 , 7693 ,7697 ,7698 ,7702 ,7703 ,7715 ,7719 ,7722 ,7724 ,7726 ,7728 ,7729 ,7730 ,7733 ,7739 ,7745 ,7748 ,7751 ,7784 ,7786 ,7804 ,7805 ,7811 ,7826 ,7831 ,7837 ,7853 ,7859 ,7863 , 7873 ,7875 ,7876 ,7881 ,7882 ,7883 ,7884 ,7885 ,7892 ,7905 ,7907 ,7927 ,7929 ,7934 ,7940 ,7943 ,7944 ,7947 ,7949 ,7959 ,7963 ,7969 ,7971 ,7973 ,7977 ,7980 ,7991 ,7992 ,7993 , 7994 ,8000 ,8002 ,8006 ,8011 ,8013 ,8019 ,8020 ,8029 ,8032 ,8124 ,8125 ,8206 ,8207 ,8212 ,8218 ,8220 ,8222 ,8230 ,8238 ,8239 ,8240 ,8241 ,8244 ,8245 ,8264 ,8274 ,8275 ,8278 , 8279 ,8281 ,8282 ,8291 ,8299 ,8305 ,8306 ,8307 ,8309 ,8323 ,8328 ,8331 ,8338 ,8340 ,8348 ,8352 ,8353 ,8357 ,8359 ,8367 ,8378 ,8388 ,8401 ,8402 ,8436 ,8439 ,8442 ,8443 ,8447 , 8448 ,8450 ,8455 ,8459 ,8463 ,8478 ,8481 ,8483 ,8485 ,8486 ,8497 ,8501 ,8502 ,8503 ,8522 ,8529 ,8530 ,8531 ,8532 ,8533 ,8534 ,8541 ,8543 ,8544 ,8545 ,8546 ,8645 ,8660 ,8664 , 8665 ,8666 ,8667 ,8669 ,8676 ,8706 ,8707 ,8708 ,8713 ,8747 ,8753 ,8757 ,8759 ,8760 ,8761 ,8769 ,8774 ,8780 ,8797 ,8800 ,8808 ,8821 ,8837 ,8840 ,8841 ,8842 ,8843 ,8849 ,8867 , 8876 ,8878 ,8885 ,8889 ,8890 ,8901 ,8902 ,8903 ,8922 ,8924 ,8937 ,8942 ,8973 ,8980 ,8984 ,8989 ,9035 ,9084 ,9087 ,9092 ,9093 ,9106 ,9107 ,9112 ,9113 ,9115 ,9148 ,9180 ,9186 , 9190 ,9250 ,9261 ,9269 ,9278 ,9285 ,9295 ,9296 ,9297 ,9299 ,9314 ,9334 ,9335 ,9338 ,9339 ,9357 ,9359 ,9363 ,9364 ,9365 ,9369 ,9371 ,9374 ,9382 ,9383 ,9391 ,9393 ,9397 ,9411 , 9413 ,9414 ,9424 ,9458 ,9471 ,9479 ,9491 ,9492 ,9508 ,9512 ,9515 ,9524 ,9527 ,9546 ,9547 ,9548 ,9551 ,9560 ,9561 ,9578 ,9579 ,9581 ,9586 ,9587 ,9588 ,9594 ,9603 ,9604 ,9605 , 9633 ,9640 ,9643 ,9660 ,9661 ,9664 ,9669 ,9671 ,9674 ,9681 ,9693 ,9704 ,9709 ,9711 ,9719 ,9723 ,9818 ,9820 ,9822 ,9825 ,9826 ,9827 ,9832 ,9835 ,9837 ,9848 ,9855 ,9858 ,9861 , 9862 ,9889 ,9903 ,9907 ,9912 ,9913 ,9918 ,9919 ,9922 ,9924 ,9925 ,9951 ,9952 ,9960 ,9966 ,9967 ,9972 ,9976 ,9983 ,9989 ,9991 ,9992 ,10000 ,10021 ,10022 ,10031 ,10035 ,10040 ,10050 , 10054 ,10060 ,10066 ,10076 ,10077 ,10087 ,10091 ,10101 ,10105 ,10112 ,10114 ,10117 ,10133 ,10139 ,10140 ,10143 ,10145 ,10150 ,10250 ,10269 ,10282 ,10298 ,10369 ,10371 ,10385 ,10387 ,10396 ,10397 ,10410 , 10412 ,10415 ,10417 ,10418 ,10431 ,10437 ,10438 ,10458 ,10459 ,10502 ,10509 ,10518 ,10531 ,10543 ,10544 ,10555 ,10566 ,10571 ,10576 ,10588 ,10589 ,10590 ,10592 ,10604 ,10609 ,10626 ,10628 ,10632 ,10633 , 10635 ,10646 ,10651 ,10662 ,10672 ,10677 ,10678 ,10694 ,10704 ,10712 ,10722 ,10728 ,10747 ,10771 ,10825 ,10882 ,10899 ,10902 ,10923 ,10926 ,10934 ,10951 ,10957 ,10962 ,10974 ,10979 ,10988 ,10989 ,10990 , 10991 ,10996 ,11007 ,11013 ,11026 ,11028 ,11042 ,11066 ,11074 ,11133 ,11155 ,11167 ,11175 ,11185 ,11187 ,11188 ,11193 ,11195 ,11208 ,11236 ,11258 ,11276 ,11297 ,11306 ,11307 ,11311 ,11325 ,11335 ,11368 , 11369 ,11380 ,11387 ,11391 ,11392 ,11396 ,11412 ,11436 ,11449 ,11490 ,11495 ,11545 ,11551 ,11589 ,11638 ,11644 ,11651 ,11655 ,11667 ,11672 ,11680 ,11700 ,11721 ,11734 ,11738 ,11750 ,11751 ,11756 ,11759 , 11760 ,11761 ,11762 ,11775 ,11777 ,11784 ,11790 ,11795 ,11796 ,11798 ,11799 ,11810 ,11823 ,11832 ,11838 ,11839 ,11842 ,11849 ,11855 ,11867 ,11870 ,11888 ,11890 ,11893 ,11897 ,11906 ,11907 ,11922 ,11941 , 11950 ,11958 ,11959 ,11967 ,11973 ,11974 ,11975 ,11987 ,11997 ,12002 ,12008 ,12029 ,12034 ,12051 ,12072 ,12080 ,12081 ,12083 ,12099 ,12116 ,12127 ,12128 ,12146 ,12161 ,12164 ,12169 ,12186 ,12206 ,12215 , 12224 ,12225 ,12229 ,12245 ,12320 ,12327 ,12344 ,12345 ,12385 ,12388 ,12389 ,12390 ,12392 ,12399 ,12401 ,12408 ,12412 ,12437 ,12439 ,12450 ,12464 ,12476 ,12508 ,12541 ,12558 ,12627 ,12631 ,12636 ,12637 , 12641 ,12642 ,12643 ,12645 ,12646 ,12648 ,12650 ,12652 ,12654 ,12663 ,12664 ,12666 ,12667 ,12668 ,12670 ,12672 ,12673 ,12674 ,12675 ,12676 ,12679 ,12681 ,12682 ,12684 ,12685 ,12686 ,12687 ,12696 ,12697 , 12702 ,12704 ,12706 ,12708 ,12712 ,12713 ,12715 ,12724 ,12725 ,12727 ,12730 ,12732 ,12740 ,12741 ,12750 ,12760 ,12762 ,12770 ,12773 ,12788 ,12797 ,12798 ,12799 ,12800 ,12801 ,12804 ,12812 ,12821 ,12822 , 12831 ,12832 ,12833 ,12837 ,12838 ,12840 ,12844 ,12848 ,12849 ,12850 ,12864 ,12866 ,12867 ,12872 ,12873 ,12875 ,12876 ,12883 ,12888 ,12890 ,12892 ,12897 ,12900 ,12904 ,12905 ,12911 ,12930 ,12950 ,12952 , 12953 ,12954 ,12956 ,12957 ,12958 ,12960 ,12964 ,12965 ,12969 ,12970 ,12971 ,12972 ,12979 ,12980 ,12984 ,12985 ,12986 ,12987 ,12988 ,12989 ,12990 ,12991 ,12992 ,12993 ,12994 ,12997 ,13004 ,13005 ,13007 , 13023 ,13026 ,13084 ,13085 ,13088 ,13089 ,13090 ,13091 ,13100 ,13101 ,13102 ,13145 ,13148 ,13149 ,13150 ,13151 ,13154 ,13165 ,13167 ,13174 ,13176 ,13180 ,13181 ,13183 ,13184 ,13185 ,13186 ,13187 ,13188 , 13189 ,13190 ,13192 ,13193 ,13194 ,13195 ,13196 ,13197 ,13199 ,13201 ,13203 ,13204 ,13205 ,13208 ,13218 ,13237 ,13240 ,13241 ,13247 ,13271 ,13274 ,13276 ,13277 ,13280 ,13286 ,13287 ,13288 ,13291 ,13292 , 13298 ,13302 ,13303 ,13304 ,13306 ,13307 ,13308 ,13310 ,13311 ,13320 ,13328 ,13331 ,13333 ,13335 ,13336 ,13337 ,13338 ,13342 ,13344 ,13346 ,13349 ,13350 ,13351 ,13359 ,13360 ,13361 ,13364 ,13366 ,13368 , 13372 ,13374 ,13375 ,13378 ,13380 ,13382 ,13387 ,13391 ,13392 ,13395 ,13396 ,13400 ,13401 ,13402 ,13404 ,13407 ,13408 ,13409 ,13411 ,13413 ,13414 ,13419 ,13426 ,13439 ,13443 ,13448 ,13453 ,13458 ,13460 , 13462 ,13463 ,13466 ,13470 ,13472 ,13474 ,13477 ,13478 ,13482 ,13490 ,13491 ,13492 ,13493 ,13494 ,13495 ,13496 ,13497 ,13498 ,13502 ,13505 ,13508 ,13519 ,13520 ,13530 ,13537 ,13538 ,13542 ,13544 ,13545 , 13556 ,13559 ,13566 ,13567 ,13568 ,13569 ,13570 ,13571 ,13572 ,13575 ,13577 ,13579 ,13580 ,13581 ,13584 ,13586 ,13588 ,13590 ,13595 ,13597 ,13599 ,13601 ,13603 ,13605 ,13606 ,13607 ,13610 ,13612 ,13614 , 13629 ,13630 ,13632 ,13635 ,13636 ,13638 ,13645 ,13648 ,13651 ,13655 ,13657 ,13662 ,13667 ,13673 ,13676 ,13679 ,13683 ,13684 ,13685 ,13686 ,13699 ,13707 ,13711 ,13720 ,13722 ,13723 ,13732 ,13734 ,13735 , 13739 ,13742 ,13746 ,13749 ,13751 ,13777 ,13782 ,13783 ,13786 ,13791 ,13793 ,13797 ,13802 ,13807 ,13808 ,13809 ,13810 ,13811 ,13823 ,13835 ,13840 ,13847 ,13858 ,13859 ,13884 ,13886 ,13905 ,13910 ,13932 , 13936 ,13944 ,13970 ,14003 ,14014 ,14024 ,14031 ,14052 ,14053 ,14054 ,14055 ,14056 ,14057 ,14058 ,14059 ,14060 ,14063 ,14064 ,14065 ,14069 ,14071 ,14072 ,14082 ,14088 ,14092 ,14093 ,14098 ,14099 ,14103 , 14109 ,14112 ,14114 ,14121 ,14123 ,14125 ,14126 ,14127 ,14128 ,14129 ,14132 ,14134 ,14135 ,14190 ,14191 ,14192 ,14193 ,14194 ,14195 ,14196 ,14198 ,14199 ,14200 ,14201 ,14202 ,14203 ,14205 ,14209 ,14213 , 14229 ,14238 ,14249 ,14253 ,14254 ,14258 ,14262 ,14267 ,14268 ,14270 ,14274 ,14277 ,14281 ,14283 ,14284 ,14285 ,14286 ,14287 ,14289 ,14292 ,14293 ,14294 ,14296 ,14297 ,14298 ,14300 ,14301 ,14302 ,14305 , 14307 ,14310 ,14317 ,14319 ,14323 ,14338 ,14342 ,14343 ,14345 ,14349 ,14351 ,14353 ,14358 ,14376 ,14378 ,14386 ,14388 ,14391 ,14393 ,14401 ,14433 ,14437 ,14454 ,14460 ,14464 ,14466 ,14471 ,14478 ,14480 , 14486 ,14488 ,14489 ,14490 ,14491 ,14492 ,14493 ,14498 ,14501 ,14508 ,14509 ,14512 ,14524 ,14526 ,14535 ,14542 ,14550 ,14572 ,14574 ,14579 ,14584 ,14588 ,14594 ,14602 ,14606 ,14607 ,14614 ,14622 ,14657 , 14659 ,14660 ,14667 ,14681 ,14683 ,14684 ,14685 ,14690 ,14691 ,14693 ,14698 ,14714 ,14717 ,14720 ,14730 ,14734 ,14748 ,14751 ,14771 ,14794 ,14804 ,14815 ,14821 ,14833 ,14849 ,14853 ,14863 ,14864 ,14869 , 14873 ,14888 ,14889 ,14898 ,14899 ,14900 ,14902 ,14903 ,15283 ,15414 ,15418 ,15421 ,15425 ,15495 ,15515 ,15562 ,15569 ,15582 ,15587 ,15615 ,15632 ,15644 ,15659 ,15723 ,15731 ,15738 ,15889 ,15890 ,15919 , 15931 ,15933 ,15935 ,15975 ,15979 ,16011 ,16022 ,16023 ,16109 ,16115 ,16116 ,16117 ,16118 ,16119 ,16120 ,16124 ,16125 ,16126 ,16128 ,16129 ,16130 ,16132 ,16138 ,16139 ,16140 ,16141 ,16144 ,16147 ,16149 , 16150 ,16151 ,16152 ,16153 ,16160 ,16161 ,16162 ,16163 ,16173 ,16188 ,16193 ,16198 ,16199 ,16215 ,16216 ,16217 ,16220 ,16225 ,16227 ,16231 ,16232 ,16249 ,16250 ,16255 ,16264 ,16266 ,16270 ,16271 ,16272 , 16283 ,16286 ,16288 ,16301 ,16303 ,16306 ,16313 ,16314 ,16315 ,16316 ,16317 ,16318 ,16319 ,16320 ,16321 ,16323 ,16324 ,16325 ,16326 ,16328 ,16332 ,16334 ,16335 ,16338 ,16339 ,16341 ,16342 ,16343 ,16344 , 16346 ,16347 ,16349 ,16351 ,16352 ,16353 ,16355 ,16356 ,16357 ,16360 ,16361 ,16363 ,16369 ,16370 ,16371 ,16372 ,16373 ,16375 ,16377 ,16378 ,16379 ,16380 ,16383 ,16384 ,16385 ,16386 ,16387 ,16388 ,16390 , 16391 ,16392 ,16394 ,16395 ,16396 ,16397 ,16398 ,16400 ,16401 ,16402 ,16403 ,16404 ,16405 ,16406 ,16407 ,16408 ,16409 ,16410 ,16411 ,16427 ,16428 ,16433 ,16435 ,16436 ,16437 ,16438 ,16442 ,16443 ,16444 , 16445 ,16446 ,16464 ,16472 ,16473 ,16476 ,16477 ,16479 ,16481 ,16490 ,16499 ,16569 ,16576 ,16582 ,16584 ,16585 ,16586 ,16592 ,16593 ,16595 ,16653 ,16654 ,16656 ,16657 ,16658 ,16697 ,16702 ,16729 ,16736 , 16745 ,16757 ,16819 ,16871 ,16874 ,16875 ,16886 ,16889 ,16901 ,16916 ,16921 ,16927 ,16933 ,16945 ,16961 ,16980 ,16982 ,16990 ,17042 ,17101 ,17102 ,17103 ,17116 ,17125 ,17126 ,17127 ,17128 ,17132 ,17133 , 17134 ,17136 ,17137 ,17138 ,17139 ,17140 ,17141 ,17142 ,17143 ,17145 ,17147 ,17148 ,17160 ,17161 ,17162 ,17163 ,17164 ,17165 ,17166 ,17167 ,17168 ,17179 ,17189 ,17191) return(sprintf('PF%05d',X)) } generate.org.mat <- function(nSamples = 6, nOrg = 5, nVars = 1000, nDiffFunc = 100, CountVec = c(1e7,5e6,1e6,5e5,1e5), Sampling.Rate.Mat) { ##use to generate simulated data #nSamples = number of samples per condition #nOrg = number of organisms to simulate #nVars = number of features to generate #nDiffFunc = number of differentially expressed functions #CountVec = base library size for each organism #Sampling.Rate.Mat: library size parameters; rows = samples (i.e. 12 for 6 samples per condition), columns = organisms require(compcodeR) #error checking if (ncol(Sampling.Rate.Mat) != nOrg | nrow(Sampling.Rate.Mat) != nSamples*2 | length(CountVec) != nOrg) { cat('dimension mismatch for samples or organisms, please check.','\n') return(-1) } B.gen <- list() for (k in 1:nOrg) { B.gen[[k]] <- generateSyntheticData(dataset = "B", n.vars = nVars, samples.per.cond = nSamples, n.diffexp = nDiffFunc, repl.id = 1, seqdepth = CountVec[k], fraction.upregulated = 0.5, #fraction.upregulated = 0.01, filter.threshold.total = -1, filter.threshold.mediancpm = -1, minfact = Sampling.Rate.Mat[,k], maxfact = Sampling.Rate.Mat[,k] ) } return(B.gen) } DESeq2.norm.mat <- function(Xmat,cond,type) { #Xmat = raw count data for one species Xmat.col = ncol(Xmat) Xmat.row = nrow(Xmat) colData <- data.frame(condition = cond, type = type) Xmat = round(Xmat) storage.mode(Xmat) <- 'integer' colnames(Xmat) <- as.character(1:Xmat.col) dds <- DESeqDataSetFromMatrix(countData = Xmat, colData = colData, design = ~condition) colData(dds)$condition = factor(colData(dds)$condition, levels=unique(cond)) dds <- DESeq(dds, quiet = TRUE) #normalize the data YMat <- Xmat/rep(dds@colData@listData$sizeFactor, each = (Xmat.row)) print(dds@colData@listData$sizeFactor) return(YMat) } DESeq2.result <- function(Xmat,cond,type) { Xmat.col = ncol(Xmat) Xmat.row = nrow(Xmat) Xmat = round(Xmat) storage.mode(Xmat) <- 'integer' colnames(Xmat) <- as.character(1:Xmat.col) colData <- data.frame(condition = cond, type = type) dds <- DESeqDataSetFromMatrix(countData = Xmat, colData = colData, design = ~condition) colData(dds)$condition = factor(colData(dds)$condition, levels=unique(cond)) #stop DESeq2 from performing additional normalization normFactors <- matrix(1,ncol = Xmat.col, nrow = Xmat.row) normalizationFactors(dds) <- normFactors dds <- DESeq(dds, quiet = TRUE) res <- results(dds) return(res) } abind.matrix <- function(Xarray=NULL,XMat) { require(abind) #Xmat = raw count data for one species #Xmat rows = features, coloums = samples #Xarray = array containing data for already added species Xarray <- abind(Xarray,XMat, along = 3) return(Xarray) } DESeq2.tax.specific <- function(Xarray,cond.vec,type.vec) { require(DESeq2) #Xarray array[feature,sample,organism] #cond.vec = the different condition types of the sample #type.vec = the type of sequence, for example 'single-read' nOrg = dim(Xarray)[3] nSamples = dim(Xarray)[2] #check input if (length(cond.vec) != nSamples | length(type.vec) != nSamples) { return(-1) } Scaled.Mat <- matrix(0,ncol = ncol(Xarray), nrow = nrow(Xarray)) cat("scaling ",nOrg," different matrices.\n") for (i in 1:nOrg) { Scaled.Mat <- Scaled.Mat + DESeq2.norm.mat(Xarray[,,i], cond.vec,type.vec) } return(DESeq2.result(Scaled.Mat,cond.vec,type.vec)) } de.counts.with.direction <- function(Compcoder.List,DESeq.result,p.thresh,STACK.FLAG) { nOrg = length(Compcoder.List) #features are stacked if(STACK.FLAG) { I.ori.DE = c() I.ori.DIR.UP = c() I.ori.DIR.DOWN = c() for (i in 1:nOrg) { I.ori.DE = c(I.ori.DE,Compcoder.List[[i]]@variable.annotations$differential.expression == 1) I.ori.DIR.UP = c(I.ori.DIR.UP, Compcoder.List[[i]]@variable.annotations$upregulation == 1) I.ori.DIR.DOWN = c(I.ori.DIR.DOWN, Compcoder.List[[i]]@variable.annotations$downregulation == 1) } } else { #all have same direction... I.ori.DE = Compcoder.List[[1]]@variable.annotations$differential.expression == 1 I.ori.DIR.UP = Compcoder.List[[1]]@variable.annotations$upregulation == 1 I.ori.DIR.DOWN = Compcoder.List[[1]]@variable.annotations$downregulation == 1 } #set na values to a value... DESeq.result$padj[is.na(DESeq.result$padj)] = 1 I.DE = DESeq.result$padj <= p.thresh I.DIR.UP = DESeq.result$log2FoldChange > 0 I.DIR.DOWN = DESeq.result$log2FoldChange < 0 SUM.DE = sum(I.ori.DE & I.DE) SUM.DE.corre.dir = sum(I.ori.DE & I.DE & I.DIR.UP & I.ori.DIR.UP) + sum(I.ori.DE & I.DE & I.DIR.DOWN & I.ori.DIR.DOWN) SUM.DE.wrong.dir = sum(I.ori.DE & I.DE & I.DIR.UP & !I.ori.DIR.UP) + sum(I.ori.DE & I.DE & I.DIR.DOWN & !I.ori.DIR.DOWN) SUM.NDE = sum(!I.ori.DE & !I.DE) SUM.EXTRA = sum(!I.ori.DE & I.DE) SUM.MISSED = sum(I.ori.DE & !I.DE) ret <- c(SUM.DE,SUM.DE.corre.dir,SUM.DE.wrong.dir,SUM.NDE,SUM.EXTRA,SUM.MISSED) names(ret) <- c('SUM.DE','SUM.DE.corre.dir','SUM.DE.wrong.dir','SUM.NDE','SUM.EXTRA','SUM.MISSED') return(ret) } de.counts.with.direction.per.org <- function(Compcoder.List,DESeq.result,p.thresh) { nOrg = length(Compcoder.List) ret.Mat <- matrix(0,nOrg,6) #features are stacked I.ori.DE.Vec = c() I.ori.DIR.UP.Vec = c() I.ori.DIR.DOWN.Vec = c() for (i in 1:nOrg) { I.ori.DE.Vec = c(I.ori.DE.Vec,Compcoder.List[[i]]@variable.annotations$differential.expression == 1) I.ori.DIR.UP.Vec = c(I.ori.DIR.UP.Vec, Compcoder.List[[i]]@variable.annotations$upregulation == 1) I.ori.DIR.DOWN.Vec = c(I.ori.DIR.DOWN.Vec, Compcoder.List[[i]]@variable.annotations$downregulation == 1) } #set na values to a value... DESeq.result$padj[is.na(DESeq.result$padj)] = 1 #check each org seperatly nFeat.total = length(I.ori.DE.Vec) nFeat = length(Compcoder.List[[i]]@variable.annotations$differential.expression) Start.Vec = seq(1,nFeat.total,nFeat) End.Vec = seq(nFeat,nFeat.total,nFeat) I.DE.Vec = DESeq.result$padj <= p.thresh I.DIR.UP.Vec = DESeq.result$log2FoldChange > 0 I.DIR.DOWN.Vec = DESeq.result$log2FoldChange < 0 for (i in 1:nOrg) { I.DE = I.DE.Vec[Start.Vec[i]:End.Vec[i]] I.DIR.UP = I.DIR.UP.Vec[Start.Vec[i]:End.Vec[i]] I.DIR.DOWN = I.DIR.DOWN.Vec[Start.Vec[i]:End.Vec[i]] I.ori.DE = I.ori.DE.Vec[Start.Vec[i]:End.Vec[i]] I.ori.DIR.UP = I.ori.DIR.UP.Vec[Start.Vec[i]:End.Vec[i]] I.ori.DIR.DOWN = I.ori.DIR.DOWN.Vec[Start.Vec[i]:End.Vec[i]] SUM.DE = sum(I.ori.DE & I.DE) SUM.DE.corre.dir = sum(I.ori.DE & I.DE & I.DIR.UP & I.ori.DIR.UP) + sum(I.ori.DE & I.DE & I.DIR.DOWN & I.ori.DIR.DOWN) SUM.DE.wrong.dir = sum(I.ori.DE & I.DE & I.DIR.UP & !I.ori.DIR.UP) + sum(I.ori.DE & I.DE & I.DIR.DOWN & !I.ori.DIR.DOWN) SUM.NDE = sum(!I.ori.DE & !I.DE) SUM.EXTRA = sum(!I.ori.DE & I.DE) SUM.MISSED = sum(I.ori.DE & !I.DE) ret <- c(SUM.DE,SUM.DE.corre.dir,SUM.DE.wrong.dir,SUM.NDE,SUM.EXTRA,SUM.MISSED) names(ret) <- c('SUM.DE','SUM.DE.corre.dir','SUM.DE.wrong.dir','SUM.NDE','SUM.EXTRA','SUM.MISSED') ret.Mat[i,] = ret } return(ret.Mat) } all.cuts <- function(X.Stacked, Y.Stacked, X.Super, Y.Super, UP.Vec, Down.Vec) { nFac = length(X.Stacked)/length(X.Super) X.Super.rep = rep(X.Super,nFac) Y.Super.rep = rep(Y.Super,nFac) nCuts = 0; X.Stacked[is.infinite(X.Stacked)] = 1 Y.Stacked[is.infinite(Y.Stacked)] = 1 Loc = which(UP.Vec == 1 | Down.Vec == 1) for (i in Loc) { ret = find.cut.with.diagonal(c(X.Stacked[i],Y.Stacked[i]),c(X.Super.rep[i],Y.Super.rep[i])) if (ret) { nCuts = nCuts+1 } } return(nCuts) } find.cut.with.diagonal <- function(PointA,PointB) { Xvals = c(PointA[1],PointB[1]) Yvals = c(PointA[2],PointB[2]) minX = min(Xvals) maxX = max(Xvals) res = lm(Xvals ~ Yvals) res2 = lm(Yvals ~ Xvals) b = res2$coefficients[1] x.cut = res$coefficients[1] m = (PointA[2]-b)/PointA[1] funcA <- function(x) { m * x + b } funcB <- function(x) { x } XVs = seq(minX,maxX,0.01) YVs = funcA(XVs) YVs2 = funcB(XVs) I = YVs < YVs2 I2 = YVs > YVs2 if (sum(I) != length(I) & sum(I2) != length(I2)) { return(TRUE) } return(FALSE) } generate.samp.rates <- function(iOrg,curr.type) { Sampling.Rate.Mat = c() if (iOrg %in% seq(1,1,1)) { ##one Condition up #weaker if (curr.type == 'sim3likeweak') { Sampling.Rate.Mat = c(sample(seq(1.5,2,0.05), size = nSamples, replace = TRUE), sample(seq(0.5,0.667,0.05), size = nSamples, replace = TRUE)) } #stronger else if (curr.type == 'sim3likestrong') { Sampling.Rate.Mat = c(sample(seq(9,10,0.05), size = nSamples, replace = TRUE), sample(seq(0.1,0.3,0.05), size = nSamples, replace = TRUE)) } ##sim2-like else if (curr.type == 'sim2like') { Sampling.Rate.Mat = sample(seq(0.5,2,0.1), size = nSamples*2, replace = TRUE) } # #Sampling.Rate.Mat[,iOrg] = c(sample(seq(1.5,2,0.05), size = nSamples/2, replace = TRUE), sample(seq(0.5,0.667,0.05), size = nSamples/2, replace = TRUE), sample(seq(1.5,2,0.05), size = nSamples/2, replace = TRUE), sample(seq(0.5,0.667,0.05), size = nSamples/2, replace = TRUE)) # } else if (iOrg %in% seq(2,2,1)) { ##one Condition up #weaker if (curr.type == 'sim3likeweak') { Sampling.Rate.Mat = c(sample(seq(0.5,0.667,0.05), size = nSamples, replace = TRUE), sample(seq(1.5,2,0.05), size = nSamples, replace = TRUE)) } #stronger else if (curr.type == 'sim3likestrong') { Sampling.Rate.Mat = c(sample(seq(0.1,0.3,0.05), size = nSamples, replace = TRUE), sample(seq(9,10,0.05), size = nSamples, replace = TRUE)) } ##sim2-like else if (curr.type == 'sim2like') { Sampling.Rate.Mat = sample(seq(0.5,2,0.1), size = nSamples*2, replace = TRUE) } #Sampling.Rate.Mat[,iOrg] = c(sample(seq(0.5,0.667,0.05), size = nSamples/2, replace = TRUE), sample(seq(1.5,2,0.05), size = nSamples/2, replace = TRUE), sample(seq(0.5,0.667,0.05), size = nSamples/2, replace = TRUE), sample(seq(1.5,2,0.05), size = nSamples/2, replace = TRUE)) } else { cat('more than 2 organisms selected, not like the sim4 in the manuscript!') Sampling.Rate.Mat = sample(seq(0.5,2,0.1), size = nSamples*2, replace = TRUE) } return(Sampling.Rate.Mat) } generate.plot.data <- function(Result.List) { nIter = length(Result.List) nOrg = dim(Result.List[[1]]$T.Mat)[1] TP.Mat.tax = matrix(0, nrow = nIter, ncol = nOrg) TP.Mat.glo = TP.Mat.tax FP.Mat.tax = TP.Mat.tax FP.Mat.glo = TP.Mat.tax for (i in 1:nIter) { TP.Mat.tax[i,] = Result.List[[i]]$T.Mat[,2] FP.Mat.tax[i,] = Result.List[[i]]$T.Mat[,5] + Result.List[[i]]$T.Mat[,3] TP.Mat.glo[i,] = Result.List[[i]]$G.Mat[,2] FP.Mat.glo[i,] = Result.List[[i]]$G.Mat[,5] + Result.List[[i]]$G.Mat[,3] } df <- rbind(data.frame(x = as.vector(TP.Mat.tax), Org = rep(1:nOrg, each = nIter), type = rep('TP',nOrg*nIter), method = 'tax'), data.frame(x = as.vector(TP.Mat.glo), Org = rep(1:nOrg, each = nIter), type = rep('TP',nOrg*nIter), method = 'glo'), data.frame(x = as.vector(FP.Mat.tax), Org = rep(1:nOrg, each = nIter), type = rep('FP',nOrg*nIter), method = 'tax'), data.frame(x = as.vector(FP.Mat.glo), Org = rep(1:nOrg, each = nIter), type = rep('FP',nOrg*nIter), method = 'glo') ) return(df) } generate.roc <- function(Results, nIter) { require(ROCR) pp.1 <- list() pp.2 <- list() ll <- list() nFac = dim(Results[[1]]$T.Mat)[1] for (i in 1:nIter) { pp.1[[i]] <- 1-Results[[i]]$res.tax.padj pp.2[[i]] <- 1-Results[[i]]$res.glo.padj I = Results[[i]]$UP.reg | Results[[i]]$DOWN.reg vec = vector(mode = 'numeric', length = length(I)) vec = vec - 1 vec[I] = 1 ll[[i]] <- vec } x11() par(mfrow=c(1,2)) pred<- prediction(pp.1, ll) perf <- performance(pred, "tpr", "fpr") plot(perf, avg= "threshold", colorize=T, lwd= 3, main= "taxon specific scaling") plot(perf, lty=3, col="grey78", add=T) pred<- prediction(pp.2, ll) perf <- performance(pred, "tpr", "fpr") plot(perf, avg= "threshold", colorize=T, lwd= 3, main = "global scaling") plot(perf, lty=3, col="grey78", add=T) } plot.prepare.sim4 <- function(Results,dir = 'UP') { arrow.plot <- list() some.plot <- list() nIter = length(Results) #plot only first... for ( i in 1:nIter) { Ret <- get.positions(Results,i, dir) arrow.plot[[i]] <- generate.arrow.plot(Ret) df <- prepare.plot.first(Result = Results, iter = i) p <- ggplot(df, aes(x, y)) + scale_shape_identity() + scale_size_identity() + geom_point(aes(shape = shape, size = size)) p <- p + facet_wrap(~type, scales = 'free') + geom_abline(intercept = 0, slope = 1, alpha = 0.25, linetype = 'dashed') + theme(legend.position="bottom") + xlab('log counts Cond A') + ylab('log counts Cond B') + ggtitle('superpositioned counts') some.plot[[i]] = p } return(list(arrow.plot = arrow.plot,some.plot = some.plot)) } get.positions <- function(Results,iter, dir) { nFeatures = length(Results[[iter]]$Y.super.tax) if (dir == 'UP') { X.stacked = Results[[iter]]$X.stacked[Results[[iter]]$UP.reg == 1] Y.stacked = Results[[iter]]$Y.stacked[Results[[iter]]$UP.reg == 1] X.super.tax = Results[[iter]]$X.super.tax[Results[[iter]]$UP.reg[1:nFeatures] == 1] Y.super.tax = Results[[iter]]$Y.super.tax[Results[[iter]]$UP.reg[1:nFeatures] == 1] X.super.glo = Results[[iter]]$X.super.glo[Results[[iter]]$UP.reg[1:nFeatures] == 1] Y.super.glo = Results[[iter]]$Y.super.glo[Results[[iter]]$UP.reg[1:nFeatures] == 1] } else { X.stacked = Results[[iter]]$X.stacked[Results[[iter]]$DOWN.reg == 1] Y.stacked = Results[[iter]]$Y.stacked[Results[[iter]]$DOWN.reg == 1] X.super.tax = Results[[iter]]$X.super.tax[Results[[iter]]$DOWN.reg[1:nFeatures] == 1] Y.super.tax = Results[[iter]]$Y.super.tax[Results[[iter]]$DOWN.reg[1:nFeatures] == 1] X.super.glo = Results[[iter]]$X.super.glo[Results[[iter]]$DOWN.reg[1:nFeatures] == 1] Y.super.glo = Results[[iter]]$Y.super.glo[Results[[iter]]$DOWN.reg[1:nFeatures] == 1] } ret <- list(X.stacked,Y.stacked,X.super.tax,Y.super.tax,X.super.glo,Y.super.glo) names(ret) <- c('X.stacked','Y.stacked','X.super.tax','Y.super.tax','X.super.glo','Y.super.glo') return(ret) } prepare.plot.first <- function(Results,iter,type) { nFeatures = length(Results[[iter]]$Y.super.tax) shape.Vec = rep(16,nFeatures) shape.Vec[Results[[iter]]$DOWN.reg[1:nFeatures] == 1 | Results[[iter]]$UP.reg[1:nFeatures] == 1] = 1 size.Vec = rep(0.5,nFeatures) size.Vec[Results[[iter]]$DOWN.reg[1:nFeatures] == 1 | Results[[iter]]$UP.reg[1:nFeatures] == 1] = 1.5 tmp.df.1 <- data.frame(x = Results[[iter]]$X.super.tax, y = Results[[iter]]$Y.super.tax, shape = shape.Vec, size = size.Vec, type = 'tax' ) tmp.df.2 <- data.frame(x = Results[[iter]]$X.super.glo, y = Results[[iter]]$Y.super.glo, shape = shape.Vec, size = size.Vec, type = 'glo' ) df <- rbind.data.frame(tmp.df.1,tmp.df.2) df$type <- factor(df$type, levels = c('glo','tax')) return(df) } generate.arrow.plot <- function(Ret) { nOrg = length(Ret$X.stacked)/length(Ret$X.super.glo) Start.Vec = seq(1,length(Ret$X.stacked),length(Ret$X.super.glo)) End.Vec = seq(length(Ret$X.super.glo),length(Ret$X.stacked),length(Ret$X.super.glo)) comb.df <- data.frame() df.tax <- data.frame() df.glo <- data.frame() for (i in 1:length(Start.Vec)) { df.tmp <- data.frame(x1 = Ret$X.super.tax, x2 = Ret$X.stacked[Start.Vec[i]:End.Vec[i]], y1 = Ret$Y.super.tax, y2 = Ret$Y.stacked[Start.Vec[i]:End.Vec[i]], col = i, type = 'tax') comb.df <- rbind.data.frame(comb.df,df.tmp) df.tmp <- data.frame(x1 = Ret$X.super.glo, x2 = Ret$X.stacked[Start.Vec[i]:End.Vec[i]], y1 = Ret$Y.super.glo, y2 = Ret$Y.stacked[Start.Vec[i]:End.Vec[i]], col = i, type = 'glo') comb.df <- rbind.data.frame(comb.df,df.tmp) } p <- ggplot() for (i in 1:length(Start.Vec)) { p <- p + geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2, color = factor(col)), data = comb.df, arrow = arrow(length = unit(0.3, "cm"), ends = 'first')) p <- p + geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2, color = factor(col)), data = comb.df , arrow = arrow(length = unit(0.15, "cm"), ends = 'last', angle = 90)) } return(p + geom_abline(intercept = 0, slope = 1, alpha = 0.25, linetype = 'dashed') + facet_wrap(~type) + theme(legend.position="bottom") + xlab('log count Cond A') + ylab('log count Cond B') + ggtitle('transcriptome counts to superopsitioned counts')) }