Mercurial > repos > shians > shrnaseq

--- a/hairpinTool.R	Wed Oct 01 16:00:43 2014 +1000
+++ b/hairpinTool.R	Tue Oct 14 17:05:07 2014 +1100
@@ -1,45 +1,54 @@
 # ARGS: 1.inputType         -String specifying format of input (fastq or table)
-#    IF inputType is "fastQ":
+#    IF inputType is "fastq" or "pairedFastq:
 #       2*.fastqPath        -One or more strings specifying path to fastq files
-#       2.annoPath        -String specifying path to hairpin annotation table
+#       2.annoPath          -String specifying path to hairpin annotation table
 #       3.samplePath        -String specifying path to sample annotation table
 #       4.barStart          -Integer specifying starting position of barcode
 #       5.barEnd            -Integer specifying ending position of barcode
-#       6.hpStart           -Integer specifying startins position of hairpin
+#    ###
+#    IF inputType is "pairedFastq":
+#       6.barStartRev       -Integer specifying starting position of barcode
+#                            on reverse end
+#       7.barEndRev         -Integer specifying ending position of barcode
+#                            on reverse end
+#    ###
+#       8.hpStart           -Integer specifying startins position of hairpin
 #                            unique region
-#       7.hpEnd             -Integer specifying ending position of hairpin
+#       9.hpEnd             -Integer specifying ending position of hairpin
 #                            unique region
-#    ###
 #    IF inputType is "counts":
 #       2.countPath         -String specifying path to count table
 #       3.annoPath          -String specifying path to hairpin annotation table
 #       4.samplePath        -String specifying path to sample annotation table
 #    ###
-#       8.cpmReq            -Float specifying cpm requirement
-#       9.sampleReq         -Integer specifying cpm requirement
-#       10.fdrThresh        -Float specifying the FDR requirement
-#       11.lfcThresh        -Float specifying the log-fold-change requirement
-#       12.workMode         -String specifying exact test or GLM usage
-#       13.htmlPath         -String specifying path to HTML file
-#       14.folderPath       -STring specifying path to folder for output
+#       10.secFactName      -String specifying name of secondary factor
+#       11.cpmReq           -Float specifying cpm requirement
+#       12.sampleReq        -Integer specifying cpm requirement
+#       13.readReq          -Integer specifying read requirement
+#       14.fdrThresh        -Float specifying the FDR requirement
+#       15.lfcThresh        -Float specifying the log-fold-change requirement
+#       16.workMode         -String specifying exact test or GLM usage
+#       17.htmlPath         -String specifying path to HTML file
+#       18.folderPath       -String specifying path to folder for output
 #    IF workMode is "classic" (exact test)
-#       15.pairData[2]      -String specifying first group for exact test
-#       16.pairData[1]      -String specifying second group for exact test
+#       19.pairData[2]      -String specifying first group for exact test
+#       20.pairData[1]      -String specifying second group for exact test
 #    ###
 #    IF workMode is "glm"
-#       15.contrastData     -String specifying contrasts to be made
-#       16.roastOpt         -String specifying usage of gene-wise tests
-#       17.hairpinReq       -String specifying hairpin requirement for gene-
+#       19.contrastData     -String specifying contrasts to be made
+#       20.roastOpt         -String specifying usage of gene-wise tests
+#       21.hairpinReq       -String specifying hairpin requirement for gene-
 #                            wise test
-#       18.selectOpt        -String specifying type of selection for barcode
+#       22.selectOpt        -String specifying type of selection for barcode
 #                            plots
-#       19.selectVals       -String specifying members selected for barcode
+#       23.selectVals       -String specifying members selected for barcode
 #                            plots
 #    ###
 #
 # OUT:  Bar Plot of Counts Per Index
 #       Bar Plot of Counts Per Hairpin
 #       MDS Plot
+#       BCV Plot
 #       Smear Plot
 #       Barcode Plots (If Genewise testing was selected)
 #       Top Expression Table
@@ -58,12 +67,12 @@
 library(edgeR, quietly=TRUE, warn.conflicts=FALSE)
 library(limma, quietly=TRUE, warn.conflicts=FALSE)

