%\VignetteIndexEntry{NanoStringDiff Vignette} %\VignettePackage{NanoStringDiff} %\VignetteKeyword{NanoString} %\VignetteKeyword{Differential Expression} \documentclass[12pt]{article} \usepackage{float} \usepackage{Sweave} \usepackage{amssymb} \RequirePackage{Bioconductor} \AtBeginDocument{\bibliographystyle{unsrturl}} \renewcommand{\baselinestretch}{1.3} \SweaveOpts{keep.source=TRUE,eps=FALSE,include=TRUE,width=4,height=4} \author{Hong Wang$^{1}$, Chi Wang$^{2,3}$\footnote{to whom correspondence should be addressed} \\[1em] \small{$^{1}$Department of Statistics , University of Kentucky,Lexington, KY;}\\ \small{$^{2}$Markey Cancer Center, University of Kentucky, Lexington, KY ;}\\ \small{$^{3}$Department of Biostatistics, University of Kentucky, Lexington, KY;}\\ \small{\texttt{hong.wang@uky.edu}}} \title{\textsf{\textbf{The NanoStringDiff package}}} %\bibliographystyle{abbrv} \begin{document} \SweaveOpts{concordance=TRUE} \maketitle \begin{abstract} This vignette introduces the use of the Bioconductor package NanoStringDiff, which is designed for differential analysis based on NanoString nCouner Data. NanoStringDiff considers a generalized linear model of the negative binomial family to characterize count data and allows for multi-factor design. Data normalization is incorporated in the model framework by including data normalization parameters estimated from positive controls, negative controls and housekeeping genes embedded in the nCounter system. Present method propose an empirical Bayes shrinkage approach to estimate the dispersion parameter and a likelihood ratio test to identify differential expression genes. \end{abstract} \newpage \tableofcontents \newpage \section{Citation} The package {\tt NanoStringDiff} implements statistical methods from the following publication. If you use {\tt NanoStringDiff} in the published research, please cite: \\ Hong Wang, Craig Horbinski, Hao Wu, Yinxing Liu,Arnold J. Stromberg and Chi Wang: A Negative Binomial Model-Based Method for Differential Expression Analysis Based on NanoString nCounter Data.(Manuscript) \section{Quick Start} This section show the most basic steps for a differential expression analysis for NanoString nCounter Data: \begin{enumerate} \item Create a {\tt NanoStringSet} object using {\tt createNanoStingSet} or {\tt createNanoStringSetFromCsv }(see examples in the Data Input section). In this section we use {\tt NanoStringData} directly, which is an object of {\tt NanoStringSet} contained in the package. \item Make a design matrix to describe treatment conditions. \item Estimate norlamization factors including positive size factors, housekeeping size factors and background noise using {\tt estNormalizationFactors} \item Perform a likelihood ratio test using {\tt glm.LRT}. \end{enumerate} <>= library("Biobase") library("NanoStringDiff") data(NanoStringData) pheno=pData(NanoStringData) design.full=model.matrix(~0+pheno$group) NanoStringData=estNormalizationFactors(NanoStringData) result=glm.LRT(NanoStringData,design.full,contrast=c(1,-1)) @ Here, the data NanoStringData contained in the package is an animal data, we called {\tt MoriData}\cite{mori2014hippo}. Mori et al tried to study the possible reasons responsible for the widespread miRNA repression observed in cancer, global microRNA expression in mouse liver normal tissues and liver tumors induced by deletion of Nf2 (merlin) profiled by nCounter Mouse miRNA Expression Assays. Expressions of 599 miRNAs were profiled for two replicates in each group. \section{Data Input} NanoStringDiff works on matrix of integer read counts from NanoString nCounter analyzer with endogenes, positive control genes, negative control genes and housekeeping control genes. For the matrix, rows corresponding to genes and columns to independent samples or replicates. The counts represent the total number of reads aligning to each gene (or other genomic locus). The count values must be raw counts of NanoString nCounter data, since data normalization is incorporated in the model framework by including data normalization parameters estimated from positive controls, negative controls and housekeeping genes using function {\tt estNormalizationFactors}. Hence, please do not supply normalized counts. There must be have six positive control genes order by different concentrations from high to low, since NanoString nCounter analyzer provide six different samll positive controls with six different concentrations in the 30 uL hybridization: 128 fM, 32 fM, 8 fM, 2 fM, 0.5 fM, and 0.125 fM. No such restriction for negative control genes and housekeeping control genes.Nanostring recommends at least three housekeeping genes, but the more that are included, the more accurate the normalization will be. \subsection{Create NanoStringSet from csv.file} The data produced by the nCounter Digital Analyzer are exported as a Reporter Code Count (RCC) file. RCC files are comma-separated text(.csv) files that contain the counts for each gene in a sample and the data for each sample hybridization are contained in a separate RCC file. So before you call {\tt createNanoStringSetFromCsv} to creat a {\tt NanoStringSet} object, you should create a csv.file to combine all interesting samples together, and the first three columns shoud be "CodeClass", "Name" ane "Accession", the counts data contained from the $4^{th}$ column to the last column. Note: \begin{enumerate} \item The $1^{st}$ column "CodeClass" should specify the function of genes as "Positive", "Negative","Housekeeping" or "Endogenous". \item Some data set have "Spikein" genes, you need delete these "spikein" genes or you could just leave there if you use {\tt createNanoStringSetFromCsv} to creat {\tt NanoStringSet} object(See example in the Data Input section). But NanoStringDiff only works with "positive", "negative", housekeeping" and "endogenous". \end{enumerate} The "csv.file" should looks like as following Figure: \begin{figure}[ht] \centering \includegraphics{example.PNG} \caption{Example of csv.file pattern} \label{example} \end{figure} When you created a csv.file, you will specify a variable which points on the directory in which your csv.file is located <>= directory <- "/path/to/your/files/" @ However, for demonstration purposes only, the following line of code points to the directory for the "Mori.csv" in the NanoStringDiff package. <>= directory <- system.file("extdata", package="NanoStringDiff", mustWork=TRUE) path<-paste(directory,"Mori.csv",sep="/") @ The phenotype informations of the data should be stored as data frame. <>= designs=data.frame(group=c("Normal","Normal","Tumor","Tumor")) designs @ <>= library("NanoStringDiff") NanoStringData=createNanoStringSetFromCsv(path,header=TRUE,designs) NanoStringData @ \subsection{Create NanoStringSet from matrix} If you already read your positive control genes, negative control genes, housekeeping control genes and endogous into R session separately and stored as matrix, then you can use {\tt createNanoStringSet} to create a {\tt NanoStringSet} object. <>= endogenous=matrix(rpois(100,50),25,4) colnames(endogenous)=paste("Sample",1:4) positive=matrix(rpois(24,c(128,32,8,2,0.5,0.125)*80),6,4) colnames(positive)=paste("Sample",1:4) negative=matrix(rpois(32,10),8,4) colnames(negative)=paste("Sample",1:4) housekeeping=matrix(rpois(12,100),3,4) colnames(housekeeping)=paste("Sample",1:4) designs=data.frame(group=c("Control","Control","Treatment","Treatment"), gender=c("Male","Female","Female","Male"), age=c(20,40,39,37)) NanoStringData1=createNanoStringSet(endogenous,positive, negative,housekeeping,designs) NanoStringData1 pData(NanoStringData1) head(exprs(NanoStringData1)) @ \section{Differential Expression Analysis for Single Factor Experiment} For general experiments, once normalization factors obtained using {\tt estNormalizationFactors}, given a design matrix and contrast, we can proceed with testing procedures for determing differential expression using the generalized linear model (GLM) likelihood ratio test. The testing can be done by using the function {\tt glm.LRT} and return a list with a component is table including: logFC ,lr, pvalue and qvalue(adjust p value using the procedure of Benjamini and Hochberg). \subsection{Two Group Comparisons} The simplest and most common type of experimental design is two group comparison,like treatment group vs control group. As a brief example, consider a simple situation with control group and treatment group, each with two replicates, and the researcher wants to make comparisons between them. Here, we still use {\tt NanoSreingData1} we created above to demonstrate this example. Make design matrix using {\tt model.matrix}: <>= pheno=pData(NanoStringData1) group=pheno$group design.full=model.matrix(~0+group) design.full @ Note that there is no intercept column in the dasign matrix, each column is for each group, since we include {\tt 0+} in the model formula. If the researcher want compare Treatment to Control, that is Treatment- Control, the contrast vector is try to the comparison -1*Control+1*Treatment, so the contrast is : <>= contrast=c(-1,1) @ Estimate normalization factors and return the same object with positiveFactor, negativeFactor and housekeepingFactor slots filled or replaced: <>= NanoStringData1=estNormalizationFactors(NanoStringData1) positiveFactor(NanoStringData1) negativeFactor(NanoStringData1) housekeepingFactor(NanoStringData1) @ Generalize linear model likelihood ratio test: <>= result=glm.LRT(NanoStringData1,design.full,contrast=contrast) head(result$table) str(result) @ Note that the text Treatment compare to Control tells you that the estimates of logFC is log2(Treatment/Control). \subsection{Pairwise Comparisons} Now consider a researcher has three treatment groups say A, B and C, and researcher wants to compare each groups like: B compare to A, C compare to A, and B compare to C. First create an object of {\tt NanoStringSet} with pseudo data: <>= endogenous=matrix(rpois(300,50),25,12) colnames(endogenous)=paste("Sample", 1:12) colnames(endogenous)=paste("Sample",1:12) positive=matrix(rpois(72,c(128,32,8,2,0.5,0.125)*80),6,12) negative=matrix(rpois(96,10),8,12) housekeeping=matrix(rpois(36,100),3,12) designs=data.frame(group=c(rep("A",4),rep("B",4),rep("C",4)), gender=rep(c("Male","Male","Female","Female"),3), age=c(20,40,39,37,29,47,23,45,34,65,35,64)) NanoStringData2=createNanoStringSet(endogenous,positive, negative,housekeeping,designs) NanoStringData2 pData(NanoStringData2) @ Make design matrix only consider one experimential factor {\tt group}: <>= pheno=pData(NanoStringData2) group=pheno$group design.full=model.matrix(~0+group) design.full @ Estimate normalization factors: <>= NanoStringData2=estNormalizationFactors(NanoStringData2) @ B compare to A, the contrast should be (-1,1,0) <>= result=glm.LRT(NanoStringData2,design.full,contrast=c(-1,1,0)) @ C compare to A, the contrast should be (-1,0,1) <>= result=glm.LRT(NanoStringData2,design.full,contrast=c(-1,0,1)) @ B compare to c, the contrast should be (0,1,-1) <>= result=glm.LRT(NanoStringData2,design.full,contrast=c(0,1,-1)) @ The other way to creat a design matrix in R is to include an intercept term to represents the base level of the factor. We just omitte the {\tt 0+} in the model formula when we create design matrix using model.matrix function: <>= design.full=model.matrix(~group) design.full @ In this sitution, the first coefficient measure the log scaler of baseline mean expression level in group A, and the second and third column are now relative to the baseline, not represent the mean expression level in group B and C. So, if we want to compare B to A, now that is equivalent to test the $2^{nd}$ coefficient equal to 0. <>= result=glm.LRT(NanoStringData2,design.full,Beta=2) @ Beta=3 means compare C to A: <>= result=glm.LRT(NanoStringData2,design.full,Beta=3) @ \subsection{Multigroup Comparisons} NanoStringDiff can do multigroup comparisons, for example if we want compare group A to the average of group B and C, <>= design.full=model.matrix(~0+group) result=glm.LRT(NanoStringData2,design.full,contrast=c(1,-1/2,-1/2)) @ \section{Differential Expression Analysis for Multifactor Experiment} In this section, we still use {\tt NanoStringData2} in examples, but the above examples all cotnsider single factor {\tt treatment condition} , now we include the other experiment factor {\tt gender} to consider multifactor problems. \subsection{Nested Interactions} First, we form the design matrix use factors {\tt group} and interaction between {\tt group} and {\tt gender}: <>= design=pData(NanoStringData2)[,c("group","gender")] group=design$group gender=design$gender design.full=model.matrix(~group+group:gender) design.full @ Here, we consider the mean expression level of female in group A as the baseline expression level, if we want to test the effect of gender in group A, we can use Beta=4, <>= result=glm.LRT(NanoStringData2,design.full,Beta=4) @ Compare treatment group B to A ignore the effect the gender, <>= result=glm.LRT(NanoStringData2,design.full,Beta=2) @ Consider the interaction effect between gender and group B compare to A: <>= result=glm.LRT(NanoStringData2,design.full,Beta=c(2,5)) @ \subsection{All Interactions} We also can form a design matrix consider all interactions: <>= design.full=model.matrix(~group+gender+group:gender) @ Which is equivalent to: <>= design.full=model.matrix(~group*gender) design.full @ Test the gender effect in Treatment group B: <>= result=glm.LRT(NanoStringData2,design.full,Beta=4:5) @ \section{Session Info} <>= toLatex(sessionInfo()) @ \bibliography{reference} \end{document}