From 4ef598c18965f8410731aa208cde2a9fcda7c378 Mon Sep 17 00:00:00 2001
From: zachary-foster <zacharyfoster1989@gmail.com>
Date: Tue, 4 Feb 2020 15:38:08 -0800
Subject: [PATCH] draft amplicon extraction code

---
 DESCRIPTION             |  2 +-
 NAMESPACE               |  1 +
 R/amplicon_extraction.R | 75 +++++++++++++++++++++++++++++++++++++++++
 man/calc_amplicons.Rd   | 46 +++++++++++++++++++++++++
 4 files changed, 123 insertions(+), 1 deletion(-)
 create mode 100644 R/amplicon_extraction.R
 create mode 100644 man/calc_amplicons.Rd

diff --git a/DESCRIPTION b/DESCRIPTION
index 2b1746a4..0d64731d 100755
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -60,7 +60,7 @@ Suggests:
     BiocManager,
     phyloseq
 VignetteBuilder: knitr
-RoxygenNote: 7.0.1
+RoxygenNote: 7.0.2.9000
 Date: 2019-07-17
 Encoding: UTF-8
 biocViews: 
diff --git a/NAMESPACE b/NAMESPACE
index bfcf5d24..5c3aab59 100755
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -3,6 +3,7 @@
 S3method(heat_tree,Taxmap)
 export("%>%")
 export(as_phyloseq)
+export(calc_amplicons)
 export(calc_group_mean)
 export(calc_group_median)
 export(calc_group_rsd)
