Merge branch 'dev'

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: phaseless
 Title: Admixture and imputation for low coverage sequencing data in one goal
-Version: 0.3.0
+Version: 0.4.0
 Authors@R: 
     person("Zilong", "Li", , "[email protected]", role = c("aut", "cre", "cph"),
            comment = c(ORCID = "0000-0001-5859-2078"))

Makefile

@@ -3,7 +3,7 @@ CXX      = g++
 
 # CXXFLAGS = -std=c++17 -Wall -O3 -g -fsanitize=address
 # CXXFLAGS = -std=c++17 -Wall -O3 -march=native -DNDEBUG
-CXXFLAGS = -std=c++17 -Wall -O3 -mavx2 -fPIC -DNDEBUG
+CXXFLAGS = -std=c++17 -Wall -O3  -mavx2
 INC      = -I./src -I./inst/include -I$(HTSDIR)
 LDFLAGS  =  -L$(HTSDIR) -Wl,-rpath,$(HTSDIR)
 LIBS     = $(HTSDIR)/libhts.a -llzma -lbz2 -lm -lz -lpthread

NAMESPACE

@@ -1,5 +1,10 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(plot,admixQ)
+S3method(plot,gamma)
+S3method(plot,hapfreq)
+export(admix.alignKStephens)
+export(admix.plotQ)
 export(parse_impute_opt)
 export(parse_impute_par)
 export(parse_joint_par)

R/plot_admixture.R

@@ -0,0 +1,83 @@
+#' @export
+admix.alignKStephens <- function(qlist){
+  require(label.switching)
+  K <- unique(sapply(qlist, ncol))
+  if(length(K) > 1) stop("K in qlist should be the same")
+  # if there is only 1 run, just return it
+  if(length(qlist)==1) {
+    qlist1 <- qlist
+  } else {
+    # if num of inds or K differ, throw error
+    ninds <- unique(sapply(qlist, nrow))
+    if(length(ninds) > 1) stop("number of inds in qlist should be the same")
+    # if all runs have K=1, just return it
+    if(K==1){
+      qlist1 <- qlist
+    } else {
+      qmat <- lapply(qlist,as.matrix)
+      p <- aperm(simplify2array(qmat), c(3,1,2))
+      perm <- label.switching::stephens(p)
+      # reorder and rename columns
+      qlist1 <- lapply(seq_len(dim(p)[1]),
+                       function(i) {
+                         q_perm <- qmat[[i]][, perm$permutations[i,,drop=FALSE],drop=FALSE]
+                         q_perm <- as.data.frame(q_perm)
+                         attributes(q_perm) <- attributes(qlist[[i]])
+                         q_perm
+                       }
+      )
+      names(qlist1) <- names(qlist)
+    }
+  }
+  return(qlist1)
+}
+
+#' @export
+admix.plotQ <- function(qlist, pop, sortind = TRUE, cluster = 1, debug = FALSE, ...) {
+  N <- length(qlist)
+  K <- unique(sapply(qlist, ncol))
+  nind <- unique(sapply(qlist, nrow))
+  par(mfrow=c(N, 1))
+  sortQ <- function(Q, pop) {
+    lapply(split(Q, pop), function(p) {
+      ord <- order(p[,cluster])
+      p[ord,]
+    })
+  }
+  for(i in seq_along(qlist)){
+    Q <- qlist[[i]]
+    if(sortind) {
+      s <- sortQ(Q, pop)
+      ordpop <- order(pop)
+      namepop <- names(s)
+      Q <- t(do.call(rbind, s))
+    } else {
+      ordpop <- 1:length(pop)
+      namepop <- unique(pop)
+      Q <- t(Q)
+      colnames(Q) <- pop
+    }
+    if(i == N && !debug) {
+      med<- tapply(1:nind,pop[ordpop],median)
+      ord <- match(names(med), namepop)
+      groups <- rep(NA, ncol(Q))
+      groups[as.integer(med)] <- namepop[ord]
+      colnames(Q) <- groups
+    } else if (N!=1) {
+      colnames(Q) <- NULL
+    }
+    par(mar=c(3.1,5.1,3.1,1.1))
+    h <- barplot(Q, col = 1:K, border = NA, space = 0, ylab = "Admixture proportion", main = names(qlist)[i], xaxs = "i", ...)
+    abline(v=tapply(h,pop[order(pop)],max)+0.5,col="black",lwd=2,lty=2)
+  }
+}
+
+#' @export
+plot.admixQ <- function(qfiles, pop, ...) {
+  qlist <- lapply(qfiles, read.table)
+  names(qlist) <- names(qfiles)
+  pop <- read.table(pop)[,1]
+  a <- admix.alignKStephens(qlist)
+  admix.plotQ(a, pop, cex.main = 3, cex.lab=3, cex.names = 3,...)
+}
+

