|
| 1 | +--- |
| 2 | +title: "Getting Started with Robseq" |
| 3 | +output: rmarkdown::html_vignette |
| 4 | +vignette: > |
| 5 | + %\VignetteIndexEntry{Getting Started with Robseq} |
| 6 | + %\VignetteEngine{knitr::rmarkdown} |
| 7 | + %\VignetteEncoding{UTF-8} |
| 8 | +--- |
| 9 | + |
| 10 | +```{r setup, include=FALSE} |
| 11 | +knitr::opts_chunk$set(echo = TRUE) |
| 12 | +``` |
| 13 | + |
| 14 | +# Importing libraries |
| 15 | +```{r eval=FALSE} |
| 16 | +library(Robseq) |
| 17 | +library(edgeR) |
| 18 | +library(doParallel) |
| 19 | +library(EnhancedVolcano) |
| 20 | +``` |
| 21 | + |
| 22 | +# Loading Example Data |
| 23 | +Loading Colon cancer data |
| 24 | +```{r eval=FALSE} |
| 25 | +load("~/Current Data Path/Colon Cancer.RData") |
| 26 | +``` |
| 27 | + |
| 28 | +Extracting Colon cancer gene expression data and metadata |
| 29 | +```{r eval=FALSE} |
| 30 | +features = data$counts |
| 31 | +metadata = data$metadata |
| 32 | +``` |
| 33 | + |
| 34 | +# Snapshot of data |
| 35 | +A typical bulk RNA-seq gene expression count data frame looks something like below. Note, the genes should be in the rows and the samples in columns |
| 36 | +```{r eval=FALSE} |
| 37 | +features[1:5, 1:5] |
| 38 | +``` |
| 39 | + |
| 40 | +``` |
| 41 | +## GSM4731674 GSM4731675 GSM4731676 GSM4731677 GSM4731678 |
| 42 | +## TSPAN6 103 44 76 417 630 |
| 43 | +## TNMD 1 0 0 0 18 |
| 44 | +## DPM1 86 63 97 173 309 |
| 45 | +## SCYL3 32 40 77 56 97 |
| 46 | +## C1orf112 28 35 25 60 55 |
| 47 | +``` |
| 48 | + |
| 49 | +A typical metadata data frame looks something like below. Note, in our pipeline the user must include a column labeled as "Exposure" which should be the variable which will be used by Robseq in performing differential expression analysis. If the user has a different label for the treatment/condition or disease status variable then the user should supply that name via "expVar" argument in Robseq |
| 50 | +```{r eval=FALSE} |
| 51 | +metadata[1:5, ] |
| 52 | +``` |
| 53 | + |
| 54 | +``` |
| 55 | +## Sample Exposure |
| 56 | +## 1 GSM4731674 Tumor |
| 57 | +## 2 GSM4731675 Tumor |
| 58 | +## 3 GSM4731676 Tumor |
| 59 | +## 4 GSM4731677 Tumor |
| 60 | +## 5 GSM4731678 Tumor |
| 61 | +``` |
| 62 | + |
| 63 | +# Preprocessing |
| 64 | +A typical preprocessing step is to filter lowly abundant genes. We do so using edgeR's "filterByExpr" function |
| 65 | +```{r eval=FALSE} |
| 66 | +keep.exprs <- filterByExpr(features, group = as.factor(metadata$Exposure)) |
| 67 | +paste(length(which(!keep.exprs)), ' lowly expressed genes were filtered out', sep = '') |
| 68 | +features <- features[keep.exprs, ] |
| 69 | +``` |
| 70 | + |
| 71 | +``` |
| 72 | +## [1] "18625 lowly expressed genes were filtered out" |
| 73 | +``` |
| 74 | + |
| 75 | +# Performing differential expression analysis using Robseq |
| 76 | +To perform differential gene expression (DGE) analyses using Robseq please use the code below. Note, if your metadata has only the "Exposure" variable (such as treatment groups, conditions, disease status, and etc.) set the "coVars" argument to "NULL". However, in this working example our metadata contained three covariates such as "post-mortem interval (pmi)", "rna integrity number (rin)" and "age at death" which needed to be adjusted for in our DGE analysis. Therefore, we have supplied this three variables to the "coVars" argument |
| 77 | +```{r eval=FALSE} |
| 78 | +fit <- Robseq::robust.dge(features = features, |
| 79 | + metadata = metadata, |
| 80 | + norm.method = "RLE", |
| 81 | + expVar = "Exposure", |
| 82 | + coVars = NULL, |
| 83 | + filter = FALSE, |
| 84 | + parallel = TRUE, |
| 85 | + ncores = detectCores() - 2, |
| 86 | + verbose = FALSE) |
| 87 | +``` |
| 88 | + |
| 89 | +# Obtaining results from Robseq |
| 90 | +After performing DGE analysis you can extract the results table in the following manner |
| 91 | +```{r eval=FALSE} |
| 92 | +results <- fit$res |
| 93 | +``` |
| 94 | + |
| 95 | +The results table from Robseq should look something like the following |
| 96 | +```{r eval=FALSE} |
| 97 | +results[1:5,] |
| 98 | +``` |
| 99 | + |
| 100 | +``` |
| 101 | +## Genes log2FC SE L.CI U.CI Pval adjPval |
| 102 | +## 6872 BEST4 -6.631 0.255 -6.131209 -7.130791 1.565312e-54 1.617267e-50 |
| 103 | +## 10982 ETV4 5.184 0.202 5.579913 4.788087 2.121424e-54 1.617267e-50 |
| 104 | +## 11710 OTOP3 -6.227 0.244 -5.748769 -6.705231 1.430501e-53 7.270284e-50 |
| 105 | +## 11725 OTOP2 -7.246 0.298 -6.661931 -7.830069 2.196497e-51 8.372498e-48 |
| 106 | +## 14866 SPIB -5.385 0.239 -4.916569 -5.853431 7.877754e-48 2.402242e-44 |
| 107 | +``` |
| 108 | + |
| 109 | +# Volcano plot to visualize DGE results |
| 110 | +Volcano plot to visualize the DGE results obtained from Robseq |
| 111 | +```{r eval=FALSE} |
| 112 | + EnhancedVolcano(results, |
| 113 | + lab = results$Genes, |
| 114 | + x = 'log2FC', |
| 115 | + y = 'adjPval') |
| 116 | +``` |
| 117 | + |
| 118 | + |
0 commit comments