Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

RNAseq on Cavatica tutorial #316

Merged
merged 66 commits into from
Feb 24, 2021
Merged
Show file tree
Hide file tree
Changes from 65 commits
Commits
Show all changes
66 commits
Select commit Hold shift + click to select a range
7edd415
New RNAseq on cavatica tutorial
s-canchi Jan 29, 2021
04cac22
fixed right navigation render
s-canchi Jan 29, 2021
fc45643
fixed broken link
s-canchi Jan 29, 2021
c551522
minor table formatting
s-canchi Jan 29, 2021
4e521a6
removed instant nav
s-canchi Jan 29, 2021
a98b780
added pages and nav instant
s-canchi Jan 29, 2021
0b8bd0f
Merge branch 'dev' into scanchi-rnaseq-new
s-canchi Jan 29, 2021
cc0e08f
marisa's fixes to the overview page
s-canchi Feb 8, 2021
124bc26
edits to intro to rnaseq lesson
s-canchi Feb 8, 2021
dd915aa
fixed sentences, reorder, reformat and edits
s-canchi Feb 8, 2021
18e8af2
punctuation and grammar checks
jeremywalter Feb 9, 2021
7beb638
fixes to cavatica explore lesson, new screenshots, text edits, reformat
s-canchi Feb 9, 2021
84d0b7c
Merge branch 'scanchi-rnaseq-new' of https://github.com/nih-cfde/trai…
s-canchi Feb 9, 2021
19c1232
fixed time scale on overview page
s-canchi Feb 9, 2021
1a5de23
fixes to deseq setup lesson, new info on app versions, text edits and…
s-canchi Feb 9, 2021
4defd50
minor grammar text edits
jeremywalter Feb 9, 2021
05b40f6
minor grammar fixes
jeremywalter Feb 9, 2021
614f57d
highlighting file name
s-canchi Feb 9, 2021
72d4d94
formatting symbols
s-canchi Feb 9, 2021
dc7e66c
updated with new screeshots, text edits, reformat and restructure in …
s-canchi Feb 9, 2021
37ebe02
adding line break
s-canchi Feb 9, 2021
e2c87a5
reformat line breaks'
s-canchi Feb 9, 2021
5771608
resize file image
s-canchi Feb 9, 2021
231eefa
redid excel images
s-canchi Feb 9, 2021
8431daa
added vidlet tip, reformat text in app lesson, minor txt corrections
s-canchi Feb 9, 2021
8fc90f3
Merge branch 'dev' into scanchi-rnaseq-new
s-canchi Feb 9, 2021
ac6d5be
updated data cruncher lesson, reformat, txt edits and corrections
s-canchi Feb 12, 2021
448ef56
Merge branch 'dev' into scanchi-rnaseq-new
s-canchi Feb 15, 2021
68a2310
fix landing pg table links
marisalim Feb 18, 2021
1ce72ba
minor text, format edits. clarified some input steps, streamlined r s…
marisalim Feb 18, 2021
2035f6f
added commas and DGE to 1.md
jeremywalter Feb 19, 2021
fdb7be5
RStudio word tweak 1.md
jeremywalter Feb 19, 2021
30566ef
changes to tumor descriptions (grammatical) 2.md
jeremywalter Feb 19, 2021
3d81cf5
missed a full-stop 2.md
jeremywalter Feb 19, 2021
b09fceb
add KF portal link 3.md
jeremywalter Feb 19, 2021
01d50a0
update rna_seq4.md - add more detailed selection directions
jeremywalter Feb 19, 2021
a100658
added some granularity to filtering steps 4.md
jeremywalter Feb 19, 2021
9e74c69
add tag name to instructions
jeremywalter Feb 19, 2021
90849fc
grammar checks md.5
jeremywalter Feb 19, 2021
edf5537
text changes and clarification
jeremywalter Feb 19, 2021
231f34e
more details for tumor filtering - add refresh mention
jeremywalter Feb 19, 2021
d471de5
Merge branch 'dev' into scanchi-rnaseq-new
s-canchi Feb 21, 2021
fa78c34
Added icons and reformat txt in lesson 4
s-canchi Feb 22, 2021
f66a4fc
added icons and minor text formatting in lesson 5
s-canchi Feb 22, 2021
0f01c83
added icons to lesson 6
s-canchi Feb 22, 2021
bacd9a4
added icons and minor txt edits in lesson 7
s-canchi Feb 22, 2021
4933bcc
added icons and minor text edits in lesson 8
s-canchi Feb 22, 2021
f8005cf
Merge branch 'dev' into scanchi-rnaseq-new
s-canchi Feb 22, 2021
e3093a3
added KF save link and note for cavatica login
s-canchi Feb 22, 2021
17b4c7f
edit link to kF save lesson
s-canchi Feb 22, 2021
107c127
updated screenshots t match new query outcome and edit query link
s-canchi Feb 23, 2021
b100989
minor txt correction
s-canchi Feb 23, 2021
ff306d0
clarify excel use
s-canchi Feb 23, 2021
43221b7
Merge branch 'dev' into scanchi-rnaseq-new
s-canchi Feb 23, 2021
4877110
text formatting and add histology_type selection
jeremywalter Feb 23, 2021
814c0f7
remove "" for tumor entry instructions
jeremywalter Feb 23, 2021
15428cc
grammar and instructions fixes, added a time marker text
jeremywalter Feb 23, 2021
79e58ce
added some clarification for Steps 3, 3a, 3b
jeremywalter Feb 23, 2021
3aece26
removed extensions to test render
s-canchi Feb 23, 2021
9f51932
removed unused extension
s-canchi Feb 23, 2021
5e9ecc8
8.md last fixes (fixing the step 3 fixes)
jeremywalter Feb 24, 2021
a9bead5
add an s
jeremywalter Feb 24, 2021
5d47c7c
typos on RNAseq
marisalim Feb 24, 2021
4bfd5f5
Merge branch 'dev' into scanchi-rnaseq-new
marisalim Feb 24, 2021
4534821
clarify r script 3b instructions
marisalim Feb 24, 2021
cf1783a
final word tweaks, added ext back, added missing icons and end of les…
s-canchi Feb 24, 2021
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions docs/Bioinformatics-Skills/.pages
Original file line number Diff line number Diff line change
Expand Up @@ -9,3 +9,4 @@ nav:
- GWAS in the Cloud: GWAS-in-the-cloud
- Snakemake Workflow Management: Snakemake
- Simulate_Illumina_Reads.md
- RNAseq-on-Cavatica
11 changes: 11 additions & 0 deletions docs/Bioinformatics-Skills/RNAseq-on-Cavatica/.pages
Original file line number Diff line number Diff line change
@@ -0,0 +1,11 @@
nav:
- rna_seq_1.md
- rna_seq_2.md
- rna_seq_3.md
- rna_seq_4.md
- rna_seq_5.md
- rna_seq_6.md
- rna_seq_7.md
- rna_seq_8.md


Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Original file line number Diff line number Diff line change
@@ -0,0 +1,172 @@
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Author: Saranya Canchi
# Filename: Cancer_DGE_Analysis.R
# Purpose: R script for differential gene expression between
# pediatric cancer types using DESeq2
# Version: 1.0
# Date: 01/22/2021
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Install libraries####

