comparison hairpinTool.R @ 6:3d04308a99f9

- Added differentially expressed hairpin count output - Added running time output - Added counts table output

5:17befe9f8b03

# Record starting time
timeStart <- as.character(Sys.time())

# Loading and checking required packages
library(methods, quietly=TRUE, warn.conflicts=FALSE)

inputType <- as.character(argv[1])
if (inputType=="fastq") {
  fastqPath <- as.character(gsub("fastq::", "", argv[grepl("fastq::", argv)], 
                                 fixed=TRUE))
  argv <- argv[!grepl("fastq::", argv, fixed=TRUE)]
  hairpinPath <- 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])

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") {
  pairData <- character()
  pairData[2] <- as.character(argv[15])
  pairData[1] <- as.character(argv[16])
} else if (workMode=="glm") {

  selectOpt <- as.character(argv[18])
  selectVals <- as.character(argv[19])
}

# Read in inputs
if (inputType=="fastq") {
  samples <- read.table(samplePath, header=TRUE, sep="\t")
  hairpins <- read.table(hairpinPath, header=TRUE, sep="\t")
} else if (inputType=="counts") {
  samples <- read.table(samplePath, header=TRUE, sep="\t")
  counts <- read.table(countPath, header=TRUE, sep="\t")
  anno <- read.table(annoPath, header=TRUE, sep="\t")
}
###################### Check inputs for correctness ############################
samples$ID <- make.names(samples$ID)

if (!any(grepl("group", names(samples)))) {
  stop("'group' column not specified in sample annotation file")

if (any(table(samples$ID)>1)){
  tab <- table(samples$ID)
  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 (any(table(hairpins$ID)>1)){
    tab <- table(hairpins$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))))) {

  } # Check that a group has be specifed for each sample
  
  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
}
################################################################################

# Process arguments
if (workMode=="glm") {
  if (roastOpt=="yes") {

if (exists("pairData")) {
  pairData <- make.names(pairData)
}

# Generate output folder and paths
dir.create(folderPath)

# Generate links for outputs
imgOut("barHairpin")
imgOut("barIndex")
imgOut("mds")

    roastOut[i] <- makeOut(paste0("roast(", contrastData[i], ").tsv"))
    barcodePng[i] <- makeOut(paste0("barcode(", contrastData[i], ").png"))
    barcodePdf[i] <- makeOut(paste0("barcode(", contrastData[i], ").pdf"))
  }
}
# Initialise data for html links and images, table with the link label and
# link address
linkData <- data.frame(Label=character(), Link=character(),
                       stringsAsFactors=FALSE)
imageData <- data.frame(Label=character(), Link=character(),

    selectVals <- unlist(strsplit(selectVals, " "))
  }                                                           
}
                                                  
if (inputType=="fastq") {
# Use EdgeR hairpin process and capture outputs
hpReadout <- capture.output(
  data <- processHairpinReads(fastqPath, samplePath, hairpinPath,
                              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") {
  # Process counts information, set ID column to be row names
  rownames(counts) <- counts$ID
  counts <- counts[ , !(colnames(counts)=="ID")]
  countsRows <- nrow(counts)

  anno <- anno[sel, ]
  
  # Create DGEList
  data <- DGEList(counts=counts, lib.size=colSums(counts), 
                  norm.factors=rep(1,ncol(counts)), genes=anno, group=factors)
                  
  # Make the names of groups syntactically valid (replace spaces with periods)
  data$samples$group <- make.names(data$samples$group)
}

# Filter hairpins with low counts
sel <- rowSums(cpm(data$counts) > cpmReq) >= sampleReq
data <- data[sel, ]

# Estimate dispersions
data <- estimateDisp(data)
commonBCV <- sqrt(data$common.dispersion)

      dev.off()
    }
  }
}

# Record ending time
timeEnd <- as.character(Sys.time())
################################################################################
### HTML Generation
################################################################################
# Clear file
cat("", file=htmlPath)

if (inputType=="fastq") {
  cata("<ul>\n")
  ListItem(hpReadout[1])
  ListItem(hpReadout[2])
  cata("</ul>\n")
  cata(hpReadout[3], "<br/>\n")
  cata("<ul>\n")
  ListItem(hpReadout[4])
  ListItem(hpReadout[7])
  cata("</ul>\n")
  cata(hpReadout[8:11], sep="<br/>\n")
  cata("<br />\n")
  cata("<b>Please check that read percentages are consistent with ")
  cata("expectations.</b><br >\n")
} else if (inputType=="counts") {
  cata("<ul>\n")

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

cata("<h4>Output:</h4>\n")
cata("All images displayed have PDF copy at the bottom of the page, these can ")
cata("exported in a pdf viewer to high resolution image format. <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])
  }
}
cata("<br/>\n")

