Remove unused file

examples/edger-pipeline.R

@@ -1,247 +0,0 @@
-#!/usr/bin/Rscript
-
-source("common.R")
-
-library(rtracklayer)
-library(DiffBind)
-library(plyr)
-library(doMC)
-registerDoMC()
-options(cores=multicore:::detectCores())
-library(snow)
-library(stringr)
-library(RColorBrewer)
-library(xlsx)
-library(edgeR)
-
-tsmsg <- function(...) {
-  message(date(), ": ", ...)
-}
-
-tsmsgf <- function(...) {
-  tsmsg(sprintf(...))
-}
-
-## Additional args are passed to every call of addDataFrame
-write.xlsx.multisheet <- function(data.frames, file, sheetNames=names(data.frames), ...) {
-  if (is.null(sheetNames)) {
-    sheetNames <- str_c("Sheet", seq_along(data.frames))
-  }
-  ## Ensure correct number of sheetNames
-  stopifnot(length(sheetNames) == length(data.frames))
-  ## Fill in missing names if needed
-  sheetNames[is.na(sheetNames)] <- str_c("Sheet", seq_along(data.frames))[is.na(sheetNames)]
-  wb <- createWorkbook()
-  sheets <- llply(sheetNames, function(x) createSheet(wb, sheetName=x))
-  mlply(cbind(sheet=sheets, x=data.frames), .fun=addDataFrame, .parallel=FALSE, ...)
-  saveWorkbook(wb, file)
-}
-
-renice.self <- function(niceness=19)  {
-  system2("renice", args=c("-n", as.character(niceness), Sys.getpid()))
-}
-
-makeNiceCluster <- function(...) {
-  cl <- makeCluster(...)
-  clusterCall(cl, renice.self)
-  cl
-}
-
-select.nearest <- function(x, y) {
-  y[nearest(x,y)]
-}
-
-kgxref <- get.ucsc.table("knownGene","kgXref")
-get.kgxref <- function(kgID) {
-  x <- merge(x=DataFrame(kgID=kgID), y=kgxref,
-             all.x=TRUE, all.y=FALSE)
-  x[names(x) != "kgID"]
-}
-
-read.htseq.counts <- function(f) {
-  x <- read.table(f, header=FALSE, sep="\t")
-  names(x) <- c("ID", "count")
-  ## Discard the last 5 lines
-  exception.rows <- (nrow(x)-4):nrow(x)
-  attr(x, "exception_counts") <- x[exception.rows,]
-  x <- x[-exception.rows,]
-  x
-}
-
-get.transcript.lengths <- function(transcript.ids, exon.lengths) {
-  exons.by.transcript <- split(exon.lengths, transcript.ids)
-  sapply(exons.by.transcript, sum)
-}
-
-get.gene.lengths <- function(gene.ids, transcript.lengths, method="max") {
-  if (is.character(method)) {
-    method <- get(method)
-  }
-  transcripts.by.genes <- split(transcript.lengths, gene.ids)
-  sapply(transcripts.by.genes, method)
-}
-
-str_matches <- function(string, pattern) {
-  !is.na(str_match(string, pattern)[,1])
-}
-
-annotate.ensp.in.hsap.or.mmul <- function(ensp.ids) {
-  ensembl=useMart("ensembl")
-  datasets <- c("mmulatta_gene_ensembl", "hsapiens_gene_ensembl")
-  x <- Reduce(rbind, llply(datasets, function(x) {
-    ensembl <- useDataset(x, mart=ensembl)
-    getBM(filters="ensembl_peptide_id",
-          values=unique(ensp.ids),
-          attributes=c("hgnc_symbol", "wikigene_name", "ensembl_gene_id", "ensembl_transcript_id", "ensembl_peptide_id", "description", "wikigene_description"),
-          mart=ensembl,
-          uniqueRows=TRUE)
-  }))
-  ## Convert empty string to NA in all columns
-  x <- data.frame(llply(x, function(column) ifelse(column == "", NA, column)))
-
-  ## Unify symbols & descriptions
-  unified.symbol <- ifelse(is.na(x$hgnc_symbol), as.character(x$wikigene_name), as.character(x$hgnc_symbol))
-  unified.desc <- ifelse(is.na(x$description), as.character(x$wikigene_description), as.character(x$description))
-  x <- data.frame(symbol=unified.symbol, x[! names(x) %in% c("hgnc_symbol", "wikigene_name", "description", "wikigene_description")], description=unified.desc)
-
-  ## Reorder rows so that more-preferred rows for the same input are on top
-  x <- x[order(is.na(x$symbol), str_matches(x$symbol, "^LOC\\d+$"), is.na(x$description)),]
-
-  ## Keep only the first row for each input
-  x <- x[!duplicated(x$ensembl_peptide_id),]
-
-  data.frame(x, row.names=x$ensembl_peptide_id)
-}
-
-## TODO: Replace by biomart
-
-CE.name.annot <- {
-  x <- setNames(nm=c("CE_ensembl_peptide_id", "symbol", "description"),
-                read.table("annotation/CE.name", sep="\t", quote=""))
-  x$ensembl_peptide_id <- str_replace(x$CE_ensembl_peptide_id, "^CE_", "")
-  row.names(x) <- x$ensembl_peptide_id
-  x
-}
-
-annotation.gff <- import("cuffmerge_results/merged.gff")
-
-transcripts <- annotation.gff[annotation.gff$type == "transcript", ]
-exons <- annotation.gff[annotation.gff$type == "exon", ]
-transcript.lengths <- get.transcript.lengths(unlist(exons$Parent), width(exons))
-gene.lengths <- get.gene.lengths(transcripts$geneID, transcript.lengths[transcripts$ID])
-
-gene.ref.lists <- sapply(split(transcripts$nearest_ref, transcripts$geneID),
-                         function(x) unique(x[!is.na(x)]))
-
-## For genes with multiple ref IDs, join them with commas
-gene.refs <- sapply(gene.ref.lists, function (x) {
-  if (length(x) == 0) {
-    NA
-  } else {
-    str_c(x, collapse=",")
-  }
-})
-
-## For genes with multiple ref IDs, just pick the first one
-gene.first.refs <- sapply(gene.ref.lists, function(x) {
-  if (length(x) == 0) {
-    NA
-  } else {
-    x[[1]]
-  }
-})
-
-## Use delayed assignment so that the cluster won't be created until it is actually needed
-delayedAssign("cl", makeNiceCluster(rep(c("salomon14", "salomon18", "salomon19"), 8)))
-
-mart.raw.annot <- annotate.ensp.in.hsap.or.mmul(na.omit(gene.first.refs))
-gene.annot <- data.frame(mart.raw.annot[gene.first.refs,], row.names=names(gene.first.refs))
-
-## Fill in missing values that are available from the CE.name file
-## provided by the cyno genome paper authors
-suppl.gene.annot <- data.frame(CE.name.annot[gene.first.refs,], row.names=names(gene.first.refs))
-
-suppl.ensp <- is.na(gene.annot$ensembl_peptide_id) & !is.na(suppl.gene.annot$ensembl_peptide_id)
-suppl.symbols <- is.na(gene.annot$symbol) & !is.na(suppl.gene.annot$symbol)
-gene.annot$ensembl_peptide_id <- ifelse(suppl.ensp,
-                                        as.character(suppl.gene.annot$ensembl_peptide_id),
-                                        as.character(gene.annot$ensembl_peptide_id))
-gene.annot$symbol <- ifelse(suppl.symbols,
-                            as.character(suppl.gene.annot$symbol),
-                            as.character(gene.annot$symbol))
-## Replace description whenever we replace the symbol to keep things consistent
-gene.annot$description <- ifelse(suppl.symbols,
-                            as.character(suppl.gene.annot$description),
-                            as.character(gene.annot$description))
-gene.annot$symbol[is.na(gene.annot$ensembl_peptide_id)] <- "[No annotation]"
-gene.annot$symbol[is.na(gene.annot$symbol)] <- "[No symbol for ENSP]"
-
-expdata <- {
-  x <- read.xlsx("kenyon-cyno-expdata.xls", 1)
-  x$Animal.ID <- str_trim(x$Animal.ID)
-  x$Condition <- str_replace_all(x$Condition, "γ", "g")
-  x$Sample.ID <-
-    sprintf("%s_%s_%s",
-            str_replace_all(x$Animal.ID, "[ -]", "_"),
-            ifelse(x$Condition == "Control", "CTRL", "IFNg"),
-            x$Passage)
-  x$Sample.Name <-
-    sprintf("%s_%s_%s_L%03i",
-            str_replace_all(x$Animal.ID, "[ -]", "_"),
-            ifelse(x$Condition == "Control", 1, 2),
-            x$Index.Seq,
-            x$Lane)
-  x$FASTQ.Read1.File <- sprintf("seqprep_results/%s/%s_R1_001.fastq", x$Sample.Name, x$Sample.Name)
-  x$FASTQ.Read2.File <- sprintf("seqprep_results/%s/%s_R2_001.fastq", x$Sample.Name, x$Sample.Name)
-  x$BAM.File <- sprintf("tophat_results/%s/accepted_hits.bam", x$Sample.Name)
-  x$GFF.File <- sprintf("cufflinks_quantification/%s/transcripts.gff", x$Sample.Name)
-  x$Counts.File <- sprintf("htseq_counts/%s/counts.txt", x$Sample.Name)
-  mcparallel(write.xlsx(x, "expdata.xlsx"))
-  rownames(x) <- x$Sample.ID
-  x
-}
-
-counts.vectors <- setNames(llply(expdata$Counts.File, .parallel=TRUE, function(f) {
-  x <- read.htseq.counts(f)
-  setNames(x$count, x$ID)
-}), expdata$Sample.ID)
-
-## Put the data into a matrix, making sure we account for the
-## possibility that not all genes are listed in all samples.
-all.geneIDs <- sort(unique(unlist(llply(counts.vectors, names))))
-
-counts <- {
-  x <- matrix(data=0, ncol=length(counts.vectors), nrow=length(all.geneIDs),
-              dimnames=list(`geneID`=all.geneIDs, `sample`=names(counts.vectors)))
-  for (i in expdata$Sample.ID) {
-    cvec <- counts.vectors[[i]]
-    x[names(cvec),i] <- cbind(cvec)
-  }
-  x
-}
-
-blocked.design <- model.matrix(~Animal.ID+Passage+Condition, data=expdata)
-unblocked.design <- model.matrix(~Condition, data=expdata)
-dge <- DGEList(counts=counts,
-               group=expdata$Condition,
-               genes=gene.annot[rownames(counts),])
-dge <- calcNormFactors(dge)
-
-blocked.dge <- estimateGLMCommonDisp(dge, blocked.design, verbose=TRUE)
-blocked.dge <- estimateGLMTrendedDisp(blocked.dge, blocked.design)
-blocked.dge <- estimateGLMTagwiseDisp(blocked.dge, blocked.design)
-blocked.fit <- glmFit(blocked.dge, blocked.design)
-blocked.lrt <- glmLRT(blocked.dge, blocked.fit, coef="ConditionIFNg activated")
-blocked.tt <- topTags(blocked.lrt, n=nrow(counts))
-
-unblocked.dge <- estimateGLMCommonDisp(dge, unblocked.design, verbose=TRUE)
-unblocked.dge <- estimateGLMTrendedDisp(unblocked.dge, unblocked.design)
-unblocked.dge <- estimateGLMTagwiseDisp(unblocked.dge, unblocked.design)
-unblocked.fit <- glmFit(unblocked.dge, unblocked.design)
-unblocked.lrt <- glmLRT(unblocked.dge, unblocked.fit, coef="ConditionIFNg activated")
-unblocked.tt <- topTags(unblocked.lrt, n=nrow(counts))
-
-write.csv(blocked.tt$table, "edgeR-genes-prelim")
-write.xlsx.multisheet(list(blocked=blocked.tt$table,
-                           unblocked=unblocked.tt$table),
-                      "edgeR-genes.xlsx")