if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install(c("tximport",
"regionReport",
"org.Hs.eg.db",
"pcaExplorer",
"BiocStyle"),
update=FALSE,
ask=FALSE)

# Load libraries####

library(GenomicFeatures)
library(DESeq2)
library(tximport)
library(org.Hs.eg.db)
library(pcaExplorer)
library(regionReport)
library(ggplot2)
library(knitr)

# Generate transcript and gene name table for gene level summary####
## Using GenomicFeatures pkg to read in the reference GTF files - following steps from DESeq2 app

data_dir <- "/sbgenomics/project-files"
txdb <- makeTxDbFromGFF(file= file.path(data_dir,"Homo_sapiens.GRCh38.84.gtf"))
k <- keys(txdb, keytype = "TXNAME")

## HGNC gene name not available in list of filters. Using Ensemble gene name for one to one mapping.

tx2gene <- select(txdb, k, "GENEID", "TXNAME")
head(tx2gene)

# Read the phenotype data####

pheno_data <- read.csv(file.path(data_dir,"phenotype_filtered.csv"))
str(pheno_data)

## Converting covariates of interest to factors
### Setting Ependymoma as reference factor level for histology variable

pheno_data$histology <- factor(pheno_data$histology,
levels=c("Ependymoma", "Medulloblastoma"))
pheno_data$tumor_location <- factor(pheno_data$tumor_location)
pheno_data$diagnosis_age_range <- factor(pheno_data$diagnosis_age_range)

