deleted file mode 100644

@@ -1,134 +0,0 @@

-#' Align mutliple sequences

-#' 

-#' Perform multiple sequence alignment using one of three methods and output results to the 

-<<<<<<< HEAD

-#' console or as a pdf file.  One may perform the alignment of all amino acid or nucleotide

-#' sequences in a single sample.  Alternatively, one may search for a given sequence 

-=======

-#' consule or as a pdf file.  One may perform the alignment of all amino acid or nucleotide

-#' sequences in a single sample  Alternatively, one may search for a given sequence 

->>>>>>> Added ability to print alignment to consule

-#' within a list of samples using an edit distance threshold.

-#' 

-#' @param list A list of data frames consisting of antigen receptor sequences 

-#' imported by the LymphoSeq function readImmunoSeq.

-#' @param sample A character vector indicating the name of the sample in the productive

-#' sequence list. 

-#' @param sequence A character vector of one ore more amino acid or nucleotide 

-#' CDR3 sequences to search.

-#' @param editDistance An integer giving the minimum edit distance that the 

-#' sequence must be less than or equal to.  See details below.

-#' @param type A character vector indicating whether "aminoAcid" or "nucleotide" sequences

-#' should be aligned.  If "aminoAcid" is specified, then run productiveSeqs first.

-#' @param method A character vector indicating the multiple sequence alignment method to 

-#' be used.  Refer to the Bioconductor msa package for more details.  Options incude 

-#' "ClustalW", "ClustalOmega", and "Muscle".

-#' @param output A character vector indicating where the multiple sequence alignemnt should be

-<<<<<<< HEAD

-#' printed.  Options include "console" or "pdf".  If "pdf" is selected, the file is saved to

-=======

-#' printed.  Options include "consule" or "pdf".  If "pdf" is selected, the file is saved to

->>>>>>> Added ability to print alignment to consule

-#' the working directory.  For "pdf" to work, Texshade must be installed.  Refer to the 

-#' Bioconductor package msa installation instructions for more details.

-#' @details Edit distance is a way of quantifying how dissimilar two sequences 

-#' are to one another by counting the minimum number of operations required to 

-#' transform one sequence into the other.  For example, an edit distance of 0 

-#' means the sequences are identical and an edit distance of 1 indicates that 

-#' the sequences different by a single amino acid or nucleotide.

-#' 

-<<<<<<< HEAD

-#' @return Performs a multiple sequence alignemnt and prints to the console or saves a pdf to 

-=======

-#' @return Performs a multiple sequence alignemnt and prints to the consule or saves a pdf to 

->>>>>>> Added ability to print alignment to consule

-#' the working directory.

-#' @seealso If having trouble saving pdf files, refer to Biconductor package msa for

-#' installation instructions

-#' \url{http://bioconductor.org/packages/release/bioc/vignettes/msa/inst/doc/msa.pdf}

-#' @examples

-#' file.path <- system.file("extdata", "IGH_sequencing", package = "LymphoSeq")

-#' 

-#' file.list <- readImmunoSeq(path = file.path)

-#' 

-#' productive.nt <- productiveSeq(file.list = file.list, aggregate = "nucleotide")

-#' 

-#' alignSeq(list = productive.nt, sample = "IGH_MVQ92552A_BL", type = "nucleotide", 

-<<<<<<< HEAD

-#'          method = "ClustalW", output = "console")

-=======

-#'          method = "ClustalW", output = "consule")

->>>>>>> Added ability to print alignment to consule

-#' @export

-#' @importFrom Biostrings DNAStringSet

-#' @importFrom Biostrings AAStringSet

-#' @import msa

-<<<<<<< HEAD

