| Table of Contents |
|---|
...
Go to the following page https://posit.co/download/rstudio-desktop/ and follow the instructions provided to install first R and then Rstudio.
...
The basis for the differential expression analysis is a count table of sequence reads mapped to defined gene regions per sample. There are a variety of methods to generate this count table, but for this exercise, we will be using the output from the Nextflow nfcore/rnaseq analysis you completed in the previous workshop sessions.
To access this count table:
Go to the W:\training\rnaseq\runs\run3_RNAseq\results folder that contains the results from running the nfcore/rnaseq pipeline. The output folders from task 3 look like this:
...
The count table can be found in the star_salmon folder. A list of files and folders in the star_salmon folder will look like this:
...
Preparing your data. Two data files are needed for this analysis: a samples table and your count table
Install required R packages (only need to run once) - after installation, we only need to load the packages. NOTE: If using an rVDI virtual machine, the R packages are already installed
Load required R packages. Unlike installing the packages, this needs to be done every time you run the analysis
Import your data files (count table and samples table)into R
Checking for outliers and batch effects
PCA plot
Pairwise samples heatmap
Identify differentially expressed (DE) genes using DESeq2
Annotating your DE genes
Volcano plot
DE genes heatmap
...
a. First create a new folder in H:\workshop\RNAseq . Call it something suitable, such as ‘DE_analysis_workshop’
b. Create a sub folder subfolder here called ‘data’. This is where your two data files will be stored
...
| Code Block |
|---|
gene_id gene_name SRR20622172 SRR20622173 SRR20622174 SRR20622175 SRR20622176 SRR20622177 SRR20622178 SRR20622179 SRR20622180 ENSMUSG00000000001 Gnai3 7086 4470 2457.002 2389 6398 2744 2681 3961 4399 ENSMUSG00000000003 Pbsn 0 0 0 0 0 0 0 0 0 ENSMUSG00000000028 Cdc45 1232.999 827 42 57 1036 55 78 88 89 ENSMUSG00000000031 H19 200 139 2 0 143.622 1 17.082 24 16.077 ENSMUSG00000000037 Scml2 70 57.001 8 8 66.999 16 23 27.999 29 ENSMUSG00000000049 Apoh 0 0 1 0 2 2 1 3 0 ENSMUSG00000000056 Narf 1933 1480 519 497 1730 539 365 458 536 ENSMUSG00000000058 Cav2 6008 3417 1347.001 1344 5482 1367 2669.001 4358 4365.832 ENSMUSG00000000078 Klf6 3809 2732 4413.001 3483.978 3559 4491 3209 3980 4626 |
d. In the same W:\training\rnaseq\runs\run3_RNAseq\results\star_salmon directory there will be a file called metadata.xlsx . Copy this file to your ‘data’ folder as well. This file will normally need to be manually created by you to match your sample IDs and treatment groups, but we created this file already for you to use. This samples table needs 3 columns called ‘sample_name’, containing the sample names seen in the count table (column names), ‘sample_ID’, which is the (less messy) names you want to call the samples in this analysis workflow, and ‘group’, which contains the treatment groups each sample belongs to. The contents of this file look like this:
...
4. Import your data files into R
In this section, we will import your count table and samples table into R.
You’ll need to change the ‘setwd' line to your working directory. Click ‘Session’ → ‘Set working directory’ → ‘Choose working directory’ and then choose the analysis workshop directory you created previously , that contains your R script file and the 'Data’ directory.
| Code Block |
|---|
#### 4. Import your count data ####
# Make sure you have: a) your count table (salmon.merged.gene_counts.tsv file, if you used Nextflow nfcore/rnaseq to analyse your data). Copy this to a subdirectory called 'data'. b) your metadata file. This should be either an Excel file called 'metadata.xlsx' or a tab-separated text file called 'metadata.txt'. It needs 3 columns called 'sample_name', 'sample_ID' and 'group'. The sample names should be EXACTLY the same as the names in the count table. These names are often uninformative and long, so the 'sample_ID' is the sample labels you want to put on your plots. E.g. if you have a 'high fat' group, you might want to rename the samples HF1, HF2, HF3, etc)
## USER INPUT
# Set working directory.
# Change this to your working directory (In the RStudio menu: Session -> Set working directory -> Choose working directory)
setwd("CH:/Users/whatmorp/OneDrive - Queensland University of Technology/Desktop/Projects/RNA-Seq downstream analysisworkshop/RNAseq")
# Import your count data. make sure you've created a 'data' subdirectory and put the count table file there.
metacountdata <- read.table("./data/salmon.merged.gene_counts.tsv", header = TRUE, row.names = 1)
# Import metadata. Again, need a metadata.xlsx file in the data subdirectory.
meta <- read_excel("./data/metadata.xlsx")
# Remove 1st columns of metadata (gene_name)
counts <- metacountdata[ ,2:ncol(metacountdata)]
# Rename sample names to new sample IDs
counts <- counts[as.character(meta$sample_name)]
colnames(counts) <- meta$sample_ID
# Counts need to be rounded to integers
counts <- ceiling(counts) |
...
5. Checking for outliers and batch effects
In this section, we will create PCA plots and heatmaps to examine the relationships between samples. Outlier samples and batch effects can heavily bias your results and should be addressed (e.g. removal of outlier samples from the dataset) before any differential expression analysis is completed.
First, we need to prepare the data for plotting. Copy, paste, and run the following code in your R script. you will need to input which treatment groups you wish to plot (plotgroups <-) from the set of available treatment groups (which you can find out with unique(meta$group).
...