# Generate gene level summary using tximport####

head(list.files(data_dir))
files <- file.path(data_dir, pheno_data$name)
names(files) <- pheno_data$sample_id
head(files)
txi_sum <- tximport(files,
type="kallisto",
tx2gene=tx2gene,
ignoreTxVersion = TRUE)
names(txi_sum)
head(txi_sum$counts)

# DESeq2 import and analysis####

dds_cancer <- DESeqDataSetFromTximport(txi_sum,
colData = pheno_data,
design = ~ diagnosis_age_range + tumor_location + histology)
head(rownames(dds_cancer))

## Generating additional gene info to add to dds object

cancer_gene_map <- mapIds(org.Hs.eg.db,
keys=rownames(dds_cancer),
column="SYMBOL",
keytype="ENSEMBL",
multiVals="first")
cancer_gene_map <- stack(cancer_gene_map)
head(cancer_gene_map)
colnames(cancer_gene_map)=c("gene_symbol","ensembl_gene_id")

### Check row for row line up and add gene symbol column

all(rownames(dds_cancer) == cancer_gene_map$ensembl_gene_id)
mcols(dds_cancer) <- cbind(mcols(dds_cancer), cancer_gene_map$gene_symbol)
colnames(rowData(dds_cancer))[1]="gene_symbol"

# PCA interactive analysis - pcaExplorer####
## Select Try again in the popup window to open in a new browser tab.

pcaExplorer(dds=dds_cancer,
annotation=cancer_gene_map)

## Alternatively use Variance stabilization transformation to generate PCA plot

output_dir <- "/sbgenomics/output-files"
vsd_cancer <- vst(dds_cancer, blind=FALSE)

### Using histology and tumor_location variables

png(filename=file.path(output_dir,"PCA_histology_tumor-location.png"))
plotPCA(vsd_cancer, intgroup=c("histology","tumor_location"))
dev.off()

### Using histology and diagnosis_age_range variables

png(filename=file.path(output_dir,"PCA_histology_diagnosis-age.png"))
plotPCA(vsd_cancer, intgroup=c("histology","diagnosis_age_range"))
dev.off()

### Note: The PCA plots in the report generated by the DESeq2 app use rlog transformation.
# We can create the same plot using the following code:
#rld_cancer <- rlog(dds_cancer, blind=FALSE)
#plotPCA(rld_cancer, intgroup=c("histology"))


# DGE analysis####

dds_cancer <- DESeq(dds_cancer,
fit="parametric",
test='Wald',
betaPrior=TRUE,
parallel=TRUE)

res_cancer <- results(dds_cancer,
contrast=c("histology","Medulloblastoma","Ependymoma"),
alpha=0.05)
summary(res_cancer)
resultsNames(dds_cancer)

plotMA(res_cancer,
ylim=range(res_cancer$log2FoldChange, na.rm=TRUE))

## Write the results to a table

res_cancer_order <- res_cancer[order(res_cancer$pvalue),]
write.csv(as.data.frame(res_cancer_order),
file=file.path(output_dir,"Cancer_DESeq2_DGE_results.csv"))

## Write normalized counts to a file

norm_counts <- counts(dds_cancer,normalized=TRUE)
write.table(norm_counts,
file=file.path(output_dir,"Cancer_DESeq2_normalized_counts.txt"),
sep="\t",
quote=FALSE,
col.names = NA)

# Create a report with all the visualization from the DESeq2 vignette####

dir.create(file.path(output_dir,"DESeq2-Report"),
showWarnings = FALSE,
recursive = TRUE)