cata("<h4>Plots:</h4>\n")
for (i in 1:nrow(linkData)) {
  if (!grepl(".tsv", linkData$Link[i])) {
    HtmlLink(linkData$Link[i], linkData$Label[i])

cata("<p>alt-click any of the links to download the file, or click the name ")
cata("of this task in the galaxy history panel and click on the floppy ")
cata("disk icon to download all files in a zip archive.</p>\n")
cata("<p>.tsv files are tab seperated files that can be viewed using Excel ")
cata("or other spreadsheet programs</p>\n")
cata("<table border=\"0\">\n")

cata("<tr>\n")
TableItem("Task started at:"); TableItem(timeStart)
cata("</tr>\n")
cata("<tr>\n")
TableItem("Task ended at:"); TableItem(timeEnd)
cata("</tr>\n")

cata("</body>\n")
cata("</html>")

6:3d04308a99f9

# ARGS: 1.inputType         -String specifying format of input (fastq or table)
#    IF inputType is "fastQ":
#       2*.fastqPath        -One or more strings specifying path to fastq files
#       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
#                            unique region
#       7.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
#    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
#    ###
#    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-
#                            wise test
#       18.selectOpt        -String specifying type of selection for barcode
#                            plots
#       19.selectVals       -String specifying members selected for barcode
#                            plots
#
# OUT:  Bar Plot of Counts Per Index
#       Bar Plot of Counts Per Hairpin
#       MDS Plot
#       Smear Plot
#       Barcode Plots (If Genewise testing was selected)
#       Top Expression Table
#       HTML file linking to the ouputs
#
# Author: Shian Su - [e-mail] - Jan 2014

# Record starting time
timeStart <- as.character(Sys.time())

# Loading and checking required packages
library(methods, quietly=TRUE, warn.conflicts=FALSE)

inputType <- as.character(argv[1])
if (inputType=="fastq") {
  fastqPath <- as.character(gsub("fastq::", "", argv[grepl("fastq::", argv)], 
                                 fixed=TRUE))
  argv <- argv[!grepl("fastq::", argv, fixed=TRUE)]
  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])

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") {
  pairData <- character()
  pairData[2] <- as.character(argv[15])
  pairData[1] <- as.character(argv[16])
} else if (workMode=="glm") {

  selectOpt <- as.character(argv[18])
  selectVals <- as.character(argv[19])
}

# Read in inputs

samples <- read.table(samplePath, header=TRUE, sep="\t")
anno <- read.table(annoPath, header=TRUE, sep="\t")
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)))) {
  stop("'group' column not specified in sample annotation file")

if (any(table(samples$ID)>1)){
  tab <- table(samples$ID)
  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 (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))))) {

  } # Check that a group has be specifed for each sample
  
  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 (roastOpt == "yes") {
    if (is.na(match("Gene", colnames(anno)))) {
      tempStr <- paste("Gene-wise tests selected but'Gene' column not",
      "specified in hairpin annotation file")
      stop(tempStr)
    }
  }
}

################################################################################

# Process arguments
if (workMode=="glm") {
  if (roastOpt=="yes") {

if (exists("pairData")) {
  pairData <- make.names(pairData)
}

# Generate output folder and paths
dir.create(folderPath, showWarnings=FALSE)

# Generate links for outputs
imgOut("barHairpin")
imgOut("barIndex")
imgOut("mds")

    roastOut[i] <- makeOut(paste0("roast(", contrastData[i], ").tsv"))
    barcodePng[i] <- makeOut(paste0("barcode(", contrastData[i], ").png"))
    barcodePdf[i] <- makeOut(paste0("barcode(", contrastData[i], ").pdf"))
  }
}
countsOut <- makeOut("counts.tsv")

# Initialise data for html links and images, table with the link label and
# link address
linkData <- data.frame(Label=character(), Link=character(),
                       stringsAsFactors=FALSE)
imageData <- data.frame(Label=character(), Link=character(),

    selectVals <- unlist(strsplit(selectVals, " "))
  }                                                           
}
                                                  
