comparison deseq2.R @ 12:bd06df00180a draft

planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/deseq2 commit 62e9101c1e7b8467e395f31ccbd9321de01a6418

11:25204a860b79

  "quiet", "q", 0, "logical",
  "help", "h", 0, "logical",
  "outfile", "o", 1, "character",
  "countsfile", "n", 1, "character",
  "factors", "f", 1, "character",
  "plots" , "p", 1, "character",
  "sample_table", "s", 1, "character",
  "tximport", "i", 0, "logical",
  "txtype", "y", 1, "character",
  "tx2gene", "x", 1, "character", # a space-sep tx-to-gene map or GTF file (auto detect .gtf/.GTF)

if (!is.null(opt$factors)) {
  library("rjson")
  parser <- newJSONParser()
  parser$addData(opt$factors)
  factorList <- parser$getObject()
  factors <- sapply(factorList, function(x) x[[1]])
  primaryFactor <- factors[1]
  filenamesIn <- unname(unlist(factorList[[1]][[2]]))
  sampleTable <- data.frame(sample=basename(filenamesIn),
                            filename=filenamesIn,
                            row.names=filenamesIn,
                            stringsAsFactors=FALSE)
  for (factor in factorList) {

      files <- factor[[2]][[i]][[1]]
      sampleTable[files,factorName] <- trim(lvls[i])
    }
    sampleTable[[factorName]] <- factor(sampleTable[[factorName]], levels=lvls)
  }
  rownames(sampleTable) <- sampleTable$sample
} else {
  # read the sample_table argument
  # this table is described in ?DESeqDataSet
  # one column for the sample name, one for the filename, and
  # the remaining columns for factors in the analysis

  cat("\n\n")
}

# these are plots which are made once for each analysis
generateGenericPlots <- function(dds, factors) {
  rld <- rlog(dds)
  d=plotPCA(rld, intgroup=rev(factors), returnData=TRUE)
  labs <- paste0(seq_len(ncol(dds)), ": ", do.call(paste, as.list(colData(dds)[factors])))
  library(ggplot2)
  print(ggplot(d, aes(x=PC1,y=PC2, col=group,label=factor(labs)), environment = environment()) + geom_point() + geom_text(size=3))  
  dat <- assay(rld)
  colnames(dat) <- labs
  distsRL <- dist(t(dat))
  mat <- as.matrix(distsRL)
  hc <- hclust(distsRL)
  hmcol <- colorRampPalette(brewer.pal(9, "GnBu"))(100)
  heatmap.2(mat, Rowv=as.dendrogram(hc), symm=TRUE, trace="none", col = rev(hmcol),
            main="Sample-to-sample distances", margin=c(13,13))
  plotDispEsts(dds, main="Dispersion estimates")
}

# these are plots which can be made for each comparison, e.g.
# once for C vs A and once for B vs A

if (!useTXI) {
  # construct the object from HTSeq files
  dds <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable,
                                    directory = dir,
                                    design =  designFormula)
} else {
    # construct the object using tximport
    # first need to make the tx2gene table
    # this takes ~2-3 minutes using Bioconductor functions
    if (!is.null(gtfFile)) {

      df <- select(txdb, keys = k, keytype = "GENEID", columns = "TXNAME")
      tx2gene <- df[, 2:1]  # tx ID, then gene ID
    }
    library("tximport")
    txiFiles <- as.character(sampleTable[,2])
    names(txiFiles) <- as.character(sampleTable[,1])
    txi <- tximport(txiFiles, type=opt$txtype, tx2gene=tx2gene)
    dds <- DESeqDataSetFromTximport(txi,
                                    sampleTable[,3:ncol(sampleTable),drop=FALSE],
                                    designFormula)
}

n <- nlevels(colData(dds)[[primaryFactor]])
allLevels <- levels(colData(dds)[[primaryFactor]])

if (!is.null(opt$countsfile)) {
    labs <- paste0(seq_len(ncol(dds)), ": ", do.call(paste, as.list(colData(dds)[factors])))
    normalizedCounts<-counts(dds,normalized=TRUE)
    colnames(normalizedCounts)<-labs
    write.table(normalizedCounts, file=opt$countsfile, sep="\t", col.names=NA, quote=FALSE)
}

