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.0094728 8.524309e-03 2.5 WT>MZsox
## 2 crtap -2.5849625 1.273661e-01 2.5 WT=MZsox
## 3 rpl3 -0.1496811 4.931881e-01 2.5 WT=MZsox
## 4 zic2b -0.5252563 3.603394e-01 2.5 WT=MZsox
## 5 ccn1 -2.5004662 1.240784e-08 2.5 WT<MZsox
## 6 vipas39 0.9528740 1.621983e-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/Rtmpd7Y1rG/Rbuild1de5972f2d988e/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.