The rna-seq-tutorial-for-gene-differential-expression-analysis from amarinderthind

RNASeq tutorial for gene differential expression analysis and Functional enrichment analysis

This tutorial was created for educational purposes and was presented at a Workshop organised by Dollar Education.

Interested in exploring more applications of the RNASeq, Read here https://ro.uow.edu.au/test2021/3578/

Want to adjust for tumour purity Check https://www.nature.com/articles/s41587-022-01440-w

About the RNA-Seq analysis

The R script performs several steps in RNAseq gene differential expression analysis, including filtering, preprocessing, visualization, clustering, and Enrichment. For the analysis, several R Bioconductor packages are required to be installed (Installation commands are provided in the script. However, users can also refer to the Bioconductor website for detailed instructions).

Main Bioconductor packages to be installed

DESeq2, edgeR, biomaRt (Very useful for gene filtering and annotations), PCAtools (PCA detailed analysis), ReactomePA (enrichment analysis)

RNA-Seq-DGE.R is the R script.

RNA-Seq-DGE.rmd used to create an output of the script shown in the PDF file here.

No significant enrichment was found from the demo example, so enrichment plots are empty or commented.

Note:

If data is obtained by different batch processing please consider ~batch (batch effect in the design matrix) or use CombatSeq as defined in the new script.

Required data files

You should have a raw count and annotation/metadata file for running this analysis. Raw count files are usually obtained from tools such as featureCount, Rsem etc.

setwd("/Users/Path") #Path_to_working_directory

rawcount<-read.table ("RawCount_input.csv",header=TRUE,  sep=",",  row.names=1)

## Replace NAs by zero and change the input to the required format
rawcount <- round(rawcount) 
rawcount[is.na(rawcount)] <- 0

Loading and filtering Data annotation file

anno <-read.table ("Annotation_of_samples_12_Samples_ALL.csv",header=TRUE,  sep=",", row.names = 1) ##, In this case, we have 3 columns (a) sample (b) Condition (c) batch
#rownames(anno) <- anno$sample  ##add row names as sample name (if not already), because PCA function check row names of anno == col of the data matrix

table(anno$Condition)

library(tidyverse)
library('dplyr') ##HAS COUNT FUNCTION

### In case you want to consider the subset of samples based on some condition (when multiple e.g. >3 )
#anno <- anno %>% 
#  as.data.frame %>%
#  filter(anno$Condition =='Condition_A' |anno$Condition =='Condition_B' | anno$Condition == 'Condition_C' )  %>%
#  arrange(Condition)  	#Arrange rows by padj values

sort anno based on condition ## sorted conditional sample representation in heatmap

anno <- anno %>% 
  as.data.frame %>%
  arrange(Condition)

PCA plot for pre-DE investigation

library(PCAtools)

anno <- anno[match(colnames(rawcount), anno$Sample),] ## reordering anno rows with colnames of rawcount
lograwcount <- as.matrix(log2(rawcount +1))  ## log transformation of rawcount for PCA plot 

 top1000.order <- head(order(matrixStats::rowVars(lograwcount), decreasing = TRUE), 1000) ## taking top 1000 genes having the highest variance selected from all the genes in the input
  p <- PCAtools::pca(mat = lograwcount[top1000.order,], metadata = anno, removeVar = 0.01) ## performing PCA

  biplot(p,                                                       #visualization of PCA plot
       lab = paste0(p$metadata$Sample),
        colby = 'Batch',  #Sample #Batch #Condition #sex
        hline = 0, vline = 0,
        legendPosition = 'right',
         encircle = T )
  
  screeplot(p, axisLabSize = 18, titleLabSize = 22) ## This plot shows how much variation in the data is explained by which PC component.
  pairsplot(p) ## draw various combinations of the PCA plot