report <- DESeq2Report(dds = dds_cancer,
project = 'Pediatric Cancer DGE Analysis with DESeq2',
intgroup = c('histology','diagnosis_age_range'),
res = res_cancer,
outdir = file.path(output_dir,"DESeq2-Report"),
output = 'Cancer-DESeq2-Report',
theme = theme_bw())

save.image(file=file.path(output_dir,paste0("Cancer_DGE_",Sys.Date(), ".env.Rdata")))


Original file line number Diff line number Diff line change
@@ -0,0 +1,110 @@
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Author: Saranya Canchi
# Filename: Cancer_DGE_Analysis_Automate.R
# Purpose: R script for differential gene expression between
# pediatric cancer types using DESeq2
# Version: 1.0
# Date: 01/22/2021
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Install libraries####

if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install(c("tximport",
"regionReport",
"BiocStyle"),
update=FALSE,
ask=FALSE)

# Load libraries####

library(GenomicFeatures)
library(DESeq2)
library(tximport)
library(regionReport)
library(ggplot2)
library(knitr)

# Generate transcript and gene name table for gene level summary####
## Using GenomicFeatures pkg to read in the reference GTF files - following steps from DESeq2 app

data_dir <- "/sbgenomics/project-files"
txdb <- makeTxDbFromGFF(file= file.path(data_dir,"Homo_sapiens.GRCh38.84.gtf"))
k <- keys(txdb, keytype = "TXNAME")

## HGNC gene name not available in list of filters. Using Ensemble gene name for one to one mapping.

tx2gene <- select(txdb, k, "GENEID", "TXNAME")

# Read the phenotype data####

pheno_data <- read.csv(file.path(data_dir,"phenotype_filtered.csv"))

## Converting covariates of interest to factors
### Setting Ependymoma as reference factor level for histology variable

pheno_data$histology <- factor(pheno_data$histology, levels=c("Ependymoma", "Medulloblastoma"))
pheno_data$tumor_location <- factor(pheno_data$tumor_location)
pheno_data$diagnosis_age_range <- factor(pheno_data$diagnosis_age_range)

# Generate gene level summary using tximport####

head(list.files(data_dir))
files <- file.path(data_dir, pheno_data$name)
names(files) <- pheno_data$sample_id
txi_sum <- tximport(files,
type="kallisto",
tx2gene=tx2gene,
ignoreTxVersion = TRUE)

# DESeq2 import and analysis####

dds_cancer <- DESeqDataSetFromTximport(txi_sum,
colData = pheno_data,
design = ~ diagnosis_age_range + tumor_location + histology)

# DGE analysis####

dds_cancer <- DESeq(dds_cancer,
fit="parametric",
test='Wald',
betaPrior=TRUE,
parallel=TRUE)

res_cancer <- results(dds_cancer,
contrast=c("histology","Medulloblastoma","Ependymoma"),
alpha=0.05)

## Write the results to a table

output_dir <- "/sbgenomics/output-files"
res_cancer_order <- res_cancer[order(res_cancer$pvalue),]
write.csv(as.data.frame(res_cancer_order),
file=file.path(output_dir,"Cancer_DESeq2_DGE_results.csv"))

## Write normalized counts to a file

norm_counts <- counts(dds_cancer,normalized=TRUE)
write.table(norm_counts,
file=file.path(output_dir,"Cancer_DESeq2_normalized_counts.txt"),
sep="\t",
quote=FALSE,
col.names = NA)

# Create a report with all the visualization from the DESeq2 vignette####

dir.create(file.path(output_dir,"DESeq2-Report"),
showWarnings = FALSE,
recursive = TRUE)

report <- DESeq2Report(dds = dds_cancer,
project = 'Pediatric Cancer DGE Analysis with DESeq2',
intgroup = c('histology','diagnosis_age_range'),
res = res_cancer,
outdir = file.path(output_dir,"DESeq2-Report"),
output = 'Cancer-DESeq2-Report',
theme = theme_bw())

save.image(file=file.path(output_dir,paste0("Cancer_DGE_",Sys.Date(), ".env.Rdata")))

