4.2 Import data
-We will read in the differential expression results we will download from online. These results are from a zebrafish microarray experiment we used for differential expression analysis for two groups using limma (Ritchie et al. 2015). The table contains Ensembl gene IDs, log fold-changes for each group, and adjusted p-values (FDR in this case). We can identify differentially regulated genes by filtering these results and use this list as input to ORA.
For ORA, we only need a list of gene IDs as our input, so this example can work for any situations where you have gene list and want to know more about what biological pathways it shares genes with.
+For this example, we will read in results from a differential expression analysis that we have already performed. Rather than reading from a local file, we will download the results table directly from a URL. These results are from a zebrafish microarray experiment we used for differential expression analysis for two groups using limma (Ritchie et al. 2015). The table contains Ensembl gene IDs, log fold-changes for each group, and adjusted p-values (FDR in this case). We can identify differentially regulated genes by filtering these results and use this list as input to ORA.
Instead of using the URL below, you can use a file path to a TSV file with your desired gene list results. First we will assign the URL to its own variable called, dge_url.
# Define the url to your differential expression results file
-dge_url <- "https://refinebio-examples.s3.us-east-2.amazonaws.com/02-microarray/results/GSE71270/GSE71270_limma_results.tsv"# Define the url to your differential expression results file
+dge_url <- "https://refinebio-examples.s3.us-east-2.amazonaws.com/02-microarray/results/GSE71270/GSE71270_limma_results.tsv"Read in the file that has differential expression results. Here we are using the URL we set up above, but this can be a local file path instead i.e. you can replace dge_url in the code below with a path to file you have on your computer like: file.path("results", "GSE71270_limma_results.tsv").
# Read in the contents of your differential expression results file
-# `dge_url` can be replaced with a file path to a TSV file with your
-# desired gene list results
-dge_df <- readr::read_tsv(dge_url)## Parsed with column specification:
+# Read in the contents of your differential expression results file
+# `dge_url` can be replaced with a file path to a TSV file with your
+# desired gene list results
+dge_df <- readr::read_tsv(dge_url)
+##
+## ── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────
## cols(
## Gene = col_character(),
## logFC = col_double(),
@@ -3145,7 +4017,7 @@ 4.2 Import data
## B = col_double()
## )
read_tsv() can read TSV files online and doesn’t necessarily require you download the file first. Let’s take a look at what these contrast results from the differential expression analysis look like.
-dge_df
+
The data we’re interested in here comes from zebrafish samples, so we can obtain just the gene sets relevant to D. rerio with the species argument to msigdbr().
-dr_msigdb_df <- msigdbr(species = "Danio rerio")
+
MSigDB contains 8 different gene set collections (Subramanian et al. 2005).
H: hallmark gene sets
C1: positional gene sets
@@ -3175,21 +4047,20 @@ 4.3 Getting familiar with c
In this example, we will use pathways that are gene sets considered to be “canonical representations of a biological process compiled by domain experts” and are a subset of C2: curated gene sets (Subramanian et al. 2005).
Specifically, we will use the KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways (Kanehisa and Goto 2000).
First, let’s take a look at what information is included in this data frame.
-head(dr_msigdb_df)
+
We will need to use gs_cat and gs_subcat columns to construct a filter step that will only keep curated gene sets and KEGG pathways.
-# Filter the zebrafish data frame to the KEGG pathways that are included in the
-# curated gene sets
-dr_kegg_df <- dr_msigdb_df %>%
- dplyr::filter(
- gs_cat == "C2", # This is to filter only to the C2 curated gene sets
- gs_subcat == "CP:KEGG" # This is because we only want KEGG pathways
- )
-Note: We could have specified that we wanted the KEGG gene sets using the category and subcategory arguments of msigdbr(), but we’re going for general steps! – use ?msigdbr to see more information.
+# Filter the zebrafish data frame to the KEGG pathways that are included in the
+# curated gene sets
+dr_kegg_df <- dr_msigdb_df %>%
+ dplyr::filter(
+ gs_cat == "C2", # This is to filter only to the C2 curated gene sets
+ gs_subcat == "CP:KEGG" # This is because we only want KEGG pathways
+ )
The clusterProfiler() function we will use requires a data frame with two columns, where one column contains the term identifier or name and one column contains gene identifiers that match our gene lists we want to check for enrichment.
Our data frame with KEGG terms contains Entrez IDs and gene symbols.
In our differential expression results data frame, dge_df we have Ensembl gene identifiers. So we will need to convert our Ensembl IDs into either gene symbols or Entrez IDs.
@@ -3200,7 +4071,7 @@ 4.4 Gene identifier conversionThe annotation package org.Dr.eg.db contains information for different identifiers (Carlson 2019). org.Dr.eg.db is specific to Danio rerio – this is what the Dr in the package name is referencing.
Take a look at our other gene identifier conversion examples for examples with different species and gene ID types: the microarray example and the RNA-seq example.
We can see what types of IDs are available to us in an annotation package with keytypes().
-keytypes(org.Dr.eg.db)
+
## [1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS"
## [6] "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME"
## [11] "GO" "GOALL" "IPI" "ONTOLOGY" "ONTOLOGYALL"
@@ -3208,50 +4079,50 @@ 4.4 Gene identifier conversion
Even though we’ll use this package to convert from Ensembl gene IDs (ENSEMBL) to gene symbols (SYMBOL), we could just as easily use it to convert from an Ensembl transcript ID (ENSEMBLTRANS) to Entrez IDs (ENTREZID).
The function we will use to map from Ensembl gene IDs to gene symbols is called mapIds().
-# This returns a named vector which we can convert to a data frame, where
-# the keys (Ensembl IDs) are the names
-symbols_vector <- mapIds(org.Dr.eg.db, # Specify the annotation package
- # The vector of gene identifiers we want to map
- keys = dge_df$Gene,
- # The type of gene identifier we want returned
- column = "SYMBOL",
- # What type of gene identifiers we're starting with
- keytype = "ENSEMBL",
- # In the case of 1:many mappings, return the
- # first one. This is default behavior!
- multiVals = "first"
-)
+# This returns a named vector which we can convert to a data frame, where
+# the keys (Ensembl IDs) are the names
+symbols_vector <- mapIds(org.Dr.eg.db, # Specify the annotation package
+ # The vector of gene identifiers we want to map
+ keys = dge_df$Gene,
+ # The type of gene identifier we want returned
+ column = "SYMBOL",
+ # What type of gene identifiers we're starting with
+ keytype = "ENSEMBL",
+ # In the case of 1:many mappings, return the
+ # first one. This is default behavior!
+ multiVals = "first"
+)
## 'select()' returned 1:many mapping between keys and columns
This message is letting us know that sometimes Ensembl gene identifiers will map to multiple gene symbols. In this case, it’s also possible that a gene symbol will map to multiple Ensembl IDs. For more about how to explore this, take a look at our microarray gene ID conversion example.
Let’s create a two column data frame that shows the gene symbols and their Ensembl IDs side-by-side.
-# We would like a data frame we can join to the differential expression stats
-gene_key_df <- data.frame(
- ensembl_id = names(symbols_vector),
- gene_symbol = symbols_vector,
- stringsAsFactors = FALSE
-) %>%
- # If an Ensembl gene identifier doesn't map to a gene symbol, drop that
- # from the data frame
- dplyr::filter(!is.na(gene_symbol))
+# We would like a data frame we can join to the differential expression stats
+gene_key_df <- data.frame(
+ ensembl_id = names(symbols_vector),
+ gene_symbol = symbols_vector,
+ stringsAsFactors = FALSE
+) %>%
+ # If an Ensembl gene identifier doesn't map to a gene symbol, drop that
+ # from the data frame
+ dplyr::filter(!is.na(gene_symbol))
Let’s see a preview of gene_key_df.
-head(gene_key_df)
+
Now we are ready to add the gene_key_df to our data frame with the differential expression stats, dge_df. Here we’re using a dplyr::left_join() because we only want to retain the genes that have gene symbols and this will filter out anything in our dge_df that does not have gene symbols when we join using the Ensembl gene identifiers.
-dge_annot_df <- gene_key_df %>%
- # Using a left join removes the rows without gene symbols because those rows
- # have already been removed in `gene_symbols_df`
- dplyr::left_join(dge_df,
- # The name of the column that contains the Ensembl gene IDs
- # in the left data frame and right data frame
- by = c("ensembl_id" = "Gene")
- )
+dge_annot_df <- gene_key_df %>%
+ # Using a left join removes the rows without gene symbols because those rows
+ # have already been removed in `gene_symbols_df`
+ dplyr::left_join(dge_df,
+ # The name of the column that contains the Ensembl gene IDs
+ # in the left data frame and right data frame
+ by = c("ensembl_id" = "Gene")
+ )
Let’s take a look at what this data frame looks like.
-# Print out a preview
-head(dge_annot_df)
+
Looks like there are four KEGG sets returned as significant at FDR of 0.1.
@@ -3328,52 +4197,54 @@ 4.8 Run ORA using the enric
4.9 Visualizing results
We can use a dot plot to visualize our significant enrichment results. The enrichplot::dotplot() function will only plot gene sets that are significant according to the multiple testing corrected p values (in the p.adjust column) and the pvalueCutoff you provided in the enricher() step.
-enrich_plot <- enrichplot::dotplot(kegg_ora_results)
-
-# Print out the plot here
-enrich_plot
-
+
+
Use ?enrichplot::dotplot to see the help page for more about how to use this function.
This plot is arguably more useful when we have a large number of significant pathways.
Let’s save it to a PNG.
-ggplot2::ggsave(file.path(plots_dir, "GSE71270_ora_enrich_plot.png"),
- plot = enrich_plot
-)
+
## Saving 7 x 5 in image
We can use an UpSet plot to visualize the overlap between the gene sets that were returned as significant.
-enrichplot::upsetplot(kegg_ora_results)
+
See that KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION and KEGG_LYSOSOME have all their genes in common. Gene sets or pathways aren’t independent!
Let’s also save this to a PNG.
-ggplot2::ggsave(file.path(plots_dir, "GSE71270_ora_upset_plot.png"),
- plot = enrich_plot
-)
+
## Saving 7 x 5 in image
+
+Fill out the pop up window with your email and our Terms and Conditions:
+
+
+
+It may take a few minutes for the dataset to process.
+You will get an email when it is ready.
+
+## About the dataset we are using for this example
+
+For this example analysis, we will use this [medulloblastoma dataset](https://www.refine.bio/experiments/GSE37382) [@Northcott2012].
+The data that we downloaded from refine.bio for this analysis has 285 microarray samples obtained from patients with medulloblastoma.
+Medulloblastoma is the most common childhood brain cancer and is often categorized by subgroups.
+We will use these `subgroup` labels from our metadata to perform differential expression with our GSVA scores.
+
+## Place the dataset in your new `data/` folder
+
+refine.bio will send you a download button in the email when it is ready.
+Follow the prompt to download a zip file that has a name with a series of letters and numbers and ends in `.zip`.
+Double clicking should unzip this for you and create a folder of the same name.
+
+
+
+For more details on the contents of this folder see [these docs on refine.bio](http://docs.refine.bio/en/latest/main_text.html#downloadable-files).
+
+The `GSE37382` folder has the data and metadata TSV files you will need for this example analysis.
+Experiment accession ids usually look something like `GSE1235` or `SRP12345`.
+
+Copy and paste the `GSE37382` folder into your newly created `data/` folder.
+
+## Check out our file structure!
+
+Your new analysis folder should contain:
+
+- The example analysis `.Rmd` you downloaded
+- A folder called "data" which contains:
+ - The `GSE37382` folder which contains:
+ - The gene expression
+ - The metadata TSV
+- A folder for `plots` (currently empty)
+- A folder for `results` (currently empty)
+
+Your example analysis folder should now look something like this (except with respective experiment accession ID and analysis notebook name you are using):
+
+
+
+In order for our example here to run without a hitch, we need these files to be in these locations so we've constructed a test to check before we get started with the analysis.
+These chunks will declare your file paths and double check that your files are in the right place.
+
+First we will declare our file paths to our data and metadata files, which should be in our data directory.
+This is handy to do because if we want to switch the dataset (see next section for more on this) we are using for this analysis, we will only have to change the file path here to get started.
+
+```{r}
+# Define the file path to the data directory
+data_dir <- file.path("data", "GSE37382") # Replace with accession number which will be the name of the folder the files will be in
+
+# Declare the file path to the gene expression matrix file using the data directory saved as `data_dir`
+data_file <- file.path(data_dir, "GSE37382.tsv") # Replace with file path to your dataset
+
+# Declare the file path to the metadata file using the data directory saved as `data_dir`
+metadata_file <- file.path(data_dir, "metadata_GSE37382.tsv") # Replace with file path to your metadata
+```
+
+Now that our file paths are declared, we can use the `file.exists()` function to check that the files are where we specified above.
+
+```{r}
+# Check if the gene expression matrix file is at the file path stored in `data_file`
+file.exists(data_file)
+
+# Check if the metadata file is at the file path stored in `metadata_file`
+file.exists(metadata_file)
+```
+
+If the chunk above printed out `FALSE` to either of those tests, you won't be able to run this analysis _as is_ until those files are in the appropriate place.
+
+If the concept of a "file path" is unfamiliar to you; we recommend taking a look at our [section about file paths](https://alexslemonade.github.io/refinebio-examples/01-getting-started/getting-started.html#an-important-note-about-file-paths-and-Rmds).
+
+# Using a different refine.bio dataset with this analysis?
+
+If you'd like to adapt an example analysis to use a different dataset from [refine.bio](https://www.refine.bio/), we recommend placing the files in the `data/` directory you created and changing the filenames and paths in the notebook to match these files (we've put comments to signify where you would need to change the code).
+We suggest saving plots and results to `plots/` and `results/` directories, respectively, as these are automatically created by the notebook.
+From here you can customize this analysis example to fit your own scientific questions and preferences.
+
+***
+
+ 
+
+It may take a few minutes for the dataset to process.
+You will get an email when it is ready.
+
+## About the dataset we are using for this example
+
+For this example analysis, we will use this [acute viral bronchiolitis dataset](https://www.refine.bio/experiments/SRP140558).
+The data that we downloaded from refine.bio for this analysis has 62 paired peripheral blood mononuclear cell RNA-seq samples obtained from 31 patients.
+Samples were collected at two time points: during their first, acute bronchiolitis visit (abbreviated "AV") and their recovery, their post-convalescence visit (abbreviated "CV").
+
+## Place the dataset in your new `data/` folder
+
+refine.bio will send you a download button in the email when it is ready.
+Follow the prompt to download a zip file that has a name with a series of letters and numbers and ends in `.zip`.
+Double clicking should unzip this for you and create a folder of the same name.
+
+
+
+In this example above, the blue project is ready to be published and there are no other incomplete projects that have been merged in, so all commits from `staging` are wanted to be brought over `master`.
+
+By checking out a separate `publish` branch, we can resolve any merge conflicts in this `publish` branch so we don't do a bad thing by committing conflict resolutions directly to `staging`.
+This also provides us with a "snapshot" if merges are continuing to happen to `staging` on other PRs -- this can make review hard if it keeps changing.
+
+**Scenario 2: Only some changes from staging should be published**
+
+- Create a new branch from the most up-to-date `master` branch, call it `publish-mychanges` where `mychanges` is a short tag relevant to the changes.
+- Checkout your newly created `publish` branch.
+- For each commit that needs to be published, add it to your new branch by using `git cherry-pick 
+
+In this example above, the blue project is ready to be published, but the red project is not.
+So in this scenario, we would cherry pick both commits from the blue project but ignore the one commit from the incomplete red project.
+
+This allows us to move forward changes from `staging` that are ready to be public facing even if other changes aren't ready (or if there isn't great timing for when the other changes will be ready).
+
+_Tips for getting commits ids_
+To perform the cherry picks, you will need the commit ids for finished projects only.
+GitKraken shows commits, but sometimes I find it hard to follow which commits belong to which branch.
+If you checkout `staging` and use `git log` you can see all the most recent commits for the `staging` branch.
+If you want to save it to a file for easy browsing and copy/pasting you can use this command `git log --pretty=format:'%h was %an, %ar, message: %s' > git.log`
+
+#### Publish a change, but quickly: direct merges to master, hotfixes
+
+In (hopefully rare) scenarios where something that has already been published is noted to be broken and should be addressed quickly, [hotfix branch](https://nvie.com/posts/a-successful-git-branching-model/#hotfix-branches) pull requests are allowed.
+These PRs should only be fairly small PRs and not anything that would require intense review.
+This more for situations where "this is broken and here's a fix".
+
+- Create a new branch from the most up-to-date `master` branch and call it `hotfix-mybug`, where where `mybug` is a short tag relevant to the changes.
+- Checkout your newly created `hotfix` branch.
+- Make the hotfix change.
+- Create a [pull request with `master` branch as the base branch](https://docs.github.com/en/free-pro-team@latest/github/collaborating-with-issues-and-pull-requests/creating-a-pull-request#changing-the-branch-range-and-destination-repository).
+- After your PR to `master` is approved and merged, merge the most up-to-date `staging` branch into your `hotfix` branch.
+- File a second PR for `hotfix` but this time with [`staging` as the base branch](https://docs.github.com/en/free-pro-team@latest/github/collaborating-with-issues-and-pull-requests/creating-a-pull-request#changing-the-branch-range-and-destination-repository).
+
+
+
+In this scenario, we start by making a branch from `master` and make the hotfix change on this `hotfix` branch.
+Once the change is approved and merged into `master` from the first PR, we start the second PR and merge in `staging` to our `hotfix` branch.
+This makes sure we don't accidentally undo any changes in `staging` that have not made it to `master` at this time.
+
+#### A summary of types of PRs.
+
+- `some-branch` -> `staging` -> NOT published to user-facing content.
+This will be how most material is prepared so we can be more incremental with our changes and ultimately make sure only the very polished material is made user-facing.
+
+- Group of changes in `staging` -> `master` -> published to user-facing.
+This is for when the updates from the previous kinds of PRs are "ready for prime time".
+
+- Hotfix PR: `hotfix` -> `master` -> published to user-facing.
+For "this-is-broken" type changes that should be hastened to the user-facing content.
+This requires a follow up pull request and merge to `staging`.
+
+### Github actions
+
+This repository has a lot of little moving parts, so to help make sure changes are where they are supposed to be and that some items aren't missed, we use Github actions to automate things where we can.
+
+Github action workflows are initiated either when a pull request is started or when a merge to `master` or `staging` is initiated.
+The following sections explain more details about which Github actions happen when and what they do.
+
+#### Automatic Spell checking and styling
+
+When pull requests are initiated, spell check and styler are run by Github actions.
+
+
+
+If the styling introduces changes to the files, these changes will be committed back to your branch.
+However if there are no styling changes, no commit will be made.
+
+If spell check finds more than 2 errors, Github actions will fail.
+See the [spell check section](#spell-checking) for instructions on how to see your spelling errors and otherwise use spell check.
+
+#### Automatic Docker image and rendering
+
+After a pull request to `staging` or `master` branch is approved and a merge to one of these branches has been initiated, a sequence of Github actions makes sure that the rendered html files are pushed to the correct branch and that the updated docker image is pushed to Dockerhub.
+
+
+
+See the [Docker](#docker-for-refinebio-examples) and the next section about automatic rendering for more on how these steps are conducted.
+
+#### Automatic rendering
This repository uses [snakemake](https://snakemake.readthedocs.io/en/stable/getting_started/installation.html) to render all notebooks.
-All notebooks are automatically re-rendered by GitHub actions upon merges to master.
-The newly rendered html files are all pushed to the `gh-pages` branch which will publish the material to https://alexslemonade.github.io/refinebio-examples/.
+All notebooks are automatically re-rendered by GitHub actions upon merges to `master` or `staging`.
+The newly rendered html files are all pushed to the `gh-pages` branch (if merging to `master`) and `gh-pages-stages` branch (if merging to `staging`).
+
+The `gh-pages` branch is where material is publshed to https://alexslemonade.github.io/refinebio-examples/.
+The `gh-pages-stages` staging branch can be viewed using html preview at `http://htmlpreview.github.io/?https://github.com/AlexsLemonade/refinebio-examples/gh-pages-stages/01-getting-started/getting-started.html`
If this automatic rendering fails, you will see a failed check at at the bottom of your PR on GitHub (and probably an email).
You can see the details of this error on by going to `Actions` and clicking the workflow of PR you are working on that also says `Build Docker` underneath.
@@ -343,6 +498,7 @@ Hopefully the error message helps you track down the problem, but you can also c
### About the render-notebooks.R script
The `render-notebooks.R` script adds a `bibliography:` specification in the `.Rmd` header so all citations are automatically rendered.
+It also adds other components like CSS styling, a footer, and Google Analytics (these items are all hard-coded into the script).
**Options:**
- `--rmd`: provided by snakemake, the input `.Rmd` file to render.
@@ -377,11 +533,3 @@ Follow these steps to add the `.html` link to the navigation bar upon rendering.
6) Replace `tech-section`, `analysis_file_name` with the corresponding file names.
7) Save the file!
8) After you [render the notebook with snakemake](#rendering-notebooks), test the link to make sure it works.
-
-## Pull request status checks
-
-To require that branches are up-to-date with `master` before merging, we need to require that a status check passes before merging to `master`.
-Turning on this setting mitigates the risk that changes that have been merged will be undone by a pull request that was filed first and alters the same file.
-The status check used is a GitHub Action that test builds the docker image.
-Most of the time, this should pull cached docker layers, so this will complete in a matter of minutes.
-As a bonus, this process also checks that any changes to the Dockerfile result in a buildable image.
diff --git a/Snakefile b/Snakefile
index 4909b910..e9c26d36 100644
--- a/Snakefile
+++ b/Snakefile
@@ -8,7 +8,9 @@ rule target:
"02-microarray/dimension-reduction_microarray_01_pca.html",
"02-microarray/dimension-reduction_microarray_02_umap.html",
"02-microarray/gene-id-annotation_microarray_01_ensembl.html",
- "02-microarray/pathway-analysis_microarray_02_ora.html",
+ "02-microarray/pathway-analysis_microarray_01_ora.html",
+ "02-microarray/pathway-analysis_microarray_02_gsea.html",
+ "02-microarray/pathway-analysis_microarray_03_gsva.html",
"02-microarray/ortholog-mapping_microarray_01_ensembl.html",
"03-rnaseq/00-intro-to-rnaseq.html",
"03-rnaseq/clustering_rnaseq_01_heatmap.html",
@@ -17,7 +19,9 @@ rule target:
"03-rnaseq/dimension-reduction_rnaseq_02_umap.html",
"03-rnaseq/gene-id-annotation_rnaseq_01_ensembl.html",
"03-rnaseq/ortholog-mapping_rnaseq_01_ensembl.html",
- "04-advanced-topics/00-intro-to-advanced-topics.html"
+ "04-advanced-topics/00-intro-to-advanced-topics.html",
+ "04-advanced-topics/network-analysis_rnaseq_01_wgcna.html"
+
rule render_citations:
input:
diff --git a/components/_navbar.html b/components/_navbar.html
index 1ed6c957..fe4358c1 100644
--- a/components/_navbar.html
+++ b/components/_navbar.html
@@ -7,7 +7,7 @@
-
+