## ----style-knitr, eval=TRUE, echo=FALSE, results="asis"------------------
BiocStyle::latex(use.unsrturl=FALSE)

## ----setup_knitr, include=FALSE, cache=FALSE-----------------------------
library(knitr)
opts_chunk$set(cache = FALSE, tidy = FALSE, tidy.opts = list(blank = FALSE), 
  highlight=FALSE, out.width = "7cm", out.height = "7cm", fig.align = "center")

## ----DSpasilla1----------------------------------------------------------
library(PasillaTranscriptExpr)

data_dir  <- system.file("extdata", package = "PasillaTranscriptExpr")

metadata <- read.table(file.path(data_dir, "metadata.txt"), header = TRUE,
  as.is = TRUE)
metadata

counts <- read.table(file.path(data_dir, "counts.txt"), header = TRUE,
  as.is = TRUE)
head(counts)


## ----DSlibrary, message=FALSE--------------------------------------------
library(DRIMSeq)

## ----DSdmDSdata_create---------------------------------------------------
d <- dmDSdata(counts = counts[, metadata$SampleName], gene_id = counts$gene_id,
  feature_id = counts$feature_id, sample_id = metadata$SampleName,
  group = metadata$condition)
d
head(counts(d), 3)
head(samples(d), 3)

## ----DSmax_nr_isoforms, echo = FALSE-------------------------------------
max_nr_isoforms <- max(elementNROWS(dm_counts(d)))

## ----DSdmDSdata_plot-----------------------------------------------------
plotData(d)

## ----DSdmDSdata_subset---------------------------------------------------
gene_id_subset <- readLines(file.path(data_dir, "gene_id_subset.txt"))
d <- d[names(d) %in% gene_id_subset, ]
d
plotData(d)

## ----DSlength_d, echo = FALSE--------------------------------------------
length_d <- length(d)

## ----DSdmFilter----------------------------------------------------------
# Check what is the minimal number of replicates per condition
table(samples(d)$group)

d <- dmFilter(d, min_samps_gene_expr = 7, min_samps_feature_expr = 3,
  min_samps_feature_prop = 0)
plotData(d)

## ----DSdmDispersion------------------------------------------------------
d <- dmDispersion(d, verbose = 1, BPPARAM = BiocParallel::SerialParam())
d
head(mean_expression(d), 3)
common_dispersion(d)
head(genewise_dispersion(d), 3)

## ----DSdmDispersion_plot-------------------------------------------------
library(ggplot2)
ggp <- plotDispersion(d)
ggp + geom_point(size = 4)

## ----DSdmFit, out.width = "16cm", fig.width = 16-------------------------
d <- dmFit(d, BPPARAM = BiocParallel::SerialParam())
d

head(proportions(d), 3)
head(statistics(d), 3)

## ----DSdmTest------------------------------------------------------------
d <- dmTest(d, verbose = 1, BPPARAM = BiocParallel::SerialParam())
d

head(proportions(d), 3)
head(statistics(d), 3)


## ----DSdmTest_results----------------------------------------------------
head(results(d), 3)

## ----DSdmTest_plot-------------------------------------------------------
plotTest(d)

## ----DSdmLRT_plotFit, out.width = "16cm", fig.width = 16-----------------
res <- results(d)
res <- res[order(res$pvalue, decreasing = FALSE), ]
gene_id <- res$gene_id[1]

plotFit(d, gene_id = gene_id)
plotFit(d, gene_id = gene_id, plot_type = "lineplot")
plotFit(d, gene_id = gene_id, plot_type = "ribbonplot")


## ----SQTLgeuvadis, message=FALSE-----------------------------------------
library(GeuvadisTranscriptExpr)

counts <- GeuvadisTranscriptExpr::counts
genotypes <- GeuvadisTranscriptExpr::genotypes
gene_ranges <- GeuvadisTranscriptExpr::gene_ranges
snp_ranges <- GeuvadisTranscriptExpr::snp_ranges


## ----SQTLlibrary, message=FALSE------------------------------------------
library(DRIMSeq)

## ----SQTLdmSQTLdata_create, message=FALSE--------------------------------
# Make sure that samples in counts and genotypes are in the same order
sample_id <- colnames(counts[, -(1:2)])

d <- dmSQTLdataFromRanges(counts = counts[, sample_id],
  gene_id = counts$Gene_Symbol, feature_id = counts$TargetID,
  gene_ranges = gene_ranges, genotypes = genotypes[, sample_id],
  snp_id = genotypes$snpId, snp_ranges = snp_ranges, sample_id = sample_id,
  window = 5e3, BPPARAM = BiocParallel::SerialParam())
d

## ----multiplot, echo = FALSE---------------------------------------------
# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
      ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
        layout.pos.col = matchidx$col))
    }
  }
}

## ----SQTLdmSQTLdata_plot, out.width = "16cm", fig.width = 16-------------
multiplot(plotlist = plotData(d), cols = 3)

## ----SQTLdmFilter, out.width = "16cm", fig.width = 16--------------------
d <- dmFilter(d, min_samps_gene_expr = 70, min_samps_feature_expr = 5,
  min_samps_feature_prop = 0, minor_allele_freq = 5,
  BPPARAM = BiocParallel::SerialParam())

multiplot(plotlist = plotData(d), cols = 3)

## ----SQTLdmDispersion----------------------------------------------------
d <- dmDispersion(d, verbose = 1, BPPARAM = BiocParallel::SerialParam())
d

## ----SQTLplotDispersion--------------------------------------------------
library(ggplot2)
ggp <- plotDispersion(d)
ggp + geom_point(size = 4)

## ----SQTLdmFit-----------------------------------------------------------
d <- dmFit(d, BPPARAM = BiocParallel::SerialParam())
d

## ----SQTLdmTest----------------------------------------------------------
d <- dmTest(d, verbose = 1, BPPARAM = BiocParallel::SerialParam())
d
head(results(d), 3)
plotTest(d)

## ----SQTLplotFit, out.width = "16cm", fig.width = 16, warning = FALSE----
res <- results(d)
res <- res[order(res$pvalue, decreasing = FALSE), ]

gene_id <- res$gene_id[1]
snp_id <- res$snp_id[1]

plotFit(d, gene_id, snp_id)
plotFit(d, gene_id, snp_id, plot_type = "boxplot2", order = FALSE)
plotFit(d, gene_id, snp_id, plot_type = "ribbonplot")

## ----sessionInfo---------------------------------------------------------
sessionInfo()