examples/edger-pipeline.R.html

@@ -1,336 +0,0 @@
-<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01//EN"
-   "http://www.w3.org/TR/html4/strict.dtd">
-
-<html>
-<head>
-  <title></title>
-  <meta http-equiv="content-type" content="text/html; charset=utf-8">
-  <style type="text/css">
-td.linenos { background-color: #f0f0f0; padding-right: 10px; }
-span.lineno { background-color: #f0f0f0; padding: 0 5px 0 5px; }
-pre { line-height: 125%; }
-body .hll { background-color: #ffffcc }
-body  { background: #f8f8f8; }
-body .c { color: #408080; font-style: italic } /* Comment */
-body .err { border: 1px solid #FF0000 } /* Error */
-body .k { color: #008000; font-weight: bold } /* Keyword */
-body .o { color: #666666 } /* Operator */
-body .ch { color: #408080; font-style: italic } /* Comment.Hashbang */
-body .cm { color: #408080; font-style: italic } /* Comment.Multiline */
-body .cp { color: #BC7A00 } /* Comment.Preproc */
-body .cpf { color: #408080; font-style: italic } /* Comment.PreprocFile */
-body .c1 { color: #408080; font-style: italic } /* Comment.Single */
-body .cs { color: #408080; font-style: italic } /* Comment.Special */
-body .gd { color: #A00000 } /* Generic.Deleted */
-body .ge { font-style: italic } /* Generic.Emph */
-body .gr { color: #FF0000 } /* Generic.Error */
-body .gh { color: #000080; font-weight: bold } /* Generic.Heading */
-body .gi { color: #00A000 } /* Generic.Inserted */
-body .go { color: #888888 } /* Generic.Output */
-body .gp { color: #000080; font-weight: bold } /* Generic.Prompt */
-body .gs { font-weight: bold } /* Generic.Strong */
-body .gu { color: #800080; font-weight: bold } /* Generic.Subheading */
-body .gt { color: #0044DD } /* Generic.Traceback */
-body .kc { color: #008000; font-weight: bold } /* Keyword.Constant */
-body .kd { color: #008000; font-weight: bold } /* Keyword.Declaration */
-body .kn { color: #008000; font-weight: bold } /* Keyword.Namespace */
-body .kp { color: #008000 } /* Keyword.Pseudo */
-body .kr { color: #008000; font-weight: bold } /* Keyword.Reserved */
-body .kt { color: #B00040 } /* Keyword.Type */
-body .m { color: #666666 } /* Literal.Number */
-body .s { color: #BA2121 } /* Literal.String */
-body .na { color: #7D9029 } /* Name.Attribute */
-body .nb { color: #008000 } /* Name.Builtin */
-body .nc { color: #0000FF; font-weight: bold } /* Name.Class */
-body .no { color: #880000 } /* Name.Constant */
-body .nd { color: #AA22FF } /* Name.Decorator */
-body .ni { color: #999999; font-weight: bold } /* Name.Entity */
-body .ne { color: #D2413A; font-weight: bold } /* Name.Exception */
-body .nf { color: #0000FF } /* Name.Function */
-body .nl { color: #A0A000 } /* Name.Label */
-body .nn { color: #0000FF; font-weight: bold } /* Name.Namespace */
-body .nt { color: #008000; font-weight: bold } /* Name.Tag */
-body .nv { color: #19177C } /* Name.Variable */
-body .ow { color: #AA22FF; font-weight: bold } /* Operator.Word */
-body .w { color: #bbbbbb } /* Text.Whitespace */
-body .mb { color: #666666 } /* Literal.Number.Bin */
-body .mf { color: #666666 } /* Literal.Number.Float */
-body .mh { color: #666666 } /* Literal.Number.Hex */
-body .mi { color: #666666 } /* Literal.Number.Integer */
-body .mo { color: #666666 } /* Literal.Number.Oct */
-body .sa { color: #BA2121 } /* Literal.String.Affix */
-body .sb { color: #BA2121 } /* Literal.String.Backtick */
-body .sc { color: #BA2121 } /* Literal.String.Char */
-body .dl { color: #BA2121 } /* Literal.String.Delimiter */
-body .sd { color: #BA2121; font-style: italic } /* Literal.String.Doc */
-body .s2 { color: #BA2121 } /* Literal.String.Double */
-body .se { color: #BB6622; font-weight: bold } /* Literal.String.Escape */
-body .sh { color: #BA2121 } /* Literal.String.Heredoc */
-body .si { color: #BB6688; font-weight: bold } /* Literal.String.Interpol */
-body .sx { color: #008000 } /* Literal.String.Other */
-body .sr { color: #BB6688 } /* Literal.String.Regex */
-body .s1 { color: #BA2121 } /* Literal.String.Single */
-body .ss { color: #19177C } /* Literal.String.Symbol */
-body .bp { color: #008000 } /* Name.Builtin.Pseudo */
-body .fm { color: #0000FF } /* Name.Function.Magic */
-body .vc { color: #19177C } /* Name.Variable.Class */
-body .vg { color: #19177C } /* Name.Variable.Global */
-body .vi { color: #19177C } /* Name.Variable.Instance */
-body .vm { color: #19177C } /* Name.Variable.Magic */
-body .il { color: #666666 } /* Literal.Number.Integer.Long */
-
-  </style>
-</head>
-<body>
-<h2></h2>
-
-<div class="highlight"><pre><span></span><span class="c1">#!/usr/bin/Rscript</span>
-
-<span class="kn">source</span><span class="p">(</span><span class="s">&quot;common.R&quot;</span><span class="p">)</span>
-
-<span class="kn">library</span><span class="p">(</span>rtracklayer<span class="p">)</span>
-<span class="kn">library</span><span class="p">(</span>DiffBind<span class="p">)</span>
-<span class="kn">library</span><span class="p">(</span>plyr<span class="p">)</span>
-<span class="kn">library</span><span class="p">(</span>doMC<span class="p">)</span>
-registerDoMC<span class="p">()</span>
-<span class="kp">options</span><span class="p">(</span>cores<span class="o">=</span>multicore<span class="o">:::</span>detectCores<span class="p">())</span>
-<span class="kn">library</span><span class="p">(</span>snow<span class="p">)</span>
-<span class="kn">library</span><span class="p">(</span>stringr<span class="p">)</span>
-<span class="kn">library</span><span class="p">(</span>RColorBrewer<span class="p">)</span>
-<span class="kn">library</span><span class="p">(</span>xlsx<span class="p">)</span>
-<span class="kn">library</span><span class="p">(</span>edgeR<span class="p">)</span>
-
-tsmsg <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span><span class="kc">...</span><span class="p">)</span> <span class="p">{</span>
-  <span class="kp">message</span><span class="p">(</span><span class="kp">date</span><span class="p">(),</span> <span class="s">&quot;: &quot;</span><span class="p">,</span> <span class="kc">...</span><span class="p">)</span>
-<span class="p">}</span>
-
-tsmsgf <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span><span class="kc">...</span><span class="p">)</span> <span class="p">{</span>
-  tsmsg<span class="p">(</span><span class="kp">sprintf</span><span class="p">(</span><span class="kc">...</span><span class="p">))</span>
-<span class="p">}</span>
-
-<span class="c1">## Additional args are passed to every call of addDataFrame</span>
-write.xlsx.multisheet <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span>data.frames<span class="p">,</span> <span class="kp">file</span><span class="p">,</span> sheetNames<span class="o">=</span><span class="kp">names</span><span class="p">(</span>data.frames<span class="p">),</span> <span class="kc">...</span><span class="p">)</span> <span class="p">{</span>
-  <span class="kr">if</span> <span class="p">(</span><span class="kp">is.null</span><span class="p">(</span>sheetNames<span class="p">))</span> <span class="p">{</span>
-    sheetNames <span class="o">&lt;-</span> str_c<span class="p">(</span><span class="s">&quot;Sheet&quot;</span><span class="p">,</span> <span class="kp">seq_along</span><span class="p">(</span>data.frames<span class="p">))</span>
-  <span class="p">}</span>
-  <span class="c1">## Ensure correct number of sheetNames</span>
-  <span class="kp">stopifnot</span><span class="p">(</span><span class="kp">length</span><span class="p">(</span>sheetNames<span class="p">)</span> <span class="o">==</span> <span class="kp">length</span><span class="p">(</span>data.frames<span class="p">))</span>
-  <span class="c1">## Fill in missing names if needed</span>
-  sheetNames<span class="p">[</span><span class="kp">is.na</span><span class="p">(</span>sheetNames<span class="p">)]</span> <span class="o">&lt;-</span> str_c<span class="p">(</span><span class="s">&quot;Sheet&quot;</span><span class="p">,</span> <span class="kp">seq_along</span><span class="p">(</span>data.frames<span class="p">))[</span><span class="kp">is.na</span><span class="p">(</span>sheetNames<span class="p">)]</span>
-  wb <span class="o">&lt;-</span> createWorkbook<span class="p">()</span>
-  sheets <span class="o">&lt;-</span> llply<span class="p">(</span>sheetNames<span class="p">,</span> <span class="kr">function</span><span class="p">(</span>x<span class="p">)</span> createSheet<span class="p">(</span>wb<span class="p">,</span> sheetName<span class="o">=</span>x<span class="p">))</span>
-  mlply<span class="p">(</span><span class="kp">cbind</span><span class="p">(</span>sheet<span class="o">=</span>sheets<span class="p">,</span> x<span class="o">=</span>data.frames<span class="p">),</span> <span class="m">.</span>fun<span class="o">=</span>addDataFrame<span class="p">,</span> <span class="m">.</span>parallel<span class="o">=</span><span class="kc">FALSE</span><span class="p">,</span> <span class="kc">...</span><span class="p">)</span>
-  saveWorkbook<span class="p">(</span>wb<span class="p">,</span> <span class="kp">file</span><span class="p">)</span>
-<span class="p">}</span>
-
-renice.self <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span>niceness<span class="o">=</span><span class="m">19</span><span class="p">)</span>  <span class="p">{</span>
-  <span class="kp">system2</span><span class="p">(</span><span class="s">&quot;renice&quot;</span><span class="p">,</span> args<span class="o">=</span><span class="kt">c</span><span class="p">(</span><span class="s">&quot;-n&quot;</span><span class="p">,</span> <span class="kp">as.character</span><span class="p">(</span>niceness<span class="p">),</span> <span class="kp">Sys.getpid</span><span class="p">()))</span>
-<span class="p">}</span>
-
-makeNiceCluster <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span><span class="kc">...</span><span class="p">)</span> <span class="p">{</span>
-  cl <span class="o">&lt;-</span> makeCluster<span class="p">(</span><span class="kc">...</span><span class="p">)</span>
-  clusterCall<span class="p">(</span>cl<span class="p">,</span> renice.self<span class="p">)</span>
-  cl
-<span class="p">}</span>
-
-select.nearest <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span>x<span class="p">,</span> y<span class="p">)</span> <span class="p">{</span>
-  y<span class="p">[</span>nearest<span class="p">(</span>x<span class="p">,</span>y<span class="p">)]</span>
-<span class="p">}</span>
-
-kgxref <span class="o">&lt;-</span> get.ucsc.table<span class="p">(</span><span class="s">&quot;knownGene&quot;</span><span class="p">,</span><span class="s">&quot;kgXref&quot;</span><span class="p">)</span>
-get.kgxref <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span>kgID<span class="p">)</span> <span class="p">{</span>
-  x <span class="o">&lt;-</span> <span class="kp">merge</span><span class="p">(</span>x<span class="o">=</span>DataFrame<span class="p">(</span>kgID<span class="o">=</span>kgID<span class="p">),</span> y<span class="o">=</span>kgxref<span class="p">,</span>
-             all.x<span class="o">=</span><span class="kc">TRUE</span><span class="p">,</span> all.y<span class="o">=</span><span class="kc">FALSE</span><span class="p">)</span>
-  x<span class="p">[</span><span class="kp">names</span><span class="p">(</span>x<span class="p">)</span> <span class="o">!=</span> <span class="s">&quot;kgID&quot;</span><span class="p">]</span>
-<span class="p">}</span>
-
-read.htseq.counts <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span>f<span class="p">)</span> <span class="p">{</span>
-  x <span class="o">&lt;-</span> read.table<span class="p">(</span>f<span class="p">,</span> header<span class="o">=</span><span class="kc">FALSE</span><span class="p">,</span> sep<span class="o">=</span><span class="s">&quot;\t&quot;</span><span class="p">)</span>
-  <span class="kp">names</span><span class="p">(</span>x<span class="p">)</span> <span class="o">&lt;-</span> <span class="kt">c</span><span class="p">(</span><span class="s">&quot;ID&quot;</span><span class="p">,</span> <span class="s">&quot;count&quot;</span><span class="p">)</span>
-  <span class="c1">## Discard the last 5 lines</span>
-  exception.rows <span class="o">&lt;-</span> <span class="p">(</span><span class="kp">nrow</span><span class="p">(</span>x<span class="p">)</span><span class="m">-4</span><span class="p">)</span><span class="o">:</span><span class="kp">nrow</span><span class="p">(</span>x<span class="p">)</span>
-  <span class="kp">attr</span><span class="p">(</span>x<span class="p">,</span> <span class="s">&quot;exception_counts&quot;</span><span class="p">)</span> <span class="o">&lt;-</span> x<span class="p">[</span>exception.rows<span class="p">,]</span>
-  x <span class="o">&lt;-</span> x<span class="p">[</span><span class="o">-</span>exception.rows<span class="p">,]</span>
-  x
-<span class="p">}</span>
-
-get.transcript.lengths <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span>transcript.ids<span class="p">,</span> exon.lengths<span class="p">)</span> <span class="p">{</span>
-  exons.by.transcript <span class="o">&lt;-</span> <span class="kp">split</span><span class="p">(</span>exon.lengths<span class="p">,</span> transcript.ids<span class="p">)</span>
-  <span class="kp">sapply</span><span class="p">(</span>exons.by.transcript<span class="p">,</span> <span class="kp">sum</span><span class="p">)</span>
-<span class="p">}</span>
-
-get.gene.lengths <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span>gene.ids<span class="p">,</span> transcript.lengths<span class="p">,</span> method<span class="o">=</span><span class="s">&quot;max&quot;</span><span class="p">)</span> <span class="p">{</span>
-  <span class="kr">if</span> <span class="p">(</span><span class="kp">is.character</span><span class="p">(</span>method<span class="p">))</span> <span class="p">{</span>
-    method <span class="o">&lt;-</span> <span class="kp">get</span><span class="p">(</span>method<span class="p">)</span>
-  <span class="p">}</span>
-  transcripts.by.genes <span class="o">&lt;-</span> <span class="kp">split</span><span class="p">(</span>transcript.lengths<span class="p">,</span> gene.ids<span class="p">)</span>
-  <span class="kp">sapply</span><span class="p">(</span>transcripts.by.genes<span class="p">,</span> method<span class="p">)</span>
-<span class="p">}</span>
-
-str_matches <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span>string<span class="p">,</span> pattern<span class="p">)</span> <span class="p">{</span>
-  <span class="o">!</span><span class="kp">is.na</span><span class="p">(</span>str_match<span class="p">(</span>string<span class="p">,</span> pattern<span class="p">)[,</span><span class="m">1</span><span class="p">])</span>
-<span class="p">}</span>
-
-annotate.ensp.in.hsap.or.mmul <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span>ensp.ids<span class="p">)</span> <span class="p">{</span>
-  ensembl<span class="o">=</span>useMart<span class="p">(</span><span class="s">&quot;ensembl&quot;</span><span class="p">)</span>
-  datasets <span class="o">&lt;-</span> <span class="kt">c</span><span class="p">(</span><span class="s">&quot;mmulatta_gene_ensembl&quot;</span><span class="p">,</span> <span class="s">&quot;hsapiens_gene_ensembl&quot;</span><span class="p">)</span>
-  x <span class="o">&lt;-</span> <span class="kp">Reduce</span><span class="p">(</span><span class="kp">rbind</span><span class="p">,</span> llply<span class="p">(</span>datasets<span class="p">,</span> <span class="kr">function</span><span class="p">(</span>x<span class="p">)</span> <span class="p">{</span>
-    ensembl <span class="o">&lt;-</span> useDataset<span class="p">(</span>x<span class="p">,</span> mart<span class="o">=</span>ensembl<span class="p">)</span>
-    getBM<span class="p">(</span>filters<span class="o">=</span><span class="s">&quot;ensembl_peptide_id&quot;</span><span class="p">,</span>
-          values<span class="o">=</span><span class="kp">unique</span><span class="p">(</span>ensp.ids<span class="p">),</span>
-          attributes<span class="o">=</span><span class="kt">c</span><span class="p">(</span><span class="s">&quot;hgnc_symbol&quot;</span><span class="p">,</span> <span class="s">&quot;wikigene_name&quot;</span><span class="p">,</span> <span class="s">&quot;ensembl_gene_id&quot;</span><span class="p">,</span> <span class="s">&quot;ensembl_transcript_id&quot;</span><span class="p">,</span> <span class="s">&quot;ensembl_peptide_id&quot;</span><span class="p">,</span> <span class="s">&quot;description&quot;</span><span class="p">,</span> <span class="s">&quot;wikigene_description&quot;</span><span class="p">),</span>
-          mart<span class="o">=</span>ensembl<span class="p">,</span>
-          uniqueRows<span class="o">=</span><span class="kc">TRUE</span><span class="p">)</span>
-  <span class="p">}))</span>
-  <span class="c1">## Convert empty string to NA in all columns</span>
-  x <span class="o">&lt;-</span> <span class="kt">data.frame</span><span class="p">(</span>llply<span class="p">(</span>x<span class="p">,</span> <span class="kr">function</span><span class="p">(</span>column<span class="p">)</span> <span class="kp">ifelse</span><span class="p">(</span>column <span class="o">==</span> <span class="s">&quot;&quot;</span><span class="p">,</span> <span class="kc">NA</span><span class="p">,</span> column<span class="p">)))</span>
-
-  <span class="c1">## Unify symbols &amp; descriptions</span>
-  unified.symbol <span class="o">&lt;-</span> <span class="kp">ifelse</span><span class="p">(</span><span class="kp">is.na</span><span class="p">(</span>x<span class="o">$</span>hgnc_symbol<span class="p">),</span> <span class="kp">as.character</span><span class="p">(</span>x<span class="o">$</span>wikigene_name<span class="p">),</span> <span class="kp">as.character</span><span class="p">(</span>x<span class="o">$</span>hgnc_symbol<span class="p">))</span>
-  unified.desc <span class="o">&lt;-</span> <span class="kp">ifelse</span><span class="p">(</span><span class="kp">is.na</span><span class="p">(</span>x<span class="o">$</span>description<span class="p">),</span> <span class="kp">as.character</span><span class="p">(</span>x<span class="o">$</span>wikigene_description<span class="p">),</span> <span class="kp">as.character</span><span class="p">(</span>x<span class="o">$</span>description<span class="p">))</span>
-  x <span class="o">&lt;-</span> <span class="kt">data.frame</span><span class="p">(</span>symbol<span class="o">=</span>unified.symbol<span class="p">,</span> x<span class="p">[</span><span class="o">!</span> <span class="kp">names</span><span class="p">(</span>x<span class="p">)</span> <span class="o">%in%</span> <span class="kt">c</span><span class="p">(</span><span class="s">&quot;hgnc_symbol&quot;</span><span class="p">,</span> <span class="s">&quot;wikigene_name&quot;</span><span class="p">,</span> <span class="s">&quot;description&quot;</span><span class="p">,</span> <span class="s">&quot;wikigene_description&quot;</span><span class="p">)],</span> description<span class="o">=</span>unified.desc<span class="p">)</span>
-
-  <span class="c1">## Reorder rows so that more-preferred rows for the same input are on top</span>
-  x <span class="o">&lt;-</span> x<span class="p">[</span><span class="kp">order</span><span class="p">(</span><span class="kp">is.na</span><span class="p">(</span>x<span class="o">$</span>symbol<span class="p">),</span> str_matches<span class="p">(</span>x<span class="o">$</span>symbol<span class="p">,</span> <span class="s">&quot;^LOC\\d+$&quot;</span><span class="p">),</span> <span class="kp">is.na</span><span class="p">(</span>x<span class="o">$</span>description<span class="p">)),]</span>
-
-  <span class="c1">## Keep only the first row for each input</span>
-  x <span class="o">&lt;-</span> x<span class="p">[</span><span class="o">!</span><span class="kp">duplicated</span><span class="p">(</span>x<span class="o">$</span>ensembl_peptide_id<span class="p">),]</span>
-
-  <span class="kt">data.frame</span><span class="p">(</span>x<span class="p">,</span> row.names<span class="o">=</span>x<span class="o">$</span>ensembl_peptide_id<span class="p">)</span>
-<span class="p">}</span>
-
-<span class="c1">## TODO: Replace by biomart</span>
-
-CE.name.annot <span class="o">&lt;-</span> <span class="p">{</span>
-  x <span class="o">&lt;-</span> setNames<span class="p">(</span>nm<span class="o">=</span><span class="kt">c</span><span class="p">(</span><span class="s">&quot;CE_ensembl_peptide_id&quot;</span><span class="p">,</span> <span class="s">&quot;symbol&quot;</span><span class="p">,</span> <span class="s">&quot;description&quot;</span><span class="p">),</span>
-                read.table<span class="p">(</span><span class="s">&quot;annotation/CE.name&quot;</span><span class="p">,</span> sep<span class="o">=</span><span class="s">&quot;\t&quot;</span><span class="p">,</span> quote<span class="o">=</span><span class="s">&quot;&quot;</span><span class="p">))</span>
-  x<span class="o">$</span>ensembl_peptide_id <span class="o">&lt;-</span> str_replace<span class="p">(</span>x<span class="o">$</span>CE_ensembl_peptide_id<span class="p">,</span> <span class="s">&quot;^CE_&quot;</span><span class="p">,</span> <span class="s">&quot;&quot;</span><span class="p">)</span>
-  <span class="kp">row.names</span><span class="p">(</span>x<span class="p">)</span> <span class="o">&lt;-</span> x<span class="o">$</span>ensembl_peptide_id
-  x
-<span class="p">}</span>
-
-annotation.gff <span class="o">&lt;-</span> import<span class="p">(</span><span class="s">&quot;cuffmerge_results/merged.gff&quot;</span><span class="p">)</span>
-
-transcripts <span class="o">&lt;-</span> annotation.gff<span class="p">[</span>annotation.gff<span class="o">$</span>type <span class="o">==</span> <span class="s">&quot;transcript&quot;</span><span class="p">,</span> <span class="p">]</span>
-exons <span class="o">&lt;-</span> annotation.gff<span class="p">[</span>annotation.gff<span class="o">$</span>type <span class="o">==</span> <span class="s">&quot;exon&quot;</span><span class="p">,</span> <span class="p">]</span>
-transcript.lengths <span class="o">&lt;-</span> get.transcript.lengths<span class="p">(</span><span class="kp">unlist</span><span class="p">(</span>exons<span class="o">$</span>Parent<span class="p">),</span> width<span class="p">(</span>exons<span class="p">))</span>
-gene.lengths <span class="o">&lt;-</span> get.gene.lengths<span class="p">(</span>transcripts<span class="o">$</span>geneID<span class="p">,</span> transcript.lengths<span class="p">[</span>transcripts<span class="o">$</span>ID<span class="p">])</span>
-
-gene.ref.lists <span class="o">&lt;-</span> <span class="kp">sapply</span><span class="p">(</span><span class="kp">split</span><span class="p">(</span>transcripts<span class="o">$</span>nearest_ref<span class="p">,</span> transcripts<span class="o">$</span>geneID<span class="p">),</span>
-                         <span class="kr">function</span><span class="p">(</span>x<span class="p">)</span> <span class="kp">unique</span><span class="p">(</span>x<span class="p">[</span><span class="o">!</span><span class="kp">is.na</span><span class="p">(</span>x<span class="p">)]))</span>
-
-<span class="c1">## For genes with multiple ref IDs, join them with commas</span>
-gene.refs <span class="o">&lt;-</span> <span class="kp">sapply</span><span class="p">(</span>gene.ref.lists<span class="p">,</span> <span class="kr">function</span> <span class="p">(</span>x<span class="p">)</span> <span class="p">{</span>
-  <span class="kr">if</span> <span class="p">(</span><span class="kp">length</span><span class="p">(</span>x<span class="p">)</span> <span class="o">==</span> <span class="m">0</span><span class="p">)</span> <span class="p">{</span>
-    <span class="kc">NA</span>
-  <span class="p">}</span> <span class="kr">else</span> <span class="p">{</span>
-    str_c<span class="p">(</span>x<span class="p">,</span> collapse<span class="o">=</span><span class="s">&quot;,&quot;</span><span class="p">)</span>
-  <span class="p">}</span>
-<span class="p">})</span>
-
-<span class="c1">## For genes with multiple ref IDs, just pick the first one</span>
-gene.first.refs <span class="o">&lt;-</span> <span class="kp">sapply</span><span class="p">(</span>gene.ref.lists<span class="p">,</span> <span class="kr">function</span><span class="p">(</span>x<span class="p">)</span> <span class="p">{</span>
-  <span class="kr">if</span> <span class="p">(</span><span class="kp">length</span><span class="p">(</span>x<span class="p">)</span> <span class="o">==</span> <span class="m">0</span><span class="p">)</span> <span class="p">{</span>
-    <span class="kc">NA</span>
-  <span class="p">}</span> <span class="kr">else</span> <span class="p">{</span>
-    x<span class="p">[[</span><span class="m">1</span><span class="p">]]</span>
-  <span class="p">}</span>
-<span class="p">})</span>
-
-<span class="c1">## Use delayed assignment so that the cluster won&#39;t be created until it is actually needed</span>
-<span class="kp">delayedAssign</span><span class="p">(</span><span class="s">&quot;cl&quot;</span><span class="p">,</span> makeNiceCluster<span class="p">(</span><span class="kp">rep</span><span class="p">(</span><span class="kt">c</span><span class="p">(</span><span class="s">&quot;salomon14&quot;</span><span class="p">,</span> <span class="s">&quot;salomon18&quot;</span><span class="p">,</span> <span class="s">&quot;salomon19&quot;</span><span class="p">),</span> <span class="m">8</span><span class="p">)))</span>
-
-mart.raw.annot <span class="o">&lt;-</span> annotate.ensp.in.hsap.or.mmul<span class="p">(</span>na.omit<span class="p">(</span>gene.first.refs<span class="p">))</span>
-gene.annot <span class="o">&lt;-</span> <span class="kt">data.frame</span><span class="p">(</span>mart.raw.annot<span class="p">[</span>gene.first.refs<span class="p">,],</span> row.names<span class="o">=</span><span class="kp">names</span><span class="p">(</span>gene.first.refs<span class="p">))</span>
-
-<span class="c1">## Fill in missing values that are available from the CE.name file</span>
-<span class="c1">## provided by the cyno genome paper authors</span>
-suppl.gene.annot <span class="o">&lt;-</span> <span class="kt">data.frame</span><span class="p">(</span>CE.name.annot<span class="p">[</span>gene.first.refs<span class="p">,],</span> row.names<span class="o">=</span><span class="kp">names</span><span class="p">(</span>gene.first.refs<span class="p">))</span>
-
-suppl.ensp <span class="o">&lt;-</span> <span class="kp">is.na</span><span class="p">(</span>gene.annot<span class="o">$</span>ensembl_peptide_id<span class="p">)</span> <span class="o">&amp;</span> <span class="o">!</span><span class="kp">is.na</span><span class="p">(</span>suppl.gene.annot<span class="o">$</span>ensembl_peptide_id<span class="p">)</span>
-suppl.symbols <span class="o">&lt;-</span> <span class="kp">is.na</span><span class="p">(</span>gene.annot<span class="o">$</span>symbol<span class="p">)</span> <span class="o">&amp;</span> <span class="o">!</span><span class="kp">is.na</span><span class="p">(</span>suppl.gene.annot<span class="o">$</span>symbol<span class="p">)</span>
-gene.annot<span class="o">$</span>ensembl_peptide_id <span class="o">&lt;-</span> <span class="kp">ifelse</span><span class="p">(</span>suppl.ensp<span class="p">,</span>
-                                        <span class="kp">as.character</span><span class="p">(</span>suppl.gene.annot<span class="o">$</span>ensembl_peptide_id<span class="p">),</span>
-                                        <span class="kp">as.character</span><span class="p">(</span>gene.annot<span class="o">$</span>ensembl_peptide_id<span class="p">))</span>
-gene.annot<span class="o">$</span>symbol <span class="o">&lt;-</span> <span class="kp">ifelse</span><span class="p">(</span>suppl.symbols<span class="p">,</span>
-                            <span class="kp">as.character</span><span class="p">(</span>suppl.gene.annot<span class="o">$</span>symbol<span class="p">),</span>
-                            <span class="kp">as.character</span><span class="p">(</span>gene.annot<span class="o">$</span>symbol<span class="p">))</span>
-<span class="c1">## Replace description whenever we replace the symbol to keep things consistent</span>
-gene.annot<span class="o">$</span>description <span class="o">&lt;-</span> <span class="kp">ifelse</span><span class="p">(</span>suppl.symbols<span class="p">,</span>
-                            <span class="kp">as.character</span><span class="p">(</span>suppl.gene.annot<span class="o">$</span>description<span class="p">),</span>
-                            <span class="kp">as.character</span><span class="p">(</span>gene.annot<span class="o">$</span>description<span class="p">))</span>
-gene.annot<span class="o">$</span>symbol<span class="p">[</span><span class="kp">is.na</span><span class="p">(</span>gene.annot<span class="o">$</span>ensembl_peptide_id<span class="p">)]</span> <span class="o">&lt;-</span> <span class="s">&quot;[No annotation]&quot;</span>
-gene.annot<span class="o">$</span>symbol<span class="p">[</span><span class="kp">is.na</span><span class="p">(</span>gene.annot<span class="o">$</span>symbol<span class="p">)]</span> <span class="o">&lt;-</span> <span class="s">&quot;[No symbol for ENSP]&quot;</span>
-
-expdata <span class="o">&lt;-</span> <span class="p">{</span>
-  x <span class="o">&lt;-</span> read.xlsx<span class="p">(</span><span class="s">&quot;kenyon-cyno-expdata.xls&quot;</span><span class="p">,</span> <span class="m">1</span><span class="p">)</span>
-  x<span class="o">$</span>Animal.ID <span class="o">&lt;-</span> str_trim<span class="p">(</span>x<span class="o">$</span>Animal.ID<span class="p">)</span>
-  x<span class="o">$</span>Condition <span class="o">&lt;-</span> str_replace_all<span class="p">(</span>x<span class="o">$</span>Condition<span class="p">,</span> <span class="s">&quot;γ&quot;</span><span class="p">,</span> <span class="s">&quot;g&quot;</span><span class="p">)</span>
-  x<span class="o">$</span>Sample.ID <span class="o">&lt;-</span>
-    <span class="kp">sprintf</span><span class="p">(</span><span class="s">&quot;%s_%s_%s&quot;</span><span class="p">,</span>
-            str_replace_all<span class="p">(</span>x<span class="o">$</span>Animal.ID<span class="p">,</span> <span class="s">&quot;[ -]&quot;</span><span class="p">,</span> <span class="s">&quot;_&quot;</span><span class="p">),</span>
-            <span class="kp">ifelse</span><span class="p">(</span>x<span class="o">$</span>Condition <span class="o">==</span> <span class="s">&quot;Control&quot;</span><span class="p">,</span> <span class="s">&quot;CTRL&quot;</span><span class="p">,</span> <span class="s">&quot;IFNg&quot;</span><span class="p">),</span>
-            x<span class="o">$</span>Passage<span class="p">)</span>
-  x<span class="o">$</span>Sample.Name <span class="o">&lt;-</span>
-    <span class="kp">sprintf</span><span class="p">(</span><span class="s">&quot;%s_%s_%s_L%03i&quot;</span><span class="p">,</span>
-            str_replace_all<span class="p">(</span>x<span class="o">$</span>Animal.ID<span class="p">,</span> <span class="s">&quot;[ -]&quot;</span><span class="p">,</span> <span class="s">&quot;_&quot;</span><span class="p">),</span>
-            <span class="kp">ifelse</span><span class="p">(</span>x<span class="o">$</span>Condition <span class="o">==</span> <span class="s">&quot;Control&quot;</span><span class="p">,</span> <span class="m">1</span><span class="p">,</span> <span class="m">2</span><span class="p">),</span>
-            x<span class="o">$</span>Index.Seq<span class="p">,</span>
-            x<span class="o">$</span>Lane<span class="p">)</span>
-  x<span class="o">$</span>FASTQ.Read1.File <span class="o">&lt;-</span> <span class="kp">sprintf</span><span class="p">(</span><span class="s">&quot;seqprep_results/%s/%s_R1_001.fastq&quot;</span><span class="p">,</span> x<span class="o">$</span>Sample.Name<span class="p">,</span> x<span class="o">$</span>Sample.Name<span class="p">)</span>
-  x<span class="o">$</span>FASTQ.Read2.File <span class="o">&lt;-</span> <span class="kp">sprintf</span><span class="p">(</span><span class="s">&quot;seqprep_results/%s/%s_R2_001.fastq&quot;</span><span class="p">,</span> x<span class="o">$</span>Sample.Name<span class="p">,</span> x<span class="o">$</span>Sample.Name<span class="p">)</span>
-  x<span class="o">$</span>BAM.File <span class="o">&lt;-</span> <span class="kp">sprintf</span><span class="p">(</span><span class="s">&quot;tophat_results/%s/accepted_hits.bam&quot;</span><span class="p">,</span> x<span class="o">$</span>Sample.Name<span class="p">)</span>
-  x<span class="o">$</span>GFF.File <span class="o">&lt;-</span> <span class="kp">sprintf</span><span class="p">(</span><span class="s">&quot;cufflinks_quantification/%s/transcripts.gff&quot;</span><span class="p">,</span> x<span class="o">$</span>Sample.Name<span class="p">)</span>
-  x<span class="o">$</span>Counts.File <span class="o">&lt;-</span> <span class="kp">sprintf</span><span class="p">(</span><span class="s">&quot;htseq_counts/%s/counts.txt&quot;</span><span class="p">,</span> x<span class="o">$</span>Sample.Name<span class="p">)</span>
-  mcparallel<span class="p">(</span>write.xlsx<span class="p">(</span>x<span class="p">,</span> <span class="s">&quot;expdata.xlsx&quot;</span><span class="p">))</span>
-  <span class="kp">rownames</span><span class="p">(</span>x<span class="p">)</span> <span class="o">&lt;-</span> x<span class="o">$</span>Sample.ID
-  x
-<span class="p">}</span>
-
-counts.vectors <span class="o">&lt;-</span> setNames<span class="p">(</span>llply<span class="p">(</span>expdata<span class="o">$</span>Counts.File<span class="p">,</span> <span class="m">.</span>parallel<span class="o">=</span><span class="kc">TRUE</span><span class="p">,</span> <span class="kr">function</span><span class="p">(</span>f<span class="p">)</span> <span class="p">{</span>
-  x <span class="o">&lt;-</span> read.htseq.counts<span class="p">(</span>f<span class="p">)</span>
-  setNames<span class="p">(</span>x<span class="o">$</span>count<span class="p">,</span> x<span class="o">$</span>ID<span class="p">)</span>
-<span class="p">}),</span> expdata<span class="o">$</span>Sample.ID<span class="p">)</span>
-
-<span class="c1">## Put the data into a matrix, making sure we account for the</span>
-<span class="c1">## possibility that not all genes are listed in all samples.</span>
-all.geneIDs <span class="o">&lt;-</span> <span class="kp">sort</span><span class="p">(</span><span class="kp">unique</span><span class="p">(</span><span class="kp">unlist</span><span class="p">(</span>llply<span class="p">(</span>counts.vectors<span class="p">,</span> <span class="kp">names</span><span class="p">))))</span>
-
-counts <span class="o">&lt;-</span> <span class="p">{</span>
-  x <span class="o">&lt;-</span> <span class="kt">matrix</span><span class="p">(</span>data<span class="o">=</span><span class="m">0</span><span class="p">,</span> ncol<span class="o">=</span><span class="kp">length</span><span class="p">(</span>counts.vectors<span class="p">),</span> nrow<span class="o">=</span><span class="kp">length</span><span class="p">(</span>all.geneIDs<span class="p">),</span>
-              dimnames<span class="o">=</span><span class="kt">list</span><span class="p">(</span><span class="sb">`geneID`</span><span class="o">=</span>all.geneIDs<span class="p">,</span> <span class="sb">`sample`</span><span class="o">=</span><span class="kp">names</span><span class="p">(</span>counts.vectors<span class="p">)))</span>
-  <span class="kr">for</span> <span class="p">(</span>i <span class="kr">in</span> expdata<span class="o">$</span>Sample.ID<span class="p">)</span> <span class="p">{</span>
-    cvec <span class="o">&lt;-</span> counts.vectors<span class="p">[[</span>i<span class="p">]]</span>
-    x<span class="p">[</span><span class="kp">names</span><span class="p">(</span>cvec<span class="p">),</span>i<span class="p">]</span> <span class="o">&lt;-</span> <span class="kp">cbind</span><span class="p">(</span>cvec<span class="p">)</span>
-  <span class="p">}</span>
-  x
-<span class="p">}</span>
-
-blocked.design <span class="o">&lt;-</span> model.matrix<span class="p">(</span><span class="o">~</span>Animal.ID<span class="o">+</span>Passage<span class="o">+</span>Condition<span class="p">,</span> data<span class="o">=</span>expdata<span class="p">)</span>
-unblocked.design <span class="o">&lt;-</span> model.matrix<span class="p">(</span><span class="o">~</span>Condition<span class="p">,</span> data<span class="o">=</span>expdata<span class="p">)</span>
-dge <span class="o">&lt;-</span> DGEList<span class="p">(</span>counts<span class="o">=</span>counts<span class="p">,</span>
-               group<span class="o">=</span>expdata<span class="o">$</span>Condition<span class="p">,</span>
-               genes<span class="o">=</span>gene.annot<span class="p">[</span><span class="kp">rownames</span><span class="p">(</span>counts<span class="p">),])</span>
-dge <span class="o">&lt;-</span> calcNormFactors<span class="p">(</span>dge<span class="p">)</span>
-
-blocked.dge <span class="o">&lt;-</span> estimateGLMCommonDisp<span class="p">(</span>dge<span class="p">,</span> blocked.design<span class="p">,</span> verbose<span class="o">=</span><span class="kc">TRUE</span><span class="p">)</span>
-blocked.dge <span class="o">&lt;-</span> estimateGLMTrendedDisp<span class="p">(</span>blocked.dge<span class="p">,</span> blocked.design<span class="p">)</span>
-blocked.dge <span class="o">&lt;-</span> estimateGLMTagwiseDisp<span class="p">(</span>blocked.dge<span class="p">,</span> blocked.design<span class="p">)</span>
-blocked.fit <span class="o">&lt;-</span> glmFit<span class="p">(</span>blocked.dge<span class="p">,</span> blocked.design<span class="p">)</span>
-blocked.lrt <span class="o">&lt;-</span> glmLRT<span class="p">(</span>blocked.dge<span class="p">,</span> blocked.fit<span class="p">,</span> coef<span class="o">=</span><span class="s">&quot;ConditionIFNg activated&quot;</span><span class="p">)</span>
-blocked.tt <span class="o">&lt;-</span> topTags<span class="p">(</span>blocked.lrt<span class="p">,</span> n<span class="o">=</span><span class="kp">nrow</span><span class="p">(</span>counts<span class="p">))</span>
-
-unblocked.dge <span class="o">&lt;-</span> estimateGLMCommonDisp<span class="p">(</span>dge<span class="p">,</span> unblocked.design<span class="p">,</span> verbose<span class="o">=</span><span class="kc">TRUE</span><span class="p">)</span>
-unblocked.dge <span class="o">&lt;-</span> estimateGLMTrendedDisp<span class="p">(</span>unblocked.dge<span class="p">,</span> unblocked.design<span class="p">)</span>
-unblocked.dge <span class="o">&lt;-</span> estimateGLMTagwiseDisp<span class="p">(</span>unblocked.dge<span class="p">,</span> unblocked.design<span class="p">)</span>
-unblocked.fit <span class="o">&lt;-</span> glmFit<span class="p">(</span>unblocked.dge<span class="p">,</span> unblocked.design<span class="p">)</span>
-unblocked.lrt <span class="o">&lt;-</span> glmLRT<span class="p">(</span>unblocked.dge<span class="p">,</span> unblocked.fit<span class="p">,</span> coef<span class="o">=</span><span class="s">&quot;ConditionIFNg activated&quot;</span><span class="p">)</span>
-unblocked.tt <span class="o">&lt;-</span> topTags<span class="p">(</span>unblocked.lrt<span class="p">,</span> n<span class="o">=</span><span class="kp">nrow</span><span class="p">(</span>counts<span class="p">))</span>
-
-write.csv<span class="p">(</span>blocked.tt<span class="o">$</span><span class="kp">table</span><span class="p">,</span> <span class="s">&quot;edgeR-genes-prelim&quot;</span><span class="p">)</span>
-write.xlsx.multisheet<span class="p">(</span><span class="kt">list</span><span class="p">(</span>blocked<span class="o">=</span>blocked.tt<span class="o">$</span><span class="kp">table</span><span class="p">,</span>
-                           unblocked<span class="o">=</span>unblocked.tt<span class="o">$</span><span class="kp">table</span><span class="p">),</span>
-                      <span class="s">&quot;edgeR-genes.xlsx&quot;</span><span class="p">)</span>
-</pre></div>
-</body>
-</html>