Mercurial > repos > shians > shrnaseq

comparison hairpinTool.R @ 8:548802b3492f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Version 1.0.9 - Added R session info to output - Added email for bug reports - Changed edgeR requirement to be stopping condition - Changed names of gene-wise comparison tables to "gene_level" from "roast" - Fixed selection of multiple genes for barcode plotting
author shian_su <registertonysu@gmail.com>
date Fri, 02 May 2014 17:22:24 +1000
parents 91e411fcdecc
children f1076bfb0ed1
comparison
equal deleted inserted replaced
7:91e411fcdecc 8:548802b3492f
56 library(splines, quietly=TRUE, warn.conflicts=FALSE) 56 library(splines, quietly=TRUE, warn.conflicts=FALSE)
57 library(edgeR, quietly=TRUE, warn.conflicts=FALSE) 57 library(edgeR, quietly=TRUE, warn.conflicts=FALSE)
58 library(limma, quietly=TRUE, warn.conflicts=FALSE) 58 library(limma, quietly=TRUE, warn.conflicts=FALSE)
59 59
60 if (packageVersion("edgeR") < "3.5.23") { 60 if (packageVersion("edgeR") < "3.5.23") {
61 message("Please update 'edgeR' to version >= 3.5.23 to run this script") 61 stop("Please update 'edgeR' to version >= 3.5.23 to run this tool")
62 }
63
64 if (packageVersion("limma")<"3.19.19") {
65 message("Update 'limma' to version >= 3.19.19 to see updated barcode graphs")
62 } 66 }
63 67
64 ################################################################################ 68 ################################################################################
65 ### Function declarations 69 ### Function declarations
66 ################################################################################ 70 ################################################################################
303 barcodePdf <- character() 307 barcodePdf <- character()
304 for (i in 1:length(contrastData)) { 308 for (i in 1:length(contrastData)) {
305 smearPng[i] <- makeOut(paste0("smear(", contrastData[i], ").png")) 309 smearPng[i] <- makeOut(paste0("smear(", contrastData[i], ").png"))
306 smearPdf[i] <- makeOut(paste0("smear(", contrastData[i], ").pdf")) 310 smearPdf[i] <- makeOut(paste0("smear(", contrastData[i], ").pdf"))
307 topOut[i] <- makeOut(paste0("toptag(", contrastData[i], ").tsv")) 311 topOut[i] <- makeOut(paste0("toptag(", contrastData[i], ").tsv"))
308 roastOut[i] <- makeOut(paste0("roast(", contrastData[i], ").tsv")) 312 roastOut[i] <- makeOut(paste0("gene_level(", contrastData[i], ").tsv"))
309 barcodePng[i] <- makeOut(paste0("barcode(", contrastData[i], ").png")) 313 barcodePng[i] <- makeOut(paste0("barcode(", contrastData[i], ").png"))
310 barcodePdf[i] <- makeOut(paste0("barcode(", contrastData[i], ").pdf")) 314 barcodePdf[i] <- makeOut(paste0("barcode(", contrastData[i], ").pdf"))
311 } 315 }
312 } 316 }
313 countsOut <- makeOut("counts.tsv") 317 countsOut <- makeOut("counts.tsv")
318 sessionOut <- makeOut("session_info.txt")
314 319
315 # Initialise data for html links and images, table with the link label and 320 # Initialise data for html links and images, table with the link label and
316 # link address 321 # link address
317 linkData <- data.frame(Label=character(), Link=character(), 322 linkData <- data.frame(Label=character(), Link=character(),
318 stringsAsFactors=FALSE) 323 stringsAsFactors=FALSE)
344 } 349 }
345 } 350 }
346 selectVals <- as.numeric(unique(selectVals)) 351 selectVals <- as.numeric(unique(selectVals))
347 } else { 352 } else {
348 selectVals <- gsub(" ", "", selectVals, fixed=TRUE) 353 selectVals <- gsub(" ", "", selectVals, fixed=TRUE)
349 selectVals <- unlist(strsplit(selectVals, " ")) 354 selectVals <- unlist(strsplit(selectVals, ","))
350 } 355 }
351 } 356 }
352 357
353 if (inputType=="fastq") { 358 if (inputType=="fastq") {
354 # Use EdgeR hairpin process and capture outputs 359 # Use EdgeR hairpin process and capture outputs
591 contrast=contrasts, nrot=nrot) 596 contrast=contrasts, nrot=nrot)
592 roastData <- cbind(GeneID=rownames(roastData), roastData) 597 roastData <- cbind(GeneID=rownames(roastData), roastData)
593 write.table(roastData, file=roastOut[i], row.names=FALSE, sep="\t") 598 write.table(roastData, file=roastOut[i], row.names=FALSE, sep="\t")
594 linkName <- paste0("Gene Level Analysis Table(", contrastData[i], 599 linkName <- paste0("Gene Level Analysis Table(", contrastData[i],
595 ") (.tsv)") 600 ") (.tsv)")
596 linkAddr <- paste0("roast(", contrastData[i], ").tsv") 601 linkAddr <- paste0("gene_level(", contrastData[i], ").tsv")
597 linkData <- rbind(linkData, c(linkName, linkAddr)) 602 linkData <- rbind(linkData, c(linkName, linkAddr))
598 if (selectOpt=="rank") { 603 if (selectOpt=="rank") {
599 selectedGenes <- rownames(roastData)[selectVals] 604 selectedGenes <- rownames(roastData)[selectVals]
600 } else { 605 } else {
601 selectedGenes <- selectVals 606 selectedGenes <- selectVals
634 dev.off() 639 dev.off()
635 } 640 }
636 } 641 }
637 } 642 }
638 643
644 # Generate data frame of the significant differences
639 sigDiff <- data.frame(Up=upCount, Flat=flatCount, Down=downCount) 645 sigDiff <- data.frame(Up=upCount, Flat=flatCount, Down=downCount)
640 if (workMode == "glm") { 646 if (workMode == "glm") {
641 row.names(sigDiff) <- contrastData 647 row.names(sigDiff) <- contrastData
642 } else if (workMode == "classic") { 648 } else if (workMode == "classic") {
643 row.names(sigDiff) <- paste0(pairData[2], "-", pairData[1]) 649 row.names(sigDiff) <- paste0(pairData[2], "-", pairData[1])
644 } 650 }
645 651
652 # Output table of summarised counts
646 ID <- rownames(data$counts) 653 ID <- rownames(data$counts)
647 outputCounts <- cbind(ID, data$counts) 654 outputCounts <- cbind(ID, data$counts)
648 write.table(outputCounts, file=countsOut, row.names=FALSE, sep="\t") 655 write.table(outputCounts, file=countsOut, row.names=FALSE, sep="\t",
656 quote=FALSE)
649 linkName <- "Counts table (.tsv)" 657 linkName <- "Counts table (.tsv)"
650 linkAddr <- "counts.tsv" 658 linkAddr <- "counts.tsv"
651 linkData <- rbind(linkData, c(linkName, linkAddr)) 659 linkData <- rbind(linkData, c(linkName, linkAddr))
660
661 # Record session info
662 writeLines(capture.output(sessionInfo()), sessionOut)
663 linkData <- rbind(linkData, c("Session Info", "session_info.txt"))
652 664
653 # Record ending time and calculate total run time 665 # Record ending time and calculate total run time
654 timeEnd <- as.character(Sys.time()) 666 timeEnd <- as.character(Sys.time())
655 timeTaken <- capture.output(round(difftime(timeEnd,timeStart), digits=3)) 667 timeTaken <- capture.output(round(difftime(timeEnd,timeStart), digits=3))
656 timeTaken <- gsub("Time difference of ", "", timeTaken, fixed=TRUE) 668 timeTaken <- gsub("Time difference of ", "", timeTaken, fixed=TRUE)
728 } 740 }
729 cata("</table>") 741 cata("</table>")
730 742
731 cata("<h4>Plots:</h4>\n") 743 cata("<h4>Plots:</h4>\n")
732 for (i in 1:nrow(linkData)) { 744 for (i in 1:nrow(linkData)) {
733 if (!grepl(".tsv", linkData$Link[i])) { 745 if (grepl(".pdf", linkData$Link[i])) {
734 HtmlLink(linkData$Link[i], linkData$Label[i]) 746 HtmlLink(linkData$Link[i], linkData$Label[i])
735 } 747 }
736 } 748 }
737 749
738 cata("<h4>Tables:</h4>\n") 750 cata("<h4>Tables:</h4>\n")
754 } else if (inputType == "counts") { 766 } else if (inputType == "counts") {
755 ListItem("Data was gathered from a table of counts.") 767 ListItem("Data was gathered from a table of counts.")
756 } 768 }
757 769
758 if (cpmReq!=0 && sampleReq!=0) { 770 if (cpmReq!=0 && sampleReq!=0) {
759 tempStr <- paste("Hairpins with less than", cpmReq, 771 tempStr <- paste("Hairpins without more than", cpmReq,
760 "CPM in at least", sampleReq, "samples are insignificant", 772 "CPM in at least", sampleReq, "samples are insignificant",
761 "and filtered out.") 773 "and filtered out.")
762 ListItem(tempStr) 774 ListItem(tempStr)
763 filterProp <- round(filteredCount/preFilterCount*100, digits=2) 775 filterProp <- round(filteredCount/preFilterCount*100, digits=2)
764 tempStr <- paste0(filteredCount, " of ", preFilterCount," (", filterProp, 776 tempStr <- paste0(filteredCount, " of ", preFilterCount," (", filterProp,
805 ListItem(cit[2]) 817 ListItem(cit[2])
806 ListItem(cit[3]) 818 ListItem(cit[3])
807 ListItem(cit[4]) 819 ListItem(cit[4])
808 cata("</ol>\n") 820 cata("</ol>\n")
809 821
822 cata("<p>Report problems to: su.s@wehi.edu.au</p>\n")
823
824 for (i in 1:nrow(linkData)) {
825 if (grepl("session_info", linkData$Link[i])) {
826 HtmlLink(linkData$Link[i], linkData$Label[i])
827 }
828 }
829
810 cata("<table border=\"0\">\n") 830 cata("<table border=\"0\">\n")
811 cata("<tr>\n") 831 cata("<tr>\n")
812 TableItem("Task started at:"); TableItem(timeStart) 832 TableItem("Task started at:"); TableItem(timeStart)
813 cata("</tr>\n") 833 cata("</tr>\n")
814 cata("<tr>\n") 834 cata("<tr>\n")