Step-by-step Tutorial
First of all, data has to be provided in the right format. We made
use of the S4-class SummarizedExperiment that gives the
possibility to provide additional information for example on covariates
for the measurements that themselves are given as a numeric matrix. Here
is an example how your data should look like:
## class: SummarizedExperiment
## dim: 15775 40
## metadata(0):
## assays(1): mydata
## rownames: NULL
## rowData names(2): name genename
## colnames(40): WT_2.5_B1 WT_3_B1 ... MZsox_5.5_B1 MZsox_6_B1
## colData names(3): condition time replicateSummarizedExperiment provides a constructor to easily
bring your data into this format. When constructing your object, make
sure to provide genenames in the rowData argument and
information on condition, time point and, if available, replicate
identifier, in the colData argument.
Next, the conditions that should be analyzed are specified and a
threshold is provided which is used to exclude genes with expression
levels below this threshold for all conditions. This can be useful, if
expression levels are in the range of your detection limit. Since
RNAseq makes use of the parallel package, you
may specify a number of cores in order to speed up your computation
time.
analyzeConditions <- c("WT", "MZsox")
thCount <- 100
nrcores <- 1
library(SummarizedExperiment)
#if(Sys.info()[[1]]=="Windows"){nrcores <- 1} # use parallelization only on Linux and Mac
mydata <- mydata[seq(1,nrow(mydata), by=4),]
vec2Keep <- which(vapply(1:dim(mydata)[1],function(i)
!Reduce("&",assays(mydata)[[1]][i,]<thCount), c(TRUE)))
mydata <- mydata[vec2Keep,] # threshold is applied
times <- unique(sort(as.numeric(colData(mydata)$time))) # get measurement times from input dataAfter data preparation, fold change detection can be performed. The
function getFC internally calles functions from the
NBPSeq package to perform fold change analysis for each
gene and at each time point. The result is saved in a
data.frame with corresponding p-values.
resultFC <- getFC(dataset = mydata,
myanalyzeConditions = analyzeConditions,
cores = nrcores,
mytimes = times)
head(resultFC)## name logFoldChange pValue time FCdetect
## 1 enc3 1.0047255 8.726364e-03 2.5 WT>MZsox
## 2 crtap -2.1699250 2.137050e-01 2.5 WT=MZsox
## 3 rpl3 -0.1489668 4.951466e-01 2.5 WT=MZsox
## 4 zic2b -0.5849625 3.171989e-01 2.5 WT=MZsox
## 5 ccn1 -2.5027544 1.206961e-08 2.5 WT<MZsox
## 6 vipas39 0.9485284 1.667527e-02 2.5 WT=MZsoxNote that each gene appears in the data frame as often as the number of time points. The result of the fold change analysis can be visualized as a vulcano plot:
library(ggplot2)
ggplot(subset(resultFC, FCdetect!="none"),
aes(x=logFoldChange, y=-log10(pValue), color=FCdetect)) +
xlab("log2(Fold Change)") + geom_point(shape=20)
Next, the gene expression profiles are analyzed for switches. A
switch appears if the profile shows a statistically significant up- or
downregulation at a specific time point. If up- or downregulation is
detected at multiple time points, the time point with the best
likelihood value (for the one-step model) is chosen. The result is saved
in a data.frame with information on whether a switch has
been detected, at which time point and with which p-value based on the
likelihood ratio test.
resultSwitch <- getSwitch(dataset = mydata,
experimentStepDetection = "WT",
cores = nrcores,
mytimes = times)
head(resultSwitch)## name genename timepoint switch pvalueSwitch experiment
## 7 NM_001001402 enc3 NA none 0.40773375 WT
## 71 NM_001001406 crtap NA none 0.11729573 WT
## 6 NM_001001590 rpl3 NA none 0.10784956 WT
## 5 NM_001001820 zic2b NA none 0.07115325 WT
## 2 NM_001001826 ccn1 NA none 0.16854030 WT
## 4 NM_001001836 vipas39 4 down 0.04349482 WTAfter fold change and switch analysis have been performed, results
shall be collected in one and the same data.frame using the
combineResults function. This function basically prepares
the results to be handed over to plotSSGS.
## name genename timepoint switch pvalueSwitch experiment
## 7 NM_001001402 enc3 NA none 0.40773375 WT
## 71 NM_001001406 crtap NA none 0.11729573 WT
## 6 NM_001001590 rpl3 NA none 0.10784956 WT
## 5 NM_001001820 zic2b NA none 0.07115325 WT
## 2 NM_001001826 ccn1 NA none 0.16854030 WT
## 4 NM_001001836 vipas39 4 down 0.04349482 WT
## FCdown FCup
## 7 2.5hpf 3.0hpf 4.0hpf
## 71 5.0hpf 5.5hpf 6.0hpf
## 6
## 5 3.0hpf 3.5hpf 4.0hpf
## 2 2.5hpf 3.0hpf 3.5hpf 5.0hpf 5.5hpf 6.0hpf
## 4 5.0hpf 5.5hpf 6.0hpfFinally, the plotSSGS function performs Fisher’s exact
test for each combination of fold change time point and switch time
point. Results are plotted as a heat map highlighting combinations with
a high significance.
In order to document the result of the analysis, the function
outputGeneTables provides a possibility to automatically
output switch and fold change information into table files (.txt).
## [1] "Results written to /tmp/RtmpaUe5jT/Rbuild18efad5c2659ba/RNAsense/vignettes"The function outputGeneTables generates five .txt files.
Two of them (geneNamelist) contain gene lists with gene name for genes
that switch up and down respectively. The other two (genelist) contain
exactly the same output but with gene identifiers instead of gene names
depending on what you prefer for further analysis. Each column
corresponds to a combination of switch time point, fold change direction
and time point of fold change. All genes for which fold change was
detected at the indicated time point and switch was detected at the
indicated time point are listed in the corresponding column. Note that a
single gene may appear multiple times. The fifth .txt file (switchList)
contains information on detected switches in a different format. The
output consists of table with six columns with each row corresponding to
one gene. Detected switches are indicated by 1, -1 and 0 for switch up,
switch down and no switch, respectively. If a switch was detected, the
column timepoint indicated the corresponding time point of switch
detection.