60 changes: 60 additions & 0 deletions docs/Bioinformatics-Skills/RNAseq-on-Cavatica/rna_seq_1.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,60 @@
---
layout: page
title: RNAseq Tutorial Overview
---

Differential Gene Expression Analysis on Cavatica Cloud Platform
====================================================================

**RNA sequencing (RNAseq)** is a high throughput technique that provides qualitative and quantitative information about RNA biology including transcriptome-wide expression quantification, discovery of novel genes and gene isoforms, and differential expression.

The goal of this tutorial is to enable you to: </br>

1. create virtual cancer cohorts using the NIH Common Fund-supported Gabriella Miller Kids First Data portal (KF Portal). </br>
2. analyze the differential gene expression (DGE) on Cavatica, an integrated cloud based platform.

You will learn two different approaches for DGE analysis using open access human cancer data on Cavatica: **(a)** using a public workflow app and **(b)** running code from an analysis script on an instance with RStudio computational environment.

**Table of contents**

| Est. Time| Lesson Name | Description|
| ---|--------|--------|
| 10 mins |[An Introduction to RNAseq](./rna_seq_2.md)| Background about RNAseq
| 20 mins |[Selecting Kids First Cancer Cohort](./rna_seq_3.md)| Select Kids First open access cancer RNAseq files and push to Cavatica |
| 20 mins |[Cavatica - View, Filter, Tag and Download](./rna_seq_4.md) | Filter imported data, tag and download relevant metadata from Cavatica |
| 20 mins |[Setup DESeq2 Public App](./rna_seq_5.md)| Setting up the workflow app based on DESeq2 on Cavatica |
| 15 mins |[Phenotype File and Upload to Cavatica](./rna_seq_6.md) | Reformat metadata file and upload it to Cavatica |
| 50 mins |[Analysis with DESeq2 Public App](./rna_seq_7.md) | Run the DESeq2 app with appropriate inputs and computational settings |
| 60 mins |[Analysis using Data Cruncher](./rna_seq_8.md) | Analysis on an instance in the RStudio environment |

!!! note "Learning Objectives"
* learn to build virtual cohorts on KF portal
* learn to navigate project folder and perform file operations on Cavatica
* learn to upload and download data from Cavatica
* learn to search, copy, and edit public workflow apps on Cavatica
* learn to perform differential gene expression (DGE) analysis using DESeq2 app
* learn to setup analysis environment and execute code for DGE analysis

=== "Prerequisites"
* Setup: Integrated login accounts on Kid's First Data Portal & Cavatica - Follow our lessons on account setup and connecting the two accounts.
- [Setup Kids First account](../Kids-First/Portal-Setup-And-Permissions/KF_3_KF_Registration.md)
- [Register for Cavatica](../Kids-First/Portal-Setup-And-Permissions/KF_4_Cavatica_Registration.md)
- [Connect KF and Cavatica](../Kids-First/Portal-Setup-And-Permissions/KF_5_ConnectingAccounts.md)

!!! note "Login Credentials"

You do not need **eRA Commons ID** to do the lesson!

* Background: Knowledge of biology and rudimentary genetics.
* Technology: Basic knowledge of R and command line. Familiarity with RStudio is useful.
* Financial: Pilot funds ($100) are provided to every user on Cavatica with linked KF accounts.
* Time: Initial account setup may take hours to a day for verification. Setup of eRA Commons ID may take days and is institute dependent.
=== "Est.Cost"
- DESeq2 app < $1.00
- Analysis with R < $1.00
=== "Tutorial Resources"
- [Kids First Data Portal](https://kidsfirstdrc.org)
- [Cavatica Documentation](https://docs.cavatica.org/docs/getting-started)
- [Playlist of video tutorials explaining concepts used in RNAseq analysis](https://www.youtube.com/playlist?list=PLblh5JKOoLUJo2Q6xK4tZElbIvAACEykp)
- [DESeq2 vignette](https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#how-do-i-use-vst-or-rlog-data-for-differential-testing)
- [tximport](https://bioconductor.org/packages/release/bioc/vignettes/tximport/inst/doc/tximport.html)
Loading