In aayushraman/DASC: Detecting hidden batch factors through data adaptive adjustment for biological effects

knitr::opts_chunk$set(tidy=FALSE, cache=TRUE, 
                        dev="png", message=FALSE, error=FALSE, warning=TRUE)
## Libraries
require("Biobase")
require("NMF")
require("cvxclustr")

Getting started

r Biocpkg("DASC") is an R package distributed as part of the Bioconductor project. To install the package, start R and enter:

source("http://bioconductor.org/biocLite.R")
biocLite("DASC")

Introduction

r Biocpkg("DASC") is used for identifying batches and classifying samples into different batches in a high dimensional gene expression dataset. The batch information can be further used as a covariate in conjunction with other variables of interest among standard bioinformatics analysis like differential expression analysis.

Citation info

If you use r Biocpkg("DASC") for your analysis, please cite it as here below. To cite package ‘DASC’ in publications use:

    @Manual{, 
        title = {DASC: Detecting hidden batch factors through data adaptive 
            adjustment for biological effects.},
        author = {Haidong Yi, Ayush T. Raman, Han Zhang, Genevera I. Allen and 
            Zhandong Liu},
        year = {2017},
        note = {R package version 0.1.0},
    }

Quick Example

library(DASC)
data("esGolub")
samples <- c(20,21,28,30)
dat <- exprs(esGolub)[1:100,samples]
pdat <- pData(esGolub)[samples,]

## use nrun = 50 or more for better convergence of results
res <- DASC(edata = dat, pdata = pdat, factor = pdat$Cell, method = 'ama',
                type = 3, lambda = 1, rank = 2:3, nrun = 5,
                annotation="esGolub Dataset")
#consensusmap(res)
#plot(res)

Setting up the data

The first step in using DASC package is to properly format the data. For example, in case of gene expression data, it should be a matrix with features (genes, transcripts) in the rows and samples in the columns. DASC then requires the information for the variable of interest to model the gene expression data effectively.Variable of interest could be a genotype or treatment information.

Stanford RNA-Seq Dataset

Below is an example of Stanford gene expression dataset (Chen et. al. PNAS, 2015; Gilad et. al. F1000 Research, 2015). It is a filtered raw counts dataset which was published by Gilad et al. F1000 Research. 30% of genes with the lowest expression & mitochondrial genes were removed (Gilad et al.F1000 Research).

## libraries
set.seed(99999)
library(DESeq2)
library(ggplot2)
library(pcaExplorer)

## dataset
rawCounts <- stanfordData$rawCounts
metadata <- stanfordData$metadata

## Using a smaller dataset
idx <- which(metadata$tissue %in% c("adipose", "adrenal", "sigmoid"))
rawCounts <- rawCounts[,idx]
metadata <- metadata[idx,]

head(rawCounts)
head(metadata)

Batch detection using PCA Analysis

## Normalizing the dataset using DESeq2
dds <- DESeqDataSetFromMatrix(rawCounts, metadata, design = ~ species+tissue)
dds <- estimateSizeFactors(dds)
dat <- counts(dds, normalized = TRUE)
lognormalizedCounts <- log2(dat + 1)

## PCA plot using 
rld.dds <- rlog(dds)
pcaplot(rld.dds, intgroup=c("tissue","species"), ntop=1000, pcX=1, pcY=2)

In the PCA plot, PC1 shows the differences between the species. PC2 shows the differences between the species i.e. samples clustering based on tissues.

Batch detection using DASC

res <- DASC(edata = dat, pdata = metadata, factor = metadata$tissue,
                method = 'ama', type = 3, lambda = 1, rank = 2:3, nrun = 10,
                annotation = 'Stanford Dataset')

## Consensus plot
consensusmap(res)

## Residual plot
plot(res)

## Batches -- dataset has 6 batches
sample.clust <- data.frame(sample.name = colnames(lognormalizedCounts), 
                            clust = as.vector(predict(res$fit$`2`)), 
                            batch = metadata$seqBatch)
ggplot(data = sample.clust, aes(x=c(1:6), y=clust, color=factor(clust))) + 
    geom_point(size = 4) + xlab("Sample Number") + ylab("Cluster Number")

Based on the above plots, we observe that the dataset has 2 batches. This can further be compared with the sequencing platform or metadata$seqBatch. The results suggest that differences in platform led to batch effects. Batch number can be used as another covariate, when differential expression analyses using DESeq2,edgeR or limma are performed.