if (inputType=="fastq") {
  # Use EdgeR hairpin process and capture outputs
  hpReadout <- capture.output(
  data <- processHairpinReads(fastqPath, samplePath, annoPath,
                              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") {
  # Process counts information, set ID column to be row names
  rownames(counts) <- counts$ID
  counts <- counts[ , !(colnames(counts)=="ID")]
  countsRows <- nrow(counts)

  anno <- anno[sel, ]
  
  # Create DGEList
  data <- DGEList(counts=counts, lib.size=colSums(counts), 
                  norm.factors=rep(1,ncol(counts)), genes=anno, group=factors)
  
  # Make the names of groups syntactically valid (replace spaces with periods)
  data$samples$group <- make.names(data$samples$group)
}

# Filter hairpins with low counts
preFilterCount <- nrow(data)
sel <- rowSums(cpm(data$counts) > cpmReq) >= sampleReq
data <- data[sel, ]
postFilterCount <- nrow(data)
filteredCount <- preFilterCount-postFilterCount

# Estimate dispersions
data <- estimateDisp(data)
commonBCV <- sqrt(data$common.dispersion)

      dev.off()
    }
  }
}

ID <- rownames(data$counts)
outputCounts <- cbind(ID, data$counts)
write.table(outputCounts, file=countsOut, row.names=FALSE, sep="\t")
linkName <- "Counts table (.tsv)"
linkAddr <- "counts.tsv"
linkData <- rbind(linkData, c(linkName, linkAddr))

# Record ending time and calculate total run time
timeEnd <- as.character(Sys.time())
timeTaken <- capture.output(round(difftime(timeEnd,timeStart), digits=3))
timeTaken <- gsub("Time difference of ", "", timeTaken, fixed=TRUE)
################################################################################
### HTML Generation
################################################################################
# Clear file
cat("", file=htmlPath)

if (inputType=="fastq") {
  cata("<ul>\n")
  ListItem(hpReadout[1])
  ListItem(hpReadout[2])
  cata("</ul>\n")
  cata(hpReadout[3], "<br />\n")
  cata("<ul>\n")
  ListItem(hpReadout[4])
  ListItem(hpReadout[7])
  cata("</ul>\n")
  cata(hpReadout[8:11], sep="<br />\n")
  cata("<br />\n")
  cata("<b>Please check that read percentages are consistent with ")
  cata("expectations.</b><br >\n")
} else if (inputType=="counts") {
  cata("<ul>\n")

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

cata("<h4>Output:</h4>\n")
cata("All images displayed have PDF copy at the bottom of the page, these can ")
cata("exported in a pdf viewer to high resolution image format. <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])
  }
}
cata("<br />\n")

cata("<h4>Plots:</h4>\n")
for (i in 1:nrow(linkData)) {
  if (!grepl(".tsv", linkData$Link[i])) {
    HtmlLink(linkData$Link[i], linkData$Label[i])

cata("<p>alt-click any of the links to download the file, or click the name ")
cata("of this task in the galaxy history panel and click on the floppy ")
cata("disk icon to download all files in a zip archive.</p>\n")
cata("<p>.tsv files are tab seperated files that can be viewed using Excel ")
cata("or other spreadsheet programs</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 that do not have more than", cpmReq,
                   "CPM in at least", sampleReq, "samples are considered",
                   "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.")
  ListItem(tempStr)
}

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



cit <- character()
link <-character()
link[1] <- paste0("<a href=\"",
                  "http://www.bioconductor.org/packages/release/bioc/",
                  "vignettes/limma/inst/doc/usersguide.pdf",
                  "\">", "limma User's Guide", "</a>.")
link[2] <- paste0("<a href=\"",
                  "http://www.bioconductor.org/packages/release/bioc/",
                  "vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf",
                  "\">", "edgeR User's Guide", "</a>")
                  
cit[1] <- paste("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")
cit[2] <- paste("Robinson MD and Smyth GK (2007). Moderated statistical tests",
                "for assessing differences in tag abundance. Bioinformatics",
                "23, 2881-2887")
cit[3] <- paste("Robinson MD and Smyth GK (2008). Small-sample estimation of",
                "negative binomial dispersion, with applications to SAGE data.",
                "Biostatistics, 9, 321-332")

cit[4] <- paste("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")

cata("<h4>Citations</h4>")
cata("<ol>\n")
ListItem(cit[1])
ListItem(cit[2])
ListItem(cit[3])
ListItem(cit[4])
cata("</ol>\n")

cata("<table border=\"0\">\n")
cata("<tr>\n")
TableItem("Task started at:"); TableItem(timeStart)
cata("</tr>\n")
cata("<tr>\n")
TableItem("Task ended at:"); TableItem(timeEnd)
cata("</tr>\n")
cata("<tr>\n")
TableItem("Task run time:"); TableItem(timeTaken)
cata("<tr>\n")
cata("</table>\n")

cata("</body>\n")
cata("</html>")