v 1.2

data/plots.R

@@ -0,0 +1,41 @@
+
+
+
+rPar <- rainbowPlot(parrotDF, seed=12345)
+# remove multiple outliers
+toRemove <- interceptPositionPlot(rPar, highlightOutliers=T)
+parrotDF2 <- parrotDF[!parrotDF$Position %in% toRemove,]
+rainbowPlot(parrotDF2, seed=12345)
+
+
+rHum <- rainbowPlot(humanDF, seed=12345)
+# automatic outlier detection would remove SNPs that look OK
+remHum <- interceptPositionPlot(rHum, highlightOutliers = T)
+
+# select single outlier "by hand"
+which.max(coef(summary(rHum$lmer.model))[,1])
+
+humanDF2 <- humanDF[humanDF$Position != 310,] # remove "by hand"
+rainbowPlot(humanDF2, seed=12345)
+divEst(parrotFX)
+
+# generate plots for manuscript
+#pdf("hopper.pdf", width=6, height=7)
+
+#download.file("https://tinyurl.com/4mtrbkzc", destfile = "hopper.csv")
+hopperDF <- read.table("hopper.csv")
+
+hopperFit <- rainbowPlot(hopperDF,
+                         seed = 12345, printout = FALSE, title = "Grasshopper")
+#dev.off()
+
+
+#pdf("parrot.pdf", width=6, height=7)
+rainbowPlot(parrotDF2,
+            seed = 12345, printout = FALSE, title = "Parrot")
+#dev.off()
+
+#pdf("human.pdf", width=6, height=7)
+rainbowPlot(humanDF2,
+            seed = 12345, printout = FALSE, title = "Human")
+#dev.off()

vagrantDNA/DESCRIPTION

@@ -1,7 +1,7 @@
 Package: vagrantDNA
 Type: Package
 Title: Estimating the proportion of vagrant DNA in a genome
-Version: 1.1.1
+Version: 1.2.0
 Author: Richard Nichols and Hannes Becher
 Maintainer: The package maintainer <[email protected]>
 Description: The package two functions rainbowPlot and divEst. They produce two different estimates of the proportion of DNA from a particular vagrant genome. They exploit low coverge data from multiple individuals in which the vagrant DNA has taken residence (see Becher and Nichols 2022).

vagrantDNA/R/Rainbow_Plot.R

@@ -74,7 +74,7 @@
 #'  download.file("https://tinyurl.com/4mtrbkzc", destfile = "hopper.csv")
 #'  hopperDF <- read.table("hopper.csv")}
 #'
-#' ## (the t.ly/ link is to the cvs file at
+#' ## (the t.ly/ link is to the CSV file at
 #' ## https://raw.githubusercontent.com/SBCSnicholsLab/pseudogene_quantification/
 #' ## main/data/grasshopper/transformedData.csv)
 #'

vagrantDNA/R/data.R

@@ -4,7 +4,7 @@
 #' sequencing data mapped against mitochondrial genome references.
 #'
 #' @format Each is a data.frame.
-#' [species]DF just contain variant calling data.
+#' \strong{\code{[species]DF}} just contain variant calling data.
 #' The columns are:
 #' \describe{
 #'   \item{Sample}{A character vector giving a unique name for each sample}
@@ -15,7 +15,7 @@
 #'   \item{nTot}{A numeric vector giving the total number of base pairs of the sample's mapping data (ideally after quality control, filtering, read trimming, ect.)}
 #'   \item{ylog}{A numeric vector giving the log of `AltProp`}
 #' }
-#' [species]FX were generated from species with diverged populations. These
+#' \strong{\code{[species]FX}} were generated from species with diverged populations. These
 #' contain information on sites with fixed differences. The columns are:
 #' \describe{
 #'   \item{pos}{A vector giving the site IDs}
@@ -32,6 +32,29 @@
 #'   }
 #'
 #' @examples
+#' # Generate the estimates reported in the paper
+#' \dontrun{
+#' ## Parrot
+#' rPar <- rainbowPlot(parrotDF, seed=12345)
+#' # remove multiple outliers
+#' toRemove <- interceptPositionPlot(rPar, highlightOutliers=T)
+#' parrotDF2 <- parrotDF[!parrotDF$Position %in% toRemove,]
+#' rainbowPlot(parrotDF2, seed=12345, title = "Parrot")
+#'
+#' ## Human
+#' rHum <- rainbowPlot(humanDF, seed=12345)
+#' # automatic outlier detection would remove SNPs that look OK
+#' interceptPositionPlot(rHum, highlightOutliers = T)
+#' # select single outlier "by hand"
+#' which.max(coef(summary(rHum$lmer.model))[,1])
+#' humanDF2 <- humanDF[humanDF$Position != 310,] # by hand, better
+#' rainbowPlot(humanDF2, seed=12345, title = "Human")
+#'
+#' ## Grasshopper
+#' download.file("https://tinyurl.com/4mtrbkzc", destfile = "hopper.csv")
+#' hopperDF <- read.table("hopper.csv")
+#' rainbowPlot(hopperDF, seed = 12345, title = "Grasshopper")
+#' }
 #' humanDF
 #' parrotDF
 #' hopperFX

vagrantDNA/R/interceptPositionPlot.R

@@ -1,27 +1,47 @@
 
-#' Intercept-position-plot for loci seect by rainbowPlot()
+#' Intercept-position-plot for loci selected by \code{rainbowPlot()}
 #'
-#' @param ff A `rainbowPlot()` fit.
+#' @param ff A \code{rainbowPlot()} fit.
 #' @param main Plot caption, default is "Intercept-position-plot"
-#' @param ... Arguments to be passed to `plot()`
+#' @param highlightOutliers Logical. Whether or not to highlight sites with
+#' outlier intercepts according to \code{boxplot.stats}. Default \code{F}
+#' @param ... Arguments to be passed to \code{plot()}
+#' @return Returns invisibly a vector of positions where the intercept estimate
+#' were outliers according to \code{boxplot.stats}
 #'
 #' @export
 #'
 #' @examples
 #'  \dontrun{
-#'  humanFit <- rainbowPlot(humanDF)
-#'  interceptPositionPlot(humanFit)
+#'  rPar <- rainbowPlot(parrotDF)
+#'  toRemove <- interceptPositionPlot(rPar, highlightOutliers=T)
+#'  parrotDF2 <- parrotDF[!parrotDF$Position %in% toRemove,]
+#'  rainbowPlot(parrotDF2)
 #'  }
-interceptPositionPlot <- function(ff, main="Intercept-position-plot", ...){
+interceptPositionPlot <- function(ff,
+                                  highlightOutliers=F,
+                                  main="Intercept-position-plot",
+                                  ...){
   coefs.a <- as.data.frame(coef(summary(ff$lmer.model)))
   coefs.a$site.numeric <- as.numeric(substr(rownames(coefs.a), 9, 30))
+  bounds <- boxplot.stats(coefs.a$Estimate)$stat
+  coefs.a$out <- coefs.a$Estimate < bounds[1] | coefs.a$Estimate > bounds[5]
+  if(highlightOutliers) {
   plot(Estimate ~ site.numeric,
        data=coefs.a,
+       col=out+1,
        main=main,
        ...)
+  } else {
+    plot(Estimate ~ site.numeric,
+         data=coefs.a,
+         main=main,
+         ...)
+  }
   l01 <- loess(Estimate ~ site.numeric,
                data=coefs.a)
   lines(predict(l01, coefs.a$site.numeric)~coefs.a$site.numeric)
+  return(invisible(coefs.a$site.numeric[coefs.a$out]))
 }