diff --git a/R/amplicon_extraction.R b/R/amplicon_extraction.R
new file mode 100644
index 00000000..982aa3ea
--- /dev/null
+++ b/R/amplicon_extraction.R
@@ -0,0 +1,75 @@
+#' Extract expected amplicons
+#'
+#' Extract expected amplicons from sequences using primer sequences.
+#'
+#' @param seqs The DNA sequences of the target to be amplified. Letter case is ignored but
+#'   preserved.
+#' @param forward The DNA sequence of the forward primer in the 5' to 3' orientation. IUPAC codes
+#'   can be present.
+#' @param reverse The DNA sequence of the forward primer in the 5' to 3' orientation. IUPAC codes
+#'   can be present.
+#' @param max_mismatch The maximum number of mismatching bases allowed. See
+#'   [Biostrings::matchPattern] for more information.
+#' @param min_mismatch The minimum number of mismatching bases allowed. See
+#'   [Biostrings::matchPattern] for more information.
+#' @param allow_indels If TRUE then indels are allowed. See [Biostrings::matchPattern] for more
+#'   information.
+#' @param algorithm One of: "auto", "naive-exact", "naive-inexact", "boyer-moore" or "shift-or". See
+#'   [Biostrings::matchPattern] for more information.
+#' @param infer_amps If TRUE, performe an additional step to find amplicons in sequences that might
+#'   not contain the primer binding sites due to being produced with similar primers. This is done
+#'   by aligning amplicons found using primers to sequences not found to contain the primers. For
+#'   each sequence, the region aligned to the best-matching amplicon will be included as an amplion
+#'   if:
+#'   *  The aligned region extends to the boundry of the sequence, suggesting the missing primer site is not in the sequence. 
+#'   *  If the aligned region covers the minimum proportion of the best matching amplicon defined by `infer_thresh_cov`.
+#'   *  If the aligned region has the minimum percent identity to the best matching amplicon defined by `infer_thresh_pid`.
+#' @param infer_thresh_cov See the `infer_amps` option description.
+#' @param infer_thresh_pid See the `infer_amps` option description.
+#' 
+#' @export
+calc_amplicons <- function(seqs, forward, reverse, min_coverage = 0.9, ...) {
+  
+  get_amplicon_chr <- function(seq) {
+    map2_chr(seq@ranges@start, seq@ranges@width, function(s, w) {
+      as.character(seq@subject[seq(from = s, length.out = w)])
+    })
+  }
+  
+  
+  # Get simulated amplicons using primers
+  full_amps <- future_map(seqs, function(s) {
+    get_amplicon_chr(matchProbePair(DNAString(s), Fprobe = forward, Rprobe = reverse))
+  })
+  
+  # Check for multiple possible amplicons per input
+  if (any(map_int(full_amps, length) > 1)) {
+    stop('Some inputs have more than one amplicons.')
+  }
+  full_amps <- unlist(full_amps[map_int(full_amps, length) == 1])
+  
+  # Remove primers from amplicons
+  full_amps <- substr(full_amps, start = nchar(forward) + 1, stop = nchar(full_amps) - nchar(reverse))
+  
+  # Align unamplified sequences with best matching amplicons
+  unamped <- seqs[! names(seqs) %in% names(full_amps)]
+  unamped_aligned <- future_map_chr(unamped, function(s) {
+    aligned <- pairwiseAlignment(pattern = full_amps, subject = s, 
+                                 type = 'global-local',
+                                 # type = 'overlap',
+                                 gapOpening = 10, gapExtension = 4)
+    best_score_i <- which.max(aligned@score)
+    best_align <- aligned[best_score_i]
+    aligned_subject_char <- gsub(as.character(best_align@subject), pattern = '-', replacement = '')
+    aligned_amp_char <- gsub(as.character(best_align@pattern), pattern = '-', replacement = '')
+    if (best_align@pattern@range@width < best_align@pattern@unaligned@ranges@width * min_coverage || nchar(aligned_subject_char) == 0) {
+      aligned_subject_char <- NA
+    }
+    aligned_subject_char
+  })
+  
+  # Add full amplicons to inferred amplicons from alignment
+  all_amps <- c(full_amps, unamped_aligned[! is.na(unamped_aligned)])
+  all_amps <- all_amps[is.na(all_amps) | nchar(all_amps) > 0]
+  setNames(all_amps[names(seqs)], names(seqs))
+}
\ No newline at end of file
diff --git a/man/calc_amplicons.Rd b/man/calc_amplicons.Rd
new file mode 100644
index 00000000..56513b4e
--- /dev/null
+++ b/man/calc_amplicons.Rd
@@ -0,0 +1,46 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/amplicon_extraction.R
+\name{calc_amplicons}
+\alias{calc_amplicons}
+\title{Extract expected amplicons}
+\usage{
+calc_amplicons(seqs, forward, reverse, min_coverage = 0.9, ...)
+}
+\arguments{
+\item{seqs}{The DNA sequences of the target to be amplified. Letter case is ignored but
+preserved.}
+
+\item{forward}{The DNA sequence of the forward primer in the 5' to 3' orientation. IUPAC codes
+can be present.}
+
+\item{reverse}{The DNA sequence of the forward primer in the 5' to 3' orientation. IUPAC codes
+can be present.}
+
+\item{max_mismatch}{The maximum number of mismatching bases allowed. See
+[Biostrings::matchPattern] for more information.}
+
+\item{min_mismatch}{The minimum number of mismatching bases allowed. See
+[Biostrings::matchPattern] for more information.}
+
+\item{allow_indels}{If TRUE then indels are allowed. See [Biostrings::matchPattern] for more
+information.}
+
+\item{algorithm}{One of: "auto", "naive-exact", "naive-inexact", "boyer-moore" or "shift-or". See
+[Biostrings::matchPattern] for more information.}
+
+\item{infer_amps}{If TRUE, performe an additional step to find amplicons in sequences that might
+not contain the primer binding sites due to being produced with similar primers. This is done
+by aligning amplicons found using primers to sequences not found to contain the primers. For
+each sequence, the region aligned to the best-matching amplicon will be included as an amplion
+if:
+*  The aligned region extends to the boundry of the sequence, suggesting the missing primer site is not in the sequence. 
+*  If the aligned region covers the minimum proportion of the best matching amplicon defined by `infer_thresh_cov`.
+*  If the aligned region has the minimum percent identity to the best matching amplicon defined by `infer_thresh_pid`.}
+
+\item{infer_thresh_cov}{See the `infer_amps` option description.}
+
+\item{infer_thresh_pid}{See the `infer_amps` option description.}
+}
+\description{
+Extract expected amplicons from sequences using primer sequences.
+}
