2 2. Analysis of an RNA-Seq Dataset
Here, an RNA-Seq dataset of Bacillus anthracis (ExpressionSet Baca) containing different stress conditions as sample covariates and COG annotations (Clusters of Orthologous Groups, Tatusov et al., 2000) as gene covariates will be analysed.
2.1 2.1 Preparation of the dataset
We load the covRNA package and the integrated Baca dataset, which contains the ExpressionSet Baca. The assayData contains deep sequenced RNA-Seq data of B. anthracis under four stress conditions (with four replicates per stress conditions). The raw sequence reads derive from Passalacqua et al. (2012) and are availaible at Gene Expression Omnibus (GEO, accession number GSE36506). We have already mapped, counted and DESeq2 (Love et al., 2014) normalised these counts. The phenoData assigns the stress condition, i.e. ctrl, cold, salt and alcohol stress, to the samples. The featureData contains COG annotations of the genes.
library(covRNA)## Loading required package: ade4## Loading required package: Biobase## Loading required package: BiocGenerics## Loading required package: parallel##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB## The following object is masked from 'package:ade4':
##
## score## The following objects are masked from 'package:stats':
##
## IQR, mad, xtabs## The following objects are masked from 'package:base':
##
## Filter, Find, Map, Position, Reduce, anyDuplicated, append,
## as.data.frame, cbind, colnames, do.call, duplicated, eval,
## evalq, get, grep, grepl, intersect, is.unsorted, lapply,
## lengths, mapply, match, mget, order, paste, pmax, pmax.int,
## pmin, pmin.int, rank, rbind, rownames, sapply, setdiff, sort,
## table, tapply, union, unique, unsplit, which, which.max,
## which.min## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.data(Baca)2.2 2.2 Fourthcorner Analysis with stat
We use the function stat to statistically analyse associations between gene and sample covariates.
statBaca <- stat(ExprSet = Baca, npermut = 999, padjust = "BH", nrcor = 2, exprvar = 1)statBaca is then an object of type stat. As a list, it saves all results as well as the input of the function. For instance, we can access the adjusted p-values of all covariate combinations.
ls(statBaca)
adjp <- statBaca$adj.pvalue; adjpTo visualise the results, the stat object can be plotted (Figure 1). If the plot shall be shown in high quality, we advise to use the default setting pdf=TRUE.
plot(statBaca, xnames = c('cold','ctrl','etoh','salt'), shiftx = -0.1)
Figure 1: Cross-tabulation of the adjusted p-values between the sample covariates and the gene functions of the ExpressionSet Baca. Red and blue cells show negative and positive significant associations at \(\alpha\)=0.05, respectively.
The cross-tabulation of the sample covariates and the gene functions visualises negative and positive significant associations at \(\alpha\)=0.05. Five significant associations can be discovered.
2.3 2.3 RLQ with ord
We use the function ord to visualise sample and gene covariates in one coordinate system.
ordBaca <- ord(Baca)ordBaca is then an object of type ord. Different features of this object can be plotted by using the feature argument of the plot function (see manpages for more information). For instance, we can plot the amount of variance explained by each axis (Figure 2).
plot(ordBaca, feature = "variance")
Figure 2: Barplot of the amount of variance explained by each axis of the ordination ordBaca. As the number of axes to be taken into account by ordination is by default set to 2, the bars of the first two axes are shown in black.
The first two axes of the RLQ explain a large amount of the variance of the data (93.81\(\%\) and 4.09\(\%\), respectively).
2.4 2.4 Combination of Results
The results of the functions stat and ord can be simultaneously visualised by the function vis (Figure 3).
vis(Stat = statBaca, Ord = ordBaca, rangex=4, rangey=4)
Figure 3: Simultaneous visualisation of the statistical analysis and of the ordination of Baca. Blue and red lines between the covariates represent positive and negative significant associations, respectively.
Here, covariates involved in at least one significant association are shown in black, others are shown in gray. All significant covariates are connected by lines which colour represents the character of their association. As expected, positively associated covariates are situated close to each other and at similar angles from the origin. On the contrary, negatively associated covariates are distant from each other.
We can further observe that the first axis seems to be spanned by the difference between the classes J and O. The second axis seems to be spanned between ctrl and the other treatments. Spatial proximity of cold and salt treatment in the second quadrant suggest that they have similar functional effects on the gene expression.
2.5 2.5 Comparison with Other Methods
To validate our results of the analysis of Baca, we compare them to traditional approaches like hypergeometric test (HG), Mann-Whitney rank test (RANK) und gene set enrichment analysis (GSEA, Subramanian et al., 2005) by using the R package BOG (see Park et al., 2015, for further details).
RANK and GSEA discover class J as significantly enriched (p=6.40e-11 and p=0.02, respectively). HG does not detect any significant gene functions.