if (is.null(opt$many_contrasts)) {
  # only contrast the first and second level of the primary factor

12:bd06df00180a

  "quiet", "q", 0, "logical",
  "help", "h", 0, "logical",
  "outfile", "o", 1, "character",
  "countsfile", "n", 1, "character",
  "factors", "f", 1, "character",
  "files_to_labels", "l", 1, "character",
  "plots" , "p", 1, "character",
  "sample_table", "s", 1, "character",
  "tximport", "i", 0, "logical",
  "txtype", "y", 1, "character",
  "tx2gene", "x", 1, "character", # a space-sep tx-to-gene map or GTF file (auto detect .gtf/.GTF)

if (!is.null(opt$factors)) {
  library("rjson")
  parser <- newJSONParser()
  parser$addData(opt$factors)
  factorList <- parser$getObject()
  filenames_to_labels <- fromJSON(opt$files_to_labels)
  factors <- sapply(factorList, function(x) x[[1]])
  primaryFactor <- factors[1]
  filenamesIn <- unname(unlist(factorList[[1]][[2]]))
  labs = unname(unlist(filenames_to_labels[basename(filenamesIn)]))
  sampleTable <- data.frame(sample=basename(filenamesIn),
                            filename=filenamesIn,
                            row.names=filenamesIn,
                            stringsAsFactors=FALSE)
  for (factor in factorList) {

      files <- factor[[2]][[i]][[1]]
      sampleTable[files,factorName] <- trim(lvls[i])
    }
    sampleTable[[factorName]] <- factor(sampleTable[[factorName]], levels=lvls)
  }
  rownames(sampleTable) <- labs
} else {
  # read the sample_table argument
  # this table is described in ?DESeqDataSet
  # one column for the sample name, one for the filename, and
  # the remaining columns for factors in the analysis

  cat("\n\n")
}

# these are plots which are made once for each analysis
generateGenericPlots <- function(dds, factors) {
  library("ggplot2")
  library("ggrepel")
  library("pheatmap")

  rld <- rlog(dds)
  p <- plotPCA(rld, intgroup=rev(factors))
  print(p + geom_text_repel(aes_string(x = "PC1", y = "PC2", label = factor(colnames(dds))), size=3)  + geom_point())
  dat <- assay(rld)
  distsRL <- dist(t(dat))
  mat <- as.matrix(distsRL)
  colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
  pheatmap(mat,
           clustering_distance_rows=distsRL,
           clustering_distance_cols=distsRL,
           col=colors,
           main="Sample-to-sample distances")
  plotDispEsts(dds, main="Dispersion estimates")
}

# these are plots which can be made for each comparison, e.g.
# once for C vs A and once for B vs A

if (!useTXI) {
  # construct the object from HTSeq files
  dds <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable,
                                    directory = dir,
                                    design =  designFormula)
  labs <- unname(unlist(filenames_to_labels[colnames(dds)]))
  colnames(dds) <- labs
} else {
    # construct the object using tximport
    # first need to make the tx2gene table
    # this takes ~2-3 minutes using Bioconductor functions
    if (!is.null(gtfFile)) {

      df <- select(txdb, keys = k, keytype = "GENEID", columns = "TXNAME")
      tx2gene <- df[, 2:1]  # tx ID, then gene ID
    }
    library("tximport")
    txiFiles <- as.character(sampleTable[,2])
    labs <- unname(unlist(filenames_to_labels[sampleTable[,1]]))
    names(txiFiles) <- labs
    txi <- tximport(txiFiles, type=opt$txtype, tx2gene=tx2gene)
    dds <- DESeqDataSetFromTximport(txi,
                                    sampleTable[,3:ncol(sampleTable),drop=FALSE],
                                    designFormula)
}

n <- nlevels(colData(dds)[[primaryFactor]])
allLevels <- levels(colData(dds)[[primaryFactor]])

if (!is.null(opt$countsfile)) {
    normalizedCounts<-counts(dds,normalized=TRUE)
    write.table(normalizedCounts, file=opt$countsfile, sep="\t", col.names=NA, quote=FALSE)
}

if (is.null(opt$many_contrasts)) {
  # only contrast the first and second level of the primary factor