R/plot_haplotypes.R

@@ -0,0 +1,110 @@
+#' @export
+plot.gamma <- function(gamma, sites = NULL, ...) {
+  stopifnot(is.list(gamma))
+  N <- length(gamma)
+  C <- nrow(gamma[[1]])
+  M <- ncol(gamma[[1]])
+  if(!is.null(sites) & is.vector(sites) & length(sites) < M) {
+    M <- length(sites)
+  } else {
+    sites <- 1:M
+  }
+  plot(0, 0, col = "white", axes=FALSE, xlim = c(0, M), ylim = c(1, N + 1),...)
+  d <- 1
+  xleft <- 1:M - d
+  xright <- 1:M - d
+  for (i in seq(N)) {
+    ytop <- i + array(0, M)
+    ybottom <- i + array(0, M)
+    for(c in 1:C) {
+      ytop <- ytop + gamma[[i]][c, sites]
+      rect(xleft = xleft - d,  xright = xright + d, ybottom = ybottom, ytop = ytop, col = c, lwd = 0, border = NA)
+      ybottom <- ytop
+    }
+  }
+}
+
+
+#' @export
+plot.hapfreq <- function(hapfreq,
+                         pos,
+                         recomb = NULL,
+                         colors = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7"),
+                         ...) {
+  stopifnot(is.matrix(hapfreq), is.vector(pos))
+  nCols <- length(colors)
+  nGrids <- length(pos)
+  K <- nrow(hapfreq)
+  sum <- array(0, nGrids)
+  xlim <- range(pos)
+  ylim <- c(0, 1)
+  ## OK so if there are grids, use the grid points
+  plot(x = 0, y = 0, xlim = xlim, ylim = ylim, axes = FALSE,  ...)
+  x <- c(pos[1], pos, pos[length(pos):1])
+  m <- array(0, c(nGrids, K + 1))
+  for(i in 1:K) {
+    m[, i + 1] <- m[, i] + hapfreq[i, ]
+  }
+  for(i in K:1) {
+    polygon(
+      x = x, y = c(m[1, i], m[, i + 1], m[nGrids:1, i]),
+      xlim = xlim, ylim = ylim, col = colors[(i %% nCols) + 1]
+    )
+  }
+  if(!is.null(recomb))
+    lines(pos[-1], recomb, type = "l", col = "red")
+}
+
+
+table(dy <- as.matrix(read.table("~/Downloads/chr16.Mkomazi.males.Ychr.txt")))
+
+table(da <- as.matrix(read.table("~/Downloads/chr16.Mkomazi.males.Achr.txt")))
+
+image(da)
+
+
+N <- 4
+M <- 7570
+d <- 1
+xleft <- 1:M - d
+xright <- 1:M - d
+
+samples <- c("NGrataE03715", "NGrataE03718", "NGrataE03717", "NGrataE03723")
+
+f <- function(da){
+  for(i in 1:N){
+    ybottom <- i + array(0, M)
+    ytop <- ybottom+1
+    rect(xleft = xleft - d,  xright = xright + d, ybottom = ybottom, ytop = ytop, col = da[i,]+1, lwd = 3, border = NA)
+    rect(xleft = xleft - d,  xright = xright + d, ybottom = ytop-0.01, ytop = ytop, col = "gray", lwd = 3, border = NA)
+    mtext(samples[i], side = 2, at = (ytop[1]+ybottom[1])/2, cex = 1.5)
+  }
+}
+
+##par(mfrow = c(2,1), cex.lab = 2, cex.main = 2)
+
+op <- par(mfrow = c(2,1),
+          oma = c(6,0,0,0) + 0.1,
+          mar = c(0,4,2,1) + 0.1,
+          cex.lab = 2, cex.main=2)
+plot(0, 0, col = "white", axes=FALSE, xlim = c(0, M), ylim = c(1, N + 1), main = "Y chromosome", xlab="",ylab="")
+f(dy)
+plot(0, 0, col = "white", axes=FALSE, xlim = c(0, M), ylim = c(1, N + 1), main = "A chromosome", xlab="SNP index", ylab="")
+f(da)
+
+
+mycols <- c("gray", "black")
+
+
+
+rect(xleft = xleft - d,  xright = xright + d, ybottom = ybottom, ytop = ytop, col = mycols[da[i,]+1], lwd = 2, border = NA)
+
+
+
+p <- read.table("~/Downloads/fypa.pyfa.pos")
+
+dev.off()
+
+
+plot.hapfreq(d,pos[1:ncol(d)], colors=1:2)
+plot.hapfreq(d, 1:ncol(d), colors=1:2)