-if (packageVersion("edgeR") < "3.5.27") {
-  stop("Please update 'edgeR' to version >= 3.5.23 to run this tool")
+if (packageVersion("edgeR") < "3.7.17") {
+  stop("Please update 'edgeR' to version >= 3.7.17 to run this tool")
 }

-if (packageVersion("limma")<"3.19.19") {
-  message("Update 'limma' to version >= 3.19.19 to see updated barcode graphs")
+if (packageVersion("limma")<"3.21.16") {
+  message("Update 'limma' to version >= 3.21.16 to see updated barcode graphs")
 }

 ################################################################################
@@ -171,103 +180,153 @@
 # Grabbing arguments from command line
 argv <- commandArgs(TRUE)

-# Remove fastq file paths after collecting from argument vector
 inputType <- as.character(argv[1])
-if (inputType=="fastq") {
+if (inputType == "fastq") {
+
   fastqPath <- as.character(gsub("fastq::", "", argv[grepl("fastq::", argv)],
                                  fixed=TRUE))
-  argv <- argv[!grepl("fastq::", argv, fixed=TRUE)]
+
+  # Remove fastq paths
+  argv <- argv[!grepl("fastq::", argv, fixed=TRUE)]
+
+  fastqPathRev <- NULL
   annoPath <- as.character(argv[2])
   samplePath <- as.character(argv[3])
   barStart <- as.numeric(argv[4])
   barEnd <- as.numeric(argv[5])
-  hpStart <- as.numeric(argv[6])
-  hpEnd <- as.numeric(argv[7])
-} else if (inputType=="counts") {
+  barStartRev <- NULL
+  barStartRev <- NULL
+  hpStart <- as.numeric(argv[8])
+  hpEnd <- as.numeric(argv[9])
+} else if (inputType=="pairedFastq") {
+
+  fastqPath <- as.character(gsub("fastq::", "", argv[grepl("fastq::", argv)],
+                                 fixed=TRUE))
+
+  fastqPathRev <- as.character(gsub("fastqRev::", "",
+                               argv[grepl("fastqRev::", argv)], fixed=TRUE))
+
+  # Remove fastq paths
+  argv <- argv[!grepl("fastq::", argv, fixed=TRUE)]
+  argv <- argv[!grepl("fastqRev::", argv, fixed=TRUE)]
+
+  annoPath <- as.character(argv[2])
+  samplePath <- as.character(argv[3])
+  barStart <- as.numeric(argv[4])
+  barEnd <- as.numeric(argv[5])
+  barStartRev <- as.numeric(argv[6])
+  barEndRev <- as.numeric(argv[7])
+  hpStart <- as.numeric(argv[8])
+  hpEnd <- as.numeric(argv[9])
+} else if (inputType == "counts") {
   countPath <- as.character(argv[2])
   annoPath <- as.character(argv[3])
   samplePath <- as.character(argv[4])
 }
-
-cpmReq <- as.numeric(argv[8])
-sampleReq <- as.numeric(argv[9])
-fdrThresh <- as.numeric(argv[10])
-lfcThresh <- as.numeric(argv[11])
-workMode <- as.character(argv[12])
-htmlPath <- as.character(argv[13])
-folderPath <- as.character(argv[14])

-if (workMode=="classic") {
+secFactName <- as.character(argv[10])
+cpmReq <- as.numeric(argv[11])
+sampleReq <- as.numeric(argv[12])
+readReq <- as.numeric(argv[13])
+fdrThresh <- as.numeric(argv[14])
+lfcThresh <- as.numeric(argv[15])
+selectDirection <- as.character(argv[16])
+workMode <- as.character(argv[17])
+htmlPath <- as.character(argv[18])
+folderPath <- as.character(argv[19])
+
+if (workMode == "classic") {
   pairData <- character()
-  pairData[2] <- as.character(argv[15])
-  pairData[1] <- as.character(argv[16])
-} else if (workMode=="glm") {
-  contrastData <- as.character(argv[15])
-  roastOpt <- as.character(argv[16])
-  hairpinReq <- as.numeric(argv[17])
-  selectOpt <- as.character(argv[18])
-  selectVals <- as.character(argv[19])
+  pairData[2] <- as.character(argv[20])
+  pairData[1] <- as.character(argv[21])
+} else if (workMode == "glm") {
+  contrastData <- as.character(argv[20])
+  roastOpt <- as.character(argv[21])
+  hairpinReq <- as.numeric(argv[22])
+  selectOpt <- as.character(argv[23])
+  selectVals <- as.character(argv[24])
 }

 # Read in inputs

 samples <- read.table(samplePath, header=TRUE, sep="\t")
+
 anno <- read.table(annoPath, header=TRUE, sep="\t")
-if (inputType=="counts") {
+
+if (inputType == "counts") {
   counts <- read.table(countPath, header=TRUE, sep="\t")
 }

 ###################### Check inputs for correctness ############################
 samples$ID <- make.names(samples$ID)

-if (!any(grepl("group", names(samples)))) {
+if ( !any(grepl("group", names(samples))) ) {
   stop("'group' column not specified in sample annotation file")
 } # Check if grouping variable has been specified

-if (any(table(samples$ID)>1)){
+if (secFactName != "none") {
+  if ( !any(grepl(secFactName, names(samples))) ) {
+  tempStr <- paste0("Second factor specified as \"", secFactName, "\" but ",
+                    "no such column name found in sample annotation file")
+  stop(tempStr)
+  } # Check if specified secondary factor is present
+}
+
+
+if ( any(table(samples$ID) > 1) ){
   tab <- table(samples$ID)
-  offenders <- paste(names(tab[tab>1]), collapse=", ")
+  offenders <- paste(names(tab[tab > 1]), collapse=", ")
   offenders <- unmake.names(offenders)
   stop("'ID' column of sample annotation must have unique values, values ",
        offenders, " are repeated")
 } # Check that IDs in sample annotation are unique

-if (inputType=="fastq") {
+if (inputType == "fastq" || inputType == "pairedFastq") {

-  if (any(table(anno$ID)>1)){
+  if ( any(table(anno$ID) > 1) ){
     tab <- table(anno$ID)
     offenders <- paste(names(tab[tab>1]), collapse=", ")
     stop("'ID' column of hairpin annotation must have unique values, values ",
     offenders, " are repeated")
   } # Check that IDs in hairpin annotation are unique

-} else if (inputType=="counts") {
-  if (any(is.na(match(samples$ID, colnames(counts))))) {
+} else if (inputType == "counts") {
+  # The first element of the colnames will be 'ID' and should not match
+  idFromSample <- samples$ID
+  idFromTable <- colnames(counts)[-1]
+  if (any(is.na(match(idFromTable, idFromSample)))) {
     stop("not all samples have groups specified")
   } # Check that a group has be specifed for each sample

-  if (any(table(counts$ID)>1)){
+  if ( any(table(counts$ID) > 1) ){
     tab <- table(counts$ID)
     offenders <- paste(names(tab[tab>1]), collapse=", ")
     stop("'ID' column of count table must have unique values, values ",
     offenders, " are repeated")
   } # Check that IDs in count table are unique
 }
-if (workMode=="glm") {
+if (workMode == "glm") {
   if (roastOpt == "yes") {
     if (is.na(match("Gene", colnames(anno)))) {
       tempStr <- paste("Gene-wise tests selected but'Gene' column not",
-      "specified in hairpin annotation file")
+                       "specified in hairpin annotation file")
       stop(tempStr)
     }
   }
 }

+if (secFactName != "none") {
+  if (workMode != "glm") {
+    tempStr <- paste("only glm analysis type possible when secondary factor",
+                     "used, please change appropriate option.")
+  }
+}
+
 ################################################################################

 # Process arguments
-if (workMode=="glm") {
-  if (roastOpt=="yes") {
+if (workMode == "glm") {
+  if (roastOpt == "yes") {
     wantRoast <- TRUE
   } else {
     wantRoast <- FALSE
@@ -299,7 +358,7 @@
   smearPng <- makeOut(paste0("smear(", pairData[2], "-", pairData[1],").png"))
   smearPdf <- makeOut(paste0("smear(", pairData[2], "-", pairData[1],").pdf"))
   topOut <- makeOut(paste0("toptag(", pairData[2], "-", pairData[1],").tsv"))
-} else if (workMode=="glm") {
+} else if (workMode == "glm") {
   smearPng <- character()
   smearPdf <- character()
   topOut <- character()
@@ -335,8 +394,8 @@
 ################################################################################

 # Transform gene selection from string into index values for mroast
-if (workMode=="glm") {
-  if (selectOpt=="rank") {
+if (workMode == "glm") {
+  if (selectOpt == "rank") {
     selectVals <- gsub(" ", "", selectVals, fixed=TRUE)
     selectVals <- unlist(strsplit(selectVals, ","))

@@ -356,33 +415,45 @@
   }
 }

-if (inputType=="fastq") {
+if (inputType == "fastq" || inputType == "pairedFastq") {
   # Use EdgeR hairpin process and capture outputs
   hpReadout <- capture.output(
-  data <- processHairpinReads(readfile=fastqPath,
-                              barcodefile=samplePath,
-                              hairpinfile=annoPath,
-                              barcodeStart=barStart, barcodeEnd=barEnd,
-                              hairpinStart=hpStart, hairpinEnd=hpEnd,
-                              verbose=TRUE)
+  data <- processAmplicons(readfile=fastqPath, readfile2=fastqPathRev,
+                            barcodefile=samplePath,
+                            hairpinfile=annoPath,
+                            barcodeStart=barStart, barcodeEnd=barEnd,
+                            barcodeStartRev=barStartRev,
+                            barcodeEndRev=barEndRev,
+                            hairpinStart=hpStart, hairpinEnd=hpEnd,
+                            verbose=TRUE)
   )
-
+
   # Remove function output entries that show processing data or is empty
   hpReadout <- hpReadout[hpReadout!=""]
   hpReadout <- hpReadout[!grepl("Processing", hpReadout)]
   hpReadout <- hpReadout[!grepl("in file", hpReadout)]
   hpReadout <- gsub(" -- ", "", hpReadout, fixed=TRUE)
-
+
   # Make the names of groups syntactically valid (replace spaces with periods)
   data$samples$group <- make.names(data$samples$group)
-} else if (inputType=="counts") {
+  if (secFactName != "none") {
+    data$samples[[secFactName]] <- make.names(data$samples[[secFactName]])
+  }
+} else if (inputType == "counts") {
   # Process counts information, set ID column to be row names
   rownames(counts) <- counts$ID
-  counts <- counts[ , !(colnames(counts)=="ID")]
+  counts <- counts[ , !(colnames(counts) == "ID")]
   countsRows <- nrow(counts)

   # Process group information
-  factors <- samples$group[match(samples$ID, colnames(counts))]
+  sampleNames <- colnames(counts)
+  matchedIndex <- match(sampleNames, samples$ID)
+  factors <- samples$group[matchedIndex]
+
+  if (secFactName != "none") {
+    secFactors <- samples[[secFactName]][matchedIndex]
+  }
+
   annoRows <- nrow(anno)
   anno <- anno[match(rownames(counts), anno$ID), ]
   annoMatched <- sum(!is.na(anno$ID))
@@ -416,14 +487,48 @@

 # Filter hairpins with low counts
 preFilterCount <- nrow(data)
-sel <- rowSums(cpm(data$counts) > cpmReq) >= sampleReq
-data <- data[sel, ]
+selRow <- rowSums(cpm(data$counts) > cpmReq) >= sampleReq
+selCol <- colSums(data$counts) > readReq
+data <- data[selRow, selCol]
+
+# Check if any data survived filtering
+if (length(data$counts) == 0) {
+  stop("no data remains after filtering, consider relaxing filters")
+}
+
+# Count number of filtered tags and samples
 postFilterCount <- nrow(data)
-filteredCount <- preFilterCount-postFilterCount
+filteredCount <- preFilterCount - postFilterCount
+sampleFilterCount <- sum(!selCol)
+
+if (secFactName == "none") {
+  # Estimate dispersions
+  data <- estimateDisp(data)
+  commonBCV <- round(sqrt(data$common.dispersion), 4)
+} else {
+  # Construct design
+  if (inputType == "counts") {
+
+    sampleNames <- colnames(counts)
+    matchedIndex <- match(sampleNames, samples$ID)
+    factors <- factor(make.names(samples$group[matchedIndex]))

-# Estimate dispersions
-data <- estimateDisp(data)
-commonBCV <- sqrt(data$common.dispersion)
+    secFactors <- factor(make.names(samples[[secFactName]][matchedIndex]))
+
+  } else if (inputType == "fastq" || inputType == "pairedFastq") {
+
+    factors <- factor(data$sample$group)
+    secFactors <- factor(data$sample[[secFactName]])
+
+  }
+
+  design <- model.matrix(~0 + factors + secFactors)
+
+  # Estimate dispersions
+  data <- estimateDisp(data, design=design)
+  commonBCV <- round(sqrt(data$common.dispersion), 4)
+}
+

 ################################################################################
 ### Output Processing
@@ -494,14 +599,23 @@
 linkData <- rbind(linkData, newEntry)
 invisible(dev.off())

-if (workMode=="classic") {
+if (workMode == "classic") {
   # Assess differential representation using classic exact testing methodology
   # in edgeR
   testData <- exactTest(data, pair=pairData)

   top <- topTags(testData, n=Inf)
+
+  if (selectDirection == "all") {
+    topIDs <- top$table[(top$table$FDR < fdrThresh) &
+                      (abs(top$table$logFC) > lfcThresh), 1]
+  } else if (selectDirection == "up") {
+    topIDs <- top$table[(top$table$FDR < fdrThresh) &
+                      (top$table$logFC > lfcThresh), 1]
+  } else if (selectDirection == "down") {
   topIDs <- top$table[(top$table$FDR < fdrThresh) &
-                      (abs(top$table$logFC) > lfcThresh), 1]
+                      (top$table$logFC < -lfcThresh), 1]
+}

   write.table(top, file=topOut, row.names=FALSE, sep="\t")

@@ -511,8 +625,10 @@
   linkData <- rbind(linkData, c(linkName, linkAddr))

   upCount[1] <- sum(top$table$FDR < fdrThresh & top$table$logFC > lfcThresh)
+
   downCount[1] <- sum(top$table$FDR < fdrThresh &
                       top$table$logFC < -lfcThresh)
+
   flatCount[1] <- sum(top$table$FDR > fdrThresh |
                       abs(top$table$logFC) < lfcThresh)

@@ -543,12 +659,40 @@
   linkData <- rbind(linkData, c(imgName, imgAddr))
   invisible(dev.off())

-} else if (workMode=="glm") {
+} else if (workMode == "glm") {
   # Generating design information
-  factors <- factor(data$sample$group)
-  design <- model.matrix(~0+factors)
+  if (secFactName == "none") {
+
+    factors <- factor(data$sample$group)
+    design <- model.matrix(~0 + factors)
+
+    colnames(design) <- gsub("factors", "", colnames(design), fixed=TRUE)
+
+  } else {
+
+    factors <- factor(data$sample$group)
+
+    if (inputType == "counts") {
+
+      sampleNames <- colnames(counts)
+      matchedIndex <- match(sampleNames, samples$ID)
+      factors <- factor(samples$group[matchedIndex])
+
+      secFactors <- factor(samples[[secFactName]][matchedIndex])
+
+    } else if (inputType == "fastq" || inputType == "pairedFastq") {
+
+      secFactors <- factor(data$sample[[secFactName]])
+
+    }
+
+    design <- model.matrix(~0 + factors + secFactors)
+
+    colnames(design) <- gsub("factors", "", colnames(design), fixed=TRUE)
+    colnames(design) <- gsub("secFactors", secFactName, colnames(design),
+                              fixed=TRUE)
+  }

-  colnames(design) <- gsub("factors", "", colnames(design), fixed=TRUE)

   # Split up contrasts seperated by comma into a vector
   contrastData <- unlist(strsplit(contrastData, split=","))
@@ -564,8 +708,18 @@

     # Select hairpins with FDR < 0.05 to highlight on plot
     top <- topTags(testData, n=Inf)
-    topIDs <- top$table[(top$table$FDR < fdrThresh) &
+
+    if (selectDirection == "all") {
+      topIDs <- top$table[(top$table$FDR < fdrThresh) &
                         (abs(top$table$logFC) > lfcThresh), 1]
+    } else if (selectDirection == "up") {
+      topIDs <- top$table[(top$table$FDR < fdrThresh) &
+                        (top$table$logFC > lfcThresh), 1]
+    } else if (selectDirection == "down") {
+      topIDs <- top$table[(top$table$FDR < fdrThresh) &
+                        (top$table$logFC < -lfcThresh), 1]
+    }
+
     write.table(top, file=topOut[i], row.names=FALSE, sep="\t")

     linkName <- paste0("Top Tags Table(", contrastData[i], ") (.tsv)")
@@ -607,8 +761,8 @@
     unq <- unq[!is.na(unq)]
     geneList <- list()
     for (gene in unq) {
-      if (length(which(genes==gene)) >= hairpinReq) {
-        geneList[[gene]] <- which(genes==gene)
+      if (length(which(genes == gene)) >= hairpinReq) {
+        geneList[[gene]] <- which(genes == gene)
       }
     }

@@ -624,7 +778,7 @@
                          ") (.tsv)")
       linkAddr <- paste0("gene_level(", contrastData[i], ").tsv")
       linkData <- rbind(linkData, c(linkName, linkAddr))
-      if (selectOpt=="rank") {
+      if (selectOpt == "rank") {
         selectedGenes <- rownames(roastData)[selectVals]
       } else {
         selectedGenes <- selectVals
@@ -668,9 +822,13 @@
 # Generate data frame of the significant differences
 sigDiff <- data.frame(Up=upCount, Flat=flatCount, Down=downCount)
 if (workMode == "glm") {
+
   row.names(sigDiff) <- contrastData
+
 } else if (workMode == "classic") {
+
   row.names(sigDiff) <- paste0(pairData[2], "-", pairData[1])
+
 }

 # Output table of summarised counts
@@ -702,7 +860,8 @@
 cata("<body>\n")
 cata("<h3>EdgeR Analysis Output:</h3>\n")
 cata("<h4>Input Summary:</h4>\n")
-if (inputType=="fastq") {
+if (inputType == "fastq" || inputType == "pairedFastq") {
+
   cata("<ul>\n")
   ListItem(hpReadout[1])
   ListItem(hpReadout[2])
@@ -716,31 +875,51 @@
   cata("<br />\n")
   cata("<b>Please check that read percentages are consistent with ")
   cata("expectations.</b><br >\n")
-} else if (inputType=="counts") {
+
+} else if (inputType == "counts") {
+
   cata("<ul>\n")
   ListItem("Number of Samples: ", ncol(data$counts))
   ListItem("Number of Hairpins: ", countsRows)
   ListItem("Number of annotations provided: ", annoRows)
   ListItem("Number of annotations matched to hairpin: ", annoMatched)
   cata("</ul>\n")
+
 }

 cata("The estimated common biological coefficient of variation (BCV) is: ",
      commonBCV, "<br />\n")

+if (secFactName == "none") {
+
+  cata("No secondary factor specified.<br />\n")
+
+} else {
+
+  cata("Secondary factor specified as: ", secFactName, "<br />\n")
+
+}
+
 cata("<h4>Output:</h4>\n")
 cata("PDF copies of JPEGS available in 'Plots' section.<br />\n")
 for (i in 1:nrow(imageData)) {
   if (grepl("barcode", imageData$Link[i])) {
+
     if (packageVersion("limma")<"3.19.19") {
+
       HtmlImage(imageData$Link[i], imageData$Label[i],
                 height=length(selectedGenes)*150)
+
     } else {
+
       HtmlImage(imageData$Link[i], imageData$Label[i],
                 height=length(selectedGenes)*300)
+
     }
   } else {
+
     HtmlImage(imageData$Link[i], imageData$Label[i])
+
   }
 }
 cata("<br />\n")
@@ -779,26 +958,42 @@
 }

 cata("<p>Alt-click links to download file.</p>\n")
-cata("<p>Click floppy disc icon associated history item to download ")
+cata("<p>Click floppy disc icon on associated history item to download ")
 cata("all files.</p>\n")
 cata("<p>.tsv files can be viewed in Excel or any spreadsheet program.</p>\n")

 cata("<h4>Additional Information:</h4>\n")

 if (inputType == "fastq") {
+
   ListItem("Data was gathered from fastq raw read file(s).")
+
 } else if (inputType == "counts") {
+
   ListItem("Data was gathered from a table of counts.")
+
 }

-if (cpmReq!=0 && sampleReq!=0) {
-  tempStr <- paste("Hairpins without more than", cpmReq,
+if (cpmReq != 0 && sampleReq != 0) {
+  tempStr <- paste("Target sequences without more than", cpmReq,
                    "CPM in at least", sampleReq, "samples are insignificant",
                    "and filtered out.")
   ListItem(tempStr)
+
   filterProp <- round(filteredCount/preFilterCount*100, digits=2)
   tempStr <- paste0(filteredCount, " of ", preFilterCount," (", filterProp,
-                   "%) hairpins were filtered out for low count-per-million.")
+                   "%) target sequences were filtered out for low ",
+                   "count-per-million.")
+  ListItem(tempStr)
+}
+
+if (sampleReq != 0) {
+  tempStr <- paste("Samples that did not produce more than", sampleReq,
+                   "counts were filtered out.")
+  ListItem(tempStr)
+
+  tempStr <- paste0(sampleFilterCount, " samples were filtered out for low ",
+                    "counts.")
   ListItem(tempStr)
 }

@@ -809,9 +1004,9 @@
 }

 if (workMode == "classic") {
-  ListItem("An exact test was performed on each hairpin.")
+  ListItem("An exact test was performed on each target sequence.")
 } else if (workMode == "glm") {
-  ListItem("A generalised linear model was fitted to each hairpin.")
+  ListItem("A generalised linear model was fitted to each target sequence.")
 }

 cit <- character()
--- a/hairpinTool.xml	Wed Oct 01 16:00:43 2014 +1000
+++ b/hairpinTool.xml	Tue Oct 14 17:05:07 2014 +1100
@@ -1,12 +1,13 @@
-<tool id="shRNAseq" name="shRNAseq Tool" version="1.0.13">
+<tool id="shRNAseq" name="shRNAseq Tool" version="1.2.0">
   <description>
-    Analyse hairpin differential representation using edgeR
+    Analyse differential representation for shRNAseq and sgRNA based procedures
+    using edgeR package from Bioconductor.
   </description>

   <requirements>
-    <requirement type="R-module" version="3.6.2">edgeR</requirement>
-    <requirement type="R-module" version="3.20.7">limma</requirement>
-    <requirement type="package" version="3.0.3">R_3_0_3</requirement>
+    <requirement type="R-module" version="3.7.17">edgeR</requirement>
+    <requirement type="R-module" version="3.21.16">limma</requirement>
+    <requirement type="package" version="3.1.1">R_3_0_3</requirement>
   </requirements>

   <stdio>
@@ -14,43 +15,90 @@
   </stdio>

   <command interpreter="Rscript">
-  hairpinTool.R $inputOpt.inputType
+  ampliconTool.R $inputOpt.inputType
                 #if $inputOpt.inputType=="fastq":
+
                   #for $i, $fas in enumerate($inputOpt.fastq):
                     fastq::$fas.file
                   #end for

                   $inputOpt.hairpin
                   $inputOpt.samples
+
+                  #if $inputOpt.positions.posOption=="yes":
+                    $inputOpt.positions.barstart
+                    $inputOpt.positions.barend
+                    0
+                    0
+                    $inputOpt.positions.hpstart
+                    $inputOpt.positions.hpend
+                  #else:
+                    1
+                    5
+                    0
+                    0
+                    37
+                    57
+                  #end if
+                #elif $inputOpt.inputType=="pairedFastq":
+
+                  #for $i, $fas in enumerate($inputOpt.fastq):
+                    fastq::$fas.file
+                  #end for
+
+                  #for $i, $fas in enumerate($inputOpt.fastq):
+                    fastqRev::$fas.fileRev
+                  #end for
+
+                  $inputOpt.hairpin
+                  $inputOpt.samples

                   #if $inputOpt.positions.posOption=="yes":
                     $inputOpt.positions.barstart
                     $inputOpt.positions.barend
+                    $inputOpt.positions.barstartRev
+                    $inputOpt.positions.barendRev
                     $inputOpt.positions.hpstart
                     $inputOpt.positions.hpend
                   #else:
                     1
                     5
+                    0
+                    0
                     37
                     57
                   #end if
-                #else:
+
+                #elif $inputOpt.inputType=="counts":
                   $inputOpt.counts
                   $inputOpt.hairpin
                   $inputOpt.samples
-                  0 0 0
+                  0
+                  0
+                  0
+                  0
+                  0
                 #end if
-
+
+                #if $inputOpt.secondaryFactor.secFactorOpt=="yes":
+                  $inputOpt.secondaryFactor.secFactName
+                #else:
+                  "none"
+                #end if
+
                 #if $filterCPM.filtOption=="yes":
                   $filterCPM.cpmReq
                   $filterCPM.sampleReq
+                  $filterCPM.readReq
                 #else:
                   -Inf
                   -Inf
+                  -Inf
                 #end if

                 $fdr
                 $lfc
+                $direction
                 $workMode.mode
                 $outFile
                 $outFile.files_path
@@ -61,6 +109,7 @@
                 #elif $workMode.mode=="glm":
                   "$workMode.contrast"
                   $workMode.roast.roastOption
+
                   #if $workMode.roast.roastOption=="yes":
                     $workMode.roast.hairpinReq
                     $workMode.roast.select.selOption
@@ -70,20 +119,22 @@
                     0
                     0
                   #end if
+
                 #end if
   </command>

   <inputs>
     <conditional name="inputOpt">
+
       <param name="inputType" type="select" label="Input File Type">
         <option value="fastq">FastQ File</option>
+        <option value="pairedFastq">Paired FastQ File</option>
         <option value="counts">Table of Counts</option>
       </param>
-
+
       <when value="fastq">
         <param name="hairpin" type="data" format="tabular"
-               label="Hairpin Annotation"/>
-
+               label="Target Annotation"/>

         <param name="samples" type="data" format="tabular"
                label="Sample Annotation"/>
@@ -91,47 +142,171 @@
         <repeat name="fastq" title="FastQ Files">
           <param name="file" type="data" format="fastq"/>
         </repeat>
+
+        <conditional name="secondaryFactor">

+          <param name="secFactorOpt" type="select"
+                 label="Include Secondary Factor">
+
+            <option value="no" selected="True">No</option>
+
+            <option value="yes">Yes</option>
+
+          </param>
+
+          <when value="yes">
+
+            <param name="secFactName" type="text" label="Secondary Factor Name"
+                   size="80"/>
+
+          </when>
+
+          <when value="no">
+          </when>
+        </conditional>
+
         <conditional name="positions">
           <param name="posOption" type="select"
-                 label="Specify Barcode and Hairpin Locations?"
-                 help="Default Positions: Barcode: 1 to 5, Hairpin: 37 to 57.">
+                 label="Specify Sample Index and Target Sequence Locations?"
+                 help="Default Positions: Index: 1 to 5, Target: 37 to 57.">
             <option value="no" selected="True">No</option>
             <option value="yes">Yes</option>
           </param>

           <when value="yes">
             <param name="barstart" type="integer" value="1"
-                   label="Barcode Starting Position"/>
+                   label="Index Starting Position"/>
             <param name="barend" type="integer" value="5"
-                   label="Barcode Ending Position"/>
+                   label="Index Ending Position"/>

             <param name="hpstart" type="integer" value="37"
-                   label="Hairpin Starting Position"/>
+                   label="Target Starting Position"/>

             <param name="hpend" type="integer" value="57"
-                   label="Hairpin Ending Position"/>
+                   label="Target Ending Position"/>
           </when>

           <when value="no"/>
         </conditional>
       </when>
-
+
+      <when value="pairedFastq">
+        <param name="hairpin" type="data" format="tabular"
+               label="Target Sequence Annotation"/>
+
+        <param name="samples" type="data" format="tabular"
+               label="Sample Annotation"/>
+
+        <repeat name="fastq" title="FastQ Files">
+          <param name="file" type="data" format="fastq"/>
+          <param name="fileRev" type="data" format="fastq"/>
+        </repeat>
+
+        <conditional name="secondaryFactor">
+
+          <param name="secFactorOpt" type="select"
+                 label="Include Secondary Factor">
+
+            <option value="no" selected="True">No</option>
+
+            <option value="yes">Yes</option>
+
+          </param>
+
+          <when value="yes">
+
+            <param name="secFactName" type="text" label="Secondary Factor Name"
+                   size="80"/>
+
+          </when>
+
+          <when value="no">
+          </when>
+        </conditional>
+
+        <conditional name="positions">
+
+          <param name="posOption" type="select"
+                 label="Specify Sample Index and Target Sequence Locations?"
+                 help="Default Positions: Index: 1 to 5, Input required for
+                       reverse end, Target: 37 to 57.">
+
+            <option value="no" selected="True">No</option>
+
+            <option value="yes">Yes</option>
+
+          </param>
+
+          <when value="yes">
+            <param name="barstart" type="integer" value="1"
+                   label="Index Starting Position"/>
+
+            <param name="barend" type="integer" value="5"
+                   label="Index Ending Position"/>
+
+            <param name="barstartRev" type="integer" value="0"
+                   label="Reverse Index Starting Position"/>
+
+            <param name="barendRev" type="integer" value="0"
+                   label="Reverse Index Ending Position"/>
+
+            <param name="hpstart" type="integer" value="37"
+                   label="Target Starting Position"/>
+
+            <param name="hpend" type="integer" value="57"
+                   label="Target Ending Position"/>
+          </when>
+
+          <when value="no">
+          </when>
+
+        </conditional>
+
+      </when>
+
       <when value="counts">
+
         <param name="counts" type="data" format="tabular" label="Counts Table"/>
+
         <param name="hairpin" type="data" format="tabular"
-               label="Hairpin Annotation"/>
+               label="Target Sequence Annotation"/>
+
         <param name="samples" type="data" format="tabular"
                label="Sample Annotation"/>
+
+        <conditional name="secondaryFactor">
+
+          <param name="secFactorOpt" type="select"
+                 label="Include Secondary Factor">
+
+            <option value="no" selected="True">No</option>
+
+            <option value="yes">Yes</option>
+
+          </param>
+
+          <when value="yes">
+
+            <param name="secFactName" type="text" label="Secondary Factor Name"
+                   size="80"/>
+
+          </when>
+
+          <when value="no">
+          </when>
+
+        </conditional>
+
       </when>
+
     </conditional>

     <conditional name="filterCPM">
       <param name="filtOption" type="select" label="Filter Low CPM?"
-       help="Ignore hairpins with very low representation when performing
-             analysis.">
+       help="Ignore target sequences with very low representation when
+             performing analysis.">
         <option value="yes">Yes</option>
-       	<option value="no">No</option>
+        <option value="no">No</option>
       </param>

         <when value="yes">
@@ -142,6 +317,12 @@
                  label="Minimum Samples"
                  help="Filter out all the genes that do not meet the minimum
                        CPM in at least this many samples."/>
+
+          <param name="readReq" type="integer" value="1000" min="0"
+                 label="Minimum Reads"
+                 help="Filter out all samples that do not have the minimum
+                       number of reads."/>
+
         </when>

         <when value="no"/>
@@ -175,17 +356,18 @@
         <conditional name="roast">
           <param name="roastOption" type="select"
                  label="Perform Gene Level Analysis?"
-                 help="Analyse LogFC tendencies for hairpins belonging
-                       to the same gene.">
+                 help="Analyse LogFC tendencies for target sequences belonging
+                       to the same gene. NOTE: this is a slow procedure that
+                       scales badly with the number of genes analysed.">
             <option value="no">No</option>
             <option value="yes">Yes</option>
           </param>

           <when value="yes">
             <param name="hairpinReq" type="integer" value="2" min="2"
-                   label="Minimum Hairpins"
-                   help="Only genes with at least this many hairpins will
-                         be analysed."/>
+                   label="Minimum Targets Found"
+                   help="Only genes with at least this many target sequences
+                         found will be analysed."/>

             <conditional name="select">
               <param name="selOption" type="select"
@@ -223,25 +405,33 @@
            label="FDR Threshold"
            help="All observations below this threshold will be highlighted
                  in the smear plot."/>
+
     <param name="lfc" type="float" value="0" min="0"
            label="Absolute LogFC Threshold"
            help="In additional to meeting the FDR requirement, the absolute
                  value of the log-fold-change of the observation must be above
                  this threshold to be highlighted."/>
+
+    <param name="direction" type="select" label="Highlight Option"
+        help="Only hightlight positive or negative fold changes in smear plot?">
+        <option value="all">Default</option>
+        <option value="up">Positive Only</option>
+        <option value="down">Negative Only</option>
+    </param>
   </inputs>

   <outputs>
-    <data format="html" name="outFile" label="shRNAseq Analysis"/>
+    <data format="html" name="outFile" label="TagSeq Analysis"/>
   </outputs>
   <help>
 .. class:: infomark

 **What it does**

-Given tables containing information about the hairpins and their associated
-barcodes, information about the samples and fastq file containing the hairpin
-reads. This tool will generate plots and tables for the analysis of differential
-representation.
+Given tables containing information about the hairpins/sgRNA and their
+associated sample indices, information about the samples and fastq file
+containing the sequencing reads. This tool will generate plots and tables for
+the analysis of differential representation.

 .. class:: infomark

@@ -257,14 +447,15 @@
 **Input File Type:**

 This tool is able to either generate counts from a raw FastQ file given the
-information regarding the samples and hairpins. Alternatively if a table of
-counts has already been generated it can also be used.
+information regarding the samples and hairpins/sgRNA. Alternatively if a table
+of counts has already been generated it can also be used.

 **Counts Table (Counts Input):**

-A tab delimited text table of information regarding the counts of hairpins.
-Should have a column 'ID' to denote the hairpins that counts correspond to. Each
-additional column should have titles corresponding to the label for the sample.
+A tab delimited text table of information regarding the counts of
+hairpins/sgRNA. Should have a column 'ID' to denote the hairpins/sgRNA that
+counts correspond to. Each additional column should have titles corresponding to
+the label for the sample.

 Example::

@@ -281,62 +472,80 @@
   Hairpin7 49501 49076 47611
   ...

-**Hairpin Annotation:**
+**Target Sequence Annotation:**

-A tab delimited text table of information regarding the hairpins. Should have
-columns 'ID', 'Sequences' and 'Gene' to uniquely identify the hairpin, align it
-with the reads to produce counts and identify which gene the hairpin acts on.
+A tab delimited text table of information regarding the targetted
+hairpins/sgRNA sequence. Should have columns 'ID', 'Sequences' and 'Gene' to
+uniquely identify the target, align it with the reads to produce counts and
+identify which gene the target acts on.

 NOTE: the column names are case sensitive and should be input exactly as they
 are shown here.

 Example::

-  ID	Sequences	Gene
-  Control1	TCTCGCTTGGGCGAGAGTAAG	2
-  Control2	CCGCCTGAAGTCTCTGATTAA	2
-  Control3	AGGAATTATAATGCTTATCTA	2
-  Hairpin1	AAGGCAGAGACTGACCACCTA	4
-  Hairpin2	GAGCGACCTGGTGTTACTCTA	4
-  Hairpin3	ATGGTGTAAATAGAGCTGTTA	4
-  Hairpin4	CAGCTCATCTTCTGTGAAGAA	4
-  Hairpin5	CAGCTCTGTGGGTCAGAAGAA	4
-  Hairpin6	CCAGGCACAGATCTCAAGATA	4
-  Hairpin7	ATGACAAGAAAGACATCTCAA	7
+  ID  Sequences Gene
+  Control1  TCTCGCTTGGGCGAGAGTAAG 2
+  Control2  CCGCCTGAAGTCTCTGATTAA 2
+  Control3  AGGAATTATAATGCTTATCTA 2
+  Hairpin1  AAGGCAGAGACTGACCACCTA 4
+  Hairpin2  GAGCGACCTGGTGTTACTCTA 4
+  Hairpin3  ATGGTGTAAATAGAGCTGTTA 4
+  Hairpin4  CAGCTCATCTTCTGTGAAGAA 4
+  Hairpin5  CAGCTCTGTGGGTCAGAAGAA 4
+  Hairpin6  CCAGGCACAGATCTCAAGATA 4
+  Hairpin7  ATGACAAGAAAGACATCTCAA 7
   ...

 **Sample Annotation (FastQ Input):**

 A tab delimited text table of information regarding the samples. Should have
 columns 'ID', 'Sequences' and 'group' to uniquely identify each sample, identify
-the sample in the reads by its barcode sequence and correctly group replicates
-for analysis. Additional columns may inserted for annotation purposes and will
-not interfere with analysis as long as the necessary columns are present.
+the sample in the reads by its sample index sequence and correctly group
+replicates for analysis. Additional columns may inserted for annotation purposes
+and will not interfere with analysis as long as the necessary columns are
+present.

-NOTE: the column names are case sensitive and should be input exactly as they
-are shown here.
+NOTE: With the exception of other_group, column names are case sensitive and
+should be input exactly as they are shown here. The other_group column can be
+named by the user and specified in the "Include Secondary Factor" option of the
+tool.

 Example::

-  ID	Sequences	group	Replicate
-  3	GAAAG	Day 2	1
-  6	GAACC	Day 10	1
-  9	GAAGA	Day 5 GFP neg	1
-  16	GAATT	Day 5 GFP pos	1
-  18	GACAC	Day 2	2
-  21	GACCA	Day 10	2
-  28	GACGT	Day 5 GFP neg	2
-  31	GACTG	Day 5 GFP pos	2
-  33	GAGAA	Day 2	3
-  40	GAGCT	Day 10	3
+  ID  Sequences group other_group Replicate
+  3 GAAAG Day 2 male 1
+  6 GAACC Day 10  female  1
+  9 GAAGA Day 5 GFP neg male 1
+  16  GAATT Day 5 GFP pos male 1
+  18  GACAC Day 2 female 2
+  21  GACCA Day 10  male  2
+  28  GACGT Day 5 GFP neg male 2
+  31  GACTG Day 5 GFP pos female 2
+  33  GAGAA Day 2 male 3
+  40  GAGCT Day 10  female  3
   ...

-**Specify Barcode and Hairpin Locations (FastQ Input):**
+**Include Secondary Factor**
+If there are two factors involved in the experiment (i.e. Age and Gender) then
+then secondary factor should be included to improve the statistical analysis.
+The secondary factor should be specified as a column in the sample annotation
+file and the corresponding column name should be input exactly as it is into
+the provided field in the tool.
+
+NOTE: Currently the secondary factor is used only to improve statistical
+analysis, comparisons can only be made in the primary factor specified as
+"group" in the sample annotation.
+
+**Specify Sample Index and Target Sequence Locations (FastQ Input):**

 It is assumed that in the sequencing reads that the first 5 bases are the
-barcodes and that bases 37-57 are the hairpins. If this is not the case then the
-values of the positions can be changed, however it still requires the barcodes
-and hairpins to be in a consistent location an in a continuous sequence.
+sample index sequence and that bases 37-57 are the hairpins/sgRNA. If this is
+not the case then the values of the positions can be changed, however it still
+requires the sample indices and hairpins/sgRNA to be in a consistent location an
+in a continuous sequence.
+
+NOTE: position values start at 1 for the first base.

 **Filter Low CPM?:**

@@ -346,21 +555,21 @@

 **Analysis Type:**

- * **Classic Exact Test:** This allows two experimental groups to be compared and
-   p-values for differential representation derivec for each hairpin. Simple and
-   fast for straightforward comparisons. In this option you will have the option of
-   "*Compare* x *To* y" which implicitly subtracts the data from y from that of x
-   to produce the comparison.
+ * **Classic Exact Test:** This allows two experimental groups to be compared
+   and p-values for differential representation derivec for each target
+   sequence. Simple and fast for straightforward comparisons. In this option you
+   will have the option of "*Compare* x *To* y" which implicitly subtracts the
+   data from y from that of x to produce the comparison.

- * **Generalised Linear Model:** This allow for complex contrasts to be specified
-   and also gene level analysis to be performed. If this option is chosen then
-   contrasts must be explicitly stated in equations and multiple contrasts can be
-   made. In addition there will be the option to analyse hairpins on a per-gene
-   basis to see if hairpins belonging to a particular gene have any overall
-   tendencies for the direction of their log-fold-change.
+ * **Generalised Linear Model:** This allow for complex contrasts to be specified
+   and also gene level analysis to be performed. If this option is chosen then
+   contrasts must be explicitly stated in equations and multiple contrasts can
+   be made. In addition there will be the option to analyse hairpins/sgRNA on a
+   per-gene basis to see if hairpins/sgRNA belonging to a particular gene have
+   any overall tendencies for the direction of their log-fold-change.

 **FDR Threshold:**
-The smear plot in the output will have hairpins highlighted to signify
+The smear plot in the output will have hairpins/sgRNA highlighted to signify
 significant differential representation. The significance is determined by
 contorlling the false discovery rate, only those with a FDR lower than the
 threshold will be highlighted in the plot.
@@ -379,10 +588,10 @@
 using.  The methodology articles are listed in Section 2.1 of the limma
 User's Guide.

-	* Smyth, GK (2005). Limma: linear models for microarray data. In:
-	  'Bioinformatics and Computational Biology Solutions using R and
-	  Bioconductor'. R. Gentleman, V. Carey, S. Dudoit, R. Irizarry,
-	  W. Huber (eds), Springer, New York, pages 397-420.
+  * Smyth, GK (2005). Limma: linear models for microarray data. In:
+    'Bioinformatics and Computational Biology Solutions using R and
+    Bioconductor'. R. Gentleman, V. Carey, S. Dudoit, R. Irizarry,
+    W. Huber (eds), Springer, New York, pages 397-420.

 .. class:: infomark

@@ -392,25 +601,24 @@
 the various original statistical methods implemented in edgeR.  See
 Section 1.2 in the User's Guide for more detail.

-	* Robinson MD, McCarthy DJ and Smyth GK (2010). edgeR: a Bioconductor
-	  package for differential expression analysis of digital gene expression
-	  data. Bioinformatics 26, 139-140
-
-	* Robinson MD and Smyth GK (2007). Moderated statistical tests for assessing
-	  differences in tag abundance. Bioinformatics 23, 2881-2887
-
-	* Robinson MD and Smyth GK (2008). Small-sample estimation of negative
-	  binomial dispersion, with applications to SAGE data.
-	  Biostatistics, 9, 321-332
-
-	* McCarthy DJ, Chen Y and Smyth GK (2012). Differential expression analysis
-	  of multifactor RNA-Seq experiments with respect to biological variation.
-	  Nucleic Acids Research 40, 4288-4297
-
+  * Robinson MD, McCarthy DJ and Smyth GK (2010). edgeR: a Bioconductor
+    package for differential expression analysis of digital gene expression
+    data. Bioinformatics 26, 139-140
+
+  * Robinson MD and Smyth GK (2007). Moderated statistical tests for assessing
+    differences in tag abundance. Bioinformatics 23, 2881-2887
+
+  * Robinson MD and Smyth GK (2008). Small-sample estimation of negative
+    binomial dispersion, with applications to SAGE data.
+    Biostatistics, 9, 321-332
+
+  * McCarthy DJ, Chen Y and Smyth GK (2012). Differential expression analysis
+    of multifactor RNA-Seq experiments with respect to biological variation.
+    Nucleic Acids Research 40, 4288-4297
+
 Report problems to: su.s@wehi.edu.au

 .. _edgeR: http://www.bioconductor.org/packages/release/bioc/html/edgeR.html
 .. _limma: http://www.bioconductor.org/packages/release/bioc/html/limma.html
   </help>
 </tool>
-