-alignSeq = function(list, sample = NULL, sequence = NULL, editDistance = 15, output = "console", type = "nucleotide", method = "ClustalOmega") {

-=======

-alignSeq = function(list, sample = NULL, sequence = NULL, editDistance = 15, output = "consule", type = "nucleotide", method = "ClustalOmega") {

->>>>>>> Added ability to print alignment to consule

-    if(!is.null(sequence) & is.null(sample)){

-        file <- searchSeq(list = list, sequence = sequence, editDistance = editDistance, type = type, match = "partial")

-        if(is.null(file)){

-            stop("There are no sequences to be aligned", call. = FALSE)

-        }

-    }

-    if(!is.null(sequence) & !is.null(sample)){

-        file <- searchSeq(list = list[sample], sequence = sequence, editDistance = editDistance, type = type, match = "partial")

-        if(is.null(file)){

-            stop("There are no sequences to be aligned", call. = FALSE)

-        }

-    }

-    if(is.null(sequence)){

-        file <- list[[sample]] 

-    }

-    if(nrow(file) > 150){

-        file <- file[1:150,]

-        message("Only the first 150 sequences will be aligned")

-    }

-    file[file == ""] <- "unresolved"

-    if(type == "nucleotide"){

-        names(file)[which(names(file) == "foundSequnece")] <- "nucleotide"

-        file <- file[nchar(file$nucleotide) > 3, ]

-        string <- Biostrings::DNAStringSet(file$nucleotide)

-    }

-    if(type == "aminoAcid"){

-        names(file)[which(names(file) == "foundSequnece")] <- "aminoAcid"

-        file <- file[nchar(file$aminoAcid) > 3, ]

-        string <- Biostrings::AAStringSet(file$aminoAcid)

-    }

-    if(nrow(file) < 3){

-        stop("There are less than 3 sequences to be aligned", call. = FALSE)

-    }

-    if(!is.null(sequence)){

-        names(string) <- paste(file$sample)

-    } else {

-        names(string) <- paste(file$vFamilyName, file$dFamilyName, file$jFamilyName, file$count)

-    }

-        names(string) <- gsub(names(string), pattern = "IGH|IGL|IGK|TCRB|TCRA", replacement = "")

-        names(string) <- gsub(names(string), pattern = "unresolved", replacement = "UNR")

-        alignment <- msa::msa(string, method = method)

-<<<<<<< HEAD

-        if(output == "console"){

-=======

-        if(output == "consule"){

->>>>>>> Added ability to print alignment to consule

-            print(alignment, show = "complete")

-        }

-        if(output == "pdf"){

-            if(!is.null(sample)){

-                file.name <- paste(names(list[sample]), ".pdf", sep = "")

-            } else {

-                file.name <- "Sequence_aligment.pdf"

-            }

-            msa::msaPrettyPrint(alignment, output = "pdf", file = file.name, 

-                                showNumbering = "right", showNames = "left", showConsensus = "top",  

-                                shadingMode = "similar", shadingColors = "blues",

-                                paperWidth = ncol(alignment)*0.1 + 5, paperHeight = nrow(alignment)*0.1 + 5,

-                                showLogo = "bottom", showLogoScale = "right", logoColors = "rasmol",

-                                psFonts = TRUE, showLegend = TRUE, askForOverwrite = FALSE,

-                                furtherCode=c("\\defconsensus{.}{lower}{upper}","\\showruler{1}{top}"))

-            message(paste(file.name, "saved to", getwd()))

-        }

-}

new file mode 100644

@@ -0,0 +1,74 @@

+#' Clonality

+#' 

+#' Creates a data frame giving the total number of sequences, number of unique 

+#' productive sequences, number of genomes, entropy, clonality, Gini 

+#' coefficient, and the frequency (\%) of the top productive sequences in a list 

+#' of sample data frames.

+#' 

+#' @param file.list A list of data frames consisting of antigen receptor 

+#' sequencing imported by the LymphoSeq function readImmunoSeq. "aminoAcid", "count", 

+#' and "frequencyCount" are required columns.  "estimatedNumberGenomes" is optional.  

+#' Note that clonality is usually calculated from productive nucleotide sequences.  

+#' Therefore, it is not recommended to run this function using a productive sequence

+#' list aggregated by amino acids.

+#' @return Returns a data frame giving the total number of sequences, number of 

+#' unique productive sequences, number of genomes, clonality, Gini coefficient, 

+#' and the frequency (\%) of the top productive sequence in each sample.

+#' @details Clonality is derived from the Shannon entropy, which is calculated 

+#' from the frequencies of all productive sequences divided by the logarithm of 

+#' the total number of unique productive sequences.  This normalized entropy 

+#' value is then inverted (1 - normalized entropy) to produce the clonality 

+#' metric.  

+#' 

+#' The Gini coefficient is an alternative metric used to calculate repertoire 

+#' diversity and is derived from the Lorenz curve.  The Lorenz curve is drawn 

+#' such that x-axis represents the cumulative percentage of unique sequences and 

+#' the y-axis represents the cumulative percentage of reads.  A line passing 

+#' through the origin with a slope of 1 reflects equal frequencies of all clones.  

+#' The Gini coefficient is the ratio of the area between the line of equality 

+#' and the observed Lorenz curve over the total area under the line of equality.  

+#' Both Gini coefficient and clonality are reported on a scale from 0 to 1 where 

+#' 0 indicates all sequences have the same frequency and 1 indicates the 

+#' repertoire is dominated by a single sequence.

+#' @examples

+#' file.path <- system.file("extdata", "TCRB_sequencing", package = "LymphoSeq")

+#' 

+#' file.list <- readImmunoSeq(path = file.path)

+#' 

+#' clonality(file.list = file.list)

+#' @seealso \code{\link{lorenzCurve}}

+#' @export

+#' @importFrom ineq Gini

+clonality <- function(file.list) {

+    table <- data.frame(samples = names(file.list))

+    i <- 1

+    for (i in 1:length(file.list)) {

+        file <- file.list[[i]]

+        total.reads <- nrow(file)

+        file$estimatedNumberGenomes <- suppressWarnings(as.integer(file$estimatedNumberGenomes))

+        total.genomes <- sum(file$estimatedNumberGenomes)

+        total.count <- sum(file$count)

+        productive <- file[!grepl("\\*", file$aminoAcid) & file$aminoAcid != "", ]

+        frequency <- productive$count/sum(productive$count)

+        entropy <- -sum(frequency * log2(frequency), na.rm = TRUE)

+        unique.productive <- nrow(productive)

+        clonality <- 1 - round(entropy/log2(unique.productive), digits = 6)

+        table$totalSequences[i] <- total.reads

+        table$uniqueProductiveSequences[i] <- unique.productive

+        table$totalGenomes[i] <- total.genomes

+        table$totalCount[i] <- total.count

+        table$clonality[i] <- clonality

+        table$giniCoefficient[i] <- ineq::Gini(frequency)

+        table$topProductiveSequence[i] <- max(frequency) * 100

+        if (any(grepl("estimatedNumberGenomes", colnames(file)))) {

+            file$estimatedNumberGenomes <- suppressWarnings(as.integer(file$estimatedNumberGenomes))

+            total.genomes <- sum(file$estimatedNumberGenomes)

+            table$totalGenomes[i] <- ifelse(total.genomes == 0, NA, total.genomes)

+        } else {

+            table$totalGenomes[i] <- NA

+        }

+    }

+    table <- table[order(table$topProductiveSequence, decreasing = TRUE), ]

+    rownames(table) <- NULL

+    return(table)

+} 

\ No newline at end of file