...

Monitor the progress:

qjobs

DESeq2

3. Loading required packages

This section needs to be run every time

Load packages

bioconductor_packages <- c("DESeq2", "EnhancedVolcano")
cran_packages <- c("ggrepel", "ggplot2", "plyr", "reshape2", "FactoMineR", "factoextra", "pheatmap")
lapply(cran_packages, require, character.only = TRUE)
lapply(bioconductor_packages, require, character.only = TRUE)

4. Import your count data

Set working directory.

Change this to your working directory

Your home wd

setwd("C:/Users/whatmorp/OneDrive - Queensland University of Technology/Desktop/Teaching/HPC_training_workshops/smRNA_seq")

HPC wd

#setwd("/work/training/smallRNAseq/")

Import your count data. make sure you've created a 'data' subdirectory and put the count table file there.

metacounts <- read.table("./data/mature_counts.txt", header = TRUE, row.names = 1)

Import metadata. Again, need a metadata_microRNA.txt file in the data subdirectory.

meta <- read.table("./data/metadata_microRNA.txt", header = TRUE)

Rename sample names to new sample IDs

counts <- metacounts[as.character(meta$sample_name)]
colnames(counts) <- meta$sample_ID

5. Outliers and batch effects

This section normalises and transforms the count data so that it can be plotted on a PCA plot and a heatmap

USER INPUT

Choose the groups you want to plot in a PCA/Heatmap. You can select any 2 or more of the groups (or all of the groups) you have in your 'groups' column of your metadata table

To see what groups are present, run the following:

unique(meta$group)

Now add which groups you want to plot (i.e. replace the groupnames below, and add more, separated by a comma and in "quotes", as needed). NOTE: R is case-sensitive, so these group names must be named EXACTLY the same as in the metadata table.

plotgroups <- c("normal", "Huntingtons_disease")

Pull out only the counts from the above groups

groupcounts <- counts[meta$group %in% plotgroups]

Normalise counts by library size, using DeSeq2's estimateSizeFactors() function. Note that DeSeq2 does this internally during DEG calling. The normalisation below is done separately for PCA and density plotting.

Set up the initial DeSeq2 experimental parameters.

condition <- factor(1:length(groupcounts))

Set up the column data. A data frame of sample ID's and conditions

coldata <- data.frame(row.names=colnames(groupcounts), condition)

Set up the DeSeq2 data set structure

f <- DESeqDataSetFromMatrix(countData = groupcounts, colData = coldata, design= ~ condition)

Estimate the size factors. See DeSeq2 manual for details

f <- estimateSizeFactors(f)

Size factors can be viewed by: sizeFactors(f)

Multiply each row (sample) by the corresponding size factor

subcount_norm <- as.matrix(groupcounts) %*% diag(sizeFactors(f))

Re-add column names

colnames(subcount_norm) <- colnames(groupcounts)

Remove low coverage transcripts (mean count < 10)

Find the mean of each row (and output as a data frame object)

means <- as.data.frame(rowMeans(subcount_norm))

Then join the means data with the counts

means <- cbind(means, subcount_norm)

Then subset out only genes with mean > 10

data <- subset(means, means[ , 1] > 10)

Remove the means column

data <- data[,-1]

Transform data

data_log <- vst(round(as.matrix(data)), nsub = nrow(data)-20)

Transformation can create some infinite values. Can't generate PCA data on these. Can see how many by: sum(sapply(data_log, is.infinite))

To remove infinite rows, use 'is.finte' or '!is.infinite'

data_log <- data_log[is.finite(rowSums(data_log)),]
colnames(data_log) <- colnames(groupcounts)

Set up the PCA plot base data

We're using the FactoMineR package to generate PCA plots (http://factominer.free.fr/index.html)

Need to transpose the data first

data_log_t <- t(data_log)

Add the group data

data_log_t_vars <- data.frame(meta$group[meta$group %in% plotgroups], data_log_t)

Generate the PCA data using FactoMineR package

res.pca <- PCA(data_log_t_vars, quali.sup = 1, graph=FALSE)

Set up the dendogram/heatmaps base data

Calculate the distance matrix:

distance_matrix <- as.matrix(dist(t(data_log)))

5a. PCA plot

Generate the PCA plot. Groups are shaded with ellipses at 95% confidence level. NOTE: at least 4 replicates need to be in a group for an ellipses to be drawn.

NOTE: change the group point colours by changing 'palette = ' below. Use the 'RColourBrewer' colour names (https://r-graph-gallery.com/38-rcolorbrewers-palettes.html). For example, if you are plotting 3 groups and choose palette = "Set1", this will use the first 3 colours from the Set1 colour palette.

p <- fviz_pca_ind(res.pca,
geom.ind = c("point", "text"), # show points only (but not "text")
col.ind = meta$group[meta$group %in% plotgroups], # color by groups
pointsize = 5, label = "all", title = "", legend.title = "Treatment groups", palette = "Dark2",
addEllipses = TRUE, ellipse.type = "t", ellipse.level = 0.95) + theme(legend.text = element_text(size = 12), legend.title = element_text(size = 14), axis.title=element_text(size=16), axis.text=element_text(size=14))

p

Output as publication quality (300dpi) tiff and pdf.

This will name your output files with the treatment groups you selected.

Create a 'results_outliers_included' subdirectory where all results_outliers_included will be output

dir.create("results_outliers_included", showWarnings = FALSE)

Create a (300dpi) tiff

ggsave(file = paste0("./results_outliers_included/PCA_", paste(plotgroups, collapse = "Vs"), ".tiff"), dpi = 300, compression = "lzw", device = "tiff", width = 10, height = 8, plot = p)

Create a pdf

ggsave(file = paste0("./results_outliers_included/PCA_", paste(plotgroups, collapse = "Vs"), ".pdf"), device = "pdf", width = 10, height = 8, plot = p)

5b. Samples heatmap and dendrogram

This section plots a heatmap and dendrogram of pairwise relationships between samples. In this way you can see if samples cluster by treatment group.

See here: https://davetang.org/muse/2018/05/15/making-a-heatmap-in-r-with-the-pheatmap-package/

Define annotation column

annot_columns <- data.frame(meta$group[meta$group %in% plotgroups])

Make the row names the sample IDs

row.names(annot_columns) <- meta$sample_ID[meta$group %in% plotgroups]
colnames(annot_columns) <- "Treatment groups"

Need to factorise it

annot_columns[[1]] <- factor(annot_columns[[1]])

Generate dendrogram and heatmap

pheatmap(distance_matrix, color=colorRampPalette(c("white", "#9999FF", "#990000"))(50), cluster_rows = TRUE, show_rownames = TRUE, treeheight_row = 0, treeheight_col = 70, fontsize_col = 12, annotation_names_col = F, annotation_col = annot_columns, filename = paste0("./results_outliers_included/Pairwise_sample_heatmap_", paste(plotgroups, collapse = "Vs"), ".tiff"))

Notes about heatmap colours.

You can change the colours used in the heatmap itself by changing the colour names (color=colorRampPalette....)

If you want to change the annotation colours, see here: https://zhiganglu.com/post/pheatmap_change_annotation_colors/

6. Differential expression analysis

In this section we use the Deseq2 package to identify differentially expressed genes.

USER INPUT

Choose the treatment groups you want to compare.

To see what groups are present, run the following:

unique(meta$group)

Enter which groups you want to compare (two groups only). BASELINE OR CONTROL GROUP SHOULD BE LISTED FIRST.

degroups <- c("normal", "Huntingtons_disease")

From the count table, pull out only the counts from the above groups

expdata <- as.matrix(counts[,meta$group %in% degroups])

Set up the experimental condition

'factor' sets up the reference level, i.e. which is the baseline group (otherwise the default baseline level is in alphabetic order)

condition <- factor(meta$group[meta$group %in% degroups], levels = degroups)

Type 'condition' in the console to see is the levels are set correctly

Set up column data (treatment groups and sample ID)

coldata <- data.frame(row.names=colnames(expdata), condition)

Create the DESeq2 dataset (dds)

dds <- DESeq2::DESeqDataSetFromMatrix(countData=expdata, colData=coldata, design=~condition)
dds$condition <- factor(dds$condition, levels = degroups)

Run DESeq2 to identify differentially expressed genes

deseq <- DESeq(dds)

Extract a results table from the DESeq analysis

res <- results(deseq)

Reorder results by adjusted p vales, so that the most signififcantly DE genes are at the top

res <- res[order(res$padj), ]

You can do a summary of the results to see how many significantly (alpha=0.05, adjust to 0.01 if needed) upregulated and downregulated DE genes were found

summary(res, alpha=0.05)

Convert from DESeq object to a data frame.

res <- data.frame(res)

Look at the top 6 DE genes

head(res)

Add normalised counts to the output table. This is so you can later plot expression trends for individual genes in R, Excel, etc.

Need to normalise the counts first, using the size factors calculated by DESeq2 (in the 'deseq' object)

expdata_norm <- as.matrix(expdata) %*% diag(deseq$sizeFactor)
colnames(expdata_norm) <- colnames(expdata)
annot_counts <- merge(x = res, y = expdata_norm, by = 0, all = TRUE)

Pull out just significant genes (change from 0.05 to 0.01 if needed)

DE_genes <- subset(annot_counts, padj < 0.05, select=colnames(annot_counts))

Export as a csv table

write.csv(DE_genes, file=paste0("./results_outliers_included/DE_genes_", paste(degroups, collapse = "Vs"), ".csv"), row.names = FALSE)

6b. Volcano plot

p <- EnhancedVolcano(res, lab = row.names(res), selectLab = row.names(res)[1:20], drawConnectors = TRUE, title = NULL, subtitle = NULL, x = 'log2FoldChange', y = 'pvalue')

p <- EnhancedVolcano(res, lab = rownames(res), pointSize = 3, drawConnectors = TRUE, title = NULL, subtitle = NULL, x = 'log2FoldChange', y = 'pvalue')

p

NOTE: the above plot shows labels for the top significantly DE (i.e. by lowest adjusted p value) genes.

Output as publication quality (300dpi) tiff and pdf.

Create a (300dpi) tiff

ggsave(file = paste0("./results_outliers_included/volcano_", paste(degroups, collapse = "Vs"), ".tiff"), dpi = 300, compression = "lzw", device = "tiff", width = 10, height = 8, plot = p)

Create a pdf

ggsave(file = paste0("./results_outliers_included/volcano_", paste(degroups, collapse = "Vs"), ".pdf"), device = "pdf", width = 10, height = 8, plot = p)

6c. DE genes heatmaps and dendrograms

sort by p-value

DE_genes <- DE_genes[order(DE_genes$padj), ]
row.names(DE_genes) <- DE_genes$Row.names

Pull out normalised counts only

siggc <- DE_genes[colnames(DE_genes) %in% colnames(expdata)]

Scale and center each row. This is important to visualise relative differences between groups and not have row-wise colouration dominated by high or low gene expression.

xts <- scale(t(siggc))
xtst <- t(xts)

Define annotation column

annot_columns <- data.frame(meta$group[meta$group %in% degroups])

Make the row names the sample IDs

row.names(annot_columns) <- meta$sample_ID[meta$group %in% degroups]
colnames(annot_columns) <- "Treatment groups"

Need to factorise it

annot_columns[[1]] <- factor(annot_columns[[1]])

Generate dendrogram and heatmap for ALL DE genes

pheatmap(xtst, color=colorRampPalette(c("#D55E00", "white", "#0072B2"))(100), annotation_col=annot_columns, annotation_names_col = F, fontsize_col = 12, fontsize_row = 7, labels_row = row.names(siggc), show_rownames = T, filename = paste0("./results_outliers_included/All_DEG_Heatmap_", paste(plotgroups, collapse = "Vs"), ".tiff"))

OUTLIER REMOVAL

This section repeats the above, but removes outliers first

REMOVE OUTLIERS FROM METADATA TABLE AN COUNT TABLE.

meta <- meta[- grep(c("WT4"), meta$sample_ID),]
counts <- counts[- grep(c("WT4"), colnames(counts))]

5. Outliers and batch effects

This section normalises and transforms the count data so that it can be plotted on a PCA plot and a heatmap

USER INPUT

Choose the groups you want to plot in a PCA/Heatmap. You can select any 2 or more of the groups (or all of the groups) you have in your 'groups' column of your metadata table

To see what groups are present, run the following:

unique(meta$group)

Now add which groups you want to plot (i.e. replace the groupnames below, and add more, separated by a comma and in "quotes", as needed). NOTE: R is case-sensitive, so these group names must be named EXACTLY the same as in the metadata table.

plotgroups <- c("normal", "Huntingtons_disease")

Pull out only the counts from the above groups

groupcounts <- counts[meta$group %in% plotgroups]

Normalise counts by library size, using DeSeq2's estimateSizeFactors() function. Note that DeSeq2 does this internally during DEG calling. The normalisation below is done separately for PCA and density plotting.

Set up the initial DeSeq2 experimental parameters.

condition <- factor(1:length(groupcounts))

Set up the column data. A data frame of sample ID's and conditions

coldata <- data.frame(row.names=colnames(groupcounts), condition)

Set up the DeSeq2 data set structure

f <- DESeqDataSetFromMatrix(countData = groupcounts, colData = coldata, design= ~ condition)

Estimate the size factors. See DeSeq2 manual for details

f <- estimateSizeFactors(f)

Size factors can be viewed by: sizeFactors(f)

Multiply each row (sample) by the corresponding size factor

subcount_norm <- as.matrix(groupcounts) %*% diag(sizeFactors(f))

Re-add column names

colnames(subcount_norm) <- colnames(groupcounts)

Remove low coverage transcripts (mean count < 10)

Find the mean of each row (and output as a data frame object)

means <- as.data.frame(rowMeans(subcount_norm))

Then join the means data with the counts

means <- cbind(means, subcount_norm)

Then subset out only genes with mean > 10

data <- subset(means, means[ , 1] > 10)

Remove the means column

data <- data[,-1]

Transform data

data_log <- vst(round(as.matrix(data)), nsub = nrow(data)-20)

Transformation can create some infinite values. Can't generate PCA data on these. Can see how many by: sum(sapply(data_log, is.infinite))

To remove infinite rows, use 'is.finte' or '!is.infinite'

data_log <- data_log[is.finite(rowSums(data_log)),]
colnames(data_log) <- colnames(groupcounts)

Set up the PCA plot base data

We're using the FactoMineR package to generate PCA plots (http://factominer.free.fr/index.html)

Need to transpose the data first

data_log_t <- t(data_log)

Add the group data

data_log_t_vars <- data.frame(meta$group[meta$group %in% plotgroups], data_log_t)

Generate the PCA data using FactoMineR package

res.pca <- PCA(data_log_t_vars, quali.sup = 1, graph=FALSE)

Set up the dendogram/heatmaps base data

Calculate the distance matrix:

distance_matrix <- as.matrix(dist(t(data_log)))

5a. PCA plot

Generate the PCA plot. Groups are shaded with ellipses at 95% confidence level. NOTE: at least 4 replicates need to be in a group for an ellipses to be drawn.

NOTE: change the group point colours by changing 'palette = ' below. Use the 'RColourBrewer' colour names (https://r-graph-gallery.com/38-rcolorbrewers-palettes.html). For example, if you are plotting 3 groups and choose palette = "Set1", this will use the first 3 colours from the Set1 colour palette.

p <- fviz_pca_ind(res.pca,
geom.ind = c("point", "text"), # show points only (but not "text")
col.ind = meta$group[meta$group %in% plotgroups], # color by groups
pointsize = 5, label = "all", title = "", legend.title = "Treatment groups", palette = "Dark2",
addEllipses = TRUE, ellipse.type = "t", ellipse.level = 0.95) + theme(legend.text = element_text(size = 12), legend.title = element_text(size = 14), axis.title=element_text(size=16), axis.text=element_text(size=14))

p

Output as publication quality (300dpi) tiff and pdf.

This will name your output files with the treatment groups you selected.

Create a 'results_outliers_removed' subdirectory where all results_outliers_removed will be output

dir.create("results_outliers_removed", showWarnings = FALSE)

Create a (300dpi) tiff

ggsave(file = paste0("./results_outliers_removed/PCA_", paste(plotgroups, collapse = "Vs"), ".tiff"), dpi = 300, compression = "lzw", device = "tiff", width = 10, height = 8, plot = p)

Create a pdf

ggsave(file = paste0("./results_outliers_removed/PCA_", paste(plotgroups, collapse = "Vs"), ".pdf"), device = "pdf", width = 10, height = 8, plot = p)

5b. Samples heatmap and dendrogram

This section plots a heatmap and dendrogram of pairwise relationships between samples. In this way you can see if samples cluster by treatment group.

See here: https://davetang.org/muse/2018/05/15/making-a-heatmap-in-r-with-the-pheatmap-package/

Define annotation column

annot_columns <- data.frame(meta$group[meta$group %in% plotgroups])

Make the row names the sample IDs

row.names(annot_columns) <- meta$sample_ID[meta$group %in% plotgroups]
colnames(annot_columns) <- "Treatment groups"

Need to factorise it

annot_columns[[1]] <- factor(annot_columns[[1]])

Generate dendrogram and heatmap

pheatmap(distance_matrix, color=colorRampPalette(c("white", "#9999FF", "#990000"))(50), cluster_rows = TRUE, show_rownames = TRUE, treeheight_row = 0, treeheight_col = 70, fontsize_col = 12, annotation_names_col = F, annotation_col = annot_columns, filename = paste0("./results_outliers_removed/Pairwise_sample_heatmap_", paste(plotgroups, collapse = "Vs"), ".tiff"))

Notes about heatmap colours.

You can change the colours used in the heatmap itself by changing the colour names (color=colorRampPalette....)

If you want to change the annotation colours, see here: https://zhiganglu.com/post/pheatmap_change_annotation_colors/

6. Differential expression analysis

In this section we use the Deseq2 package to identify differentially expressed genes.

USER INPUT

Choose the treatment groups you want to compare.

To see what groups are present, run the following:

unique(meta$group)

Enter which groups you want to compare (two groups only). BASELINE OR CONTROL GROUP SHOULD BE LISTED FIRST.

degroups <- c("normal", "Huntingtons_disease")

From the count table, pull out only the counts from the above groups

expdata <- as.matrix(counts[,meta$group %in% degroups])

Set up the experimental condition

'factor' sets up the reference level, i.e. which is the baseline group (otherwise the default baseline level is in alphabetic order)

condition <- factor(meta$group[meta$group %in% degroups], levels = degroups)

Type 'condition' in the console to see is the levels are set correctly

Set up column data (treatment groups and sample ID)

coldata <- data.frame(row.names=colnames(expdata), condition)

Create the DESeq2 dataset (dds)

dds <- DESeq2::DESeqDataSetFromMatrix(countData=expdata, colData=coldata, design=~condition)
dds$condition <- factor(dds$condition, levels = degroups)

Run DESeq2 to identify differentially expressed genes

deseq <- DESeq(dds)

Extract a results table from the DESeq analysis

res <- results(deseq)

Reorder results by adjusted p vales, so that the most signififcantly DE genes are at the top

res <- res[order(res$padj), ]

You can do a summary of the results to see how many significantly (alpha=0.05, adjust to 0.01 if needed) upregulated and downregulated DE genes were found

summary(res, alpha=0.05)

Convert from DESeq object to a data frame.

res <- data.frame(res)

Look at the top 6 DE genes

head(res)

Add normalised counts to the output table. This is so you can later plot expression trends for individual genes in R, Excel, etc.

Need to normalise the counts first, using the size factors calculated by DESeq2 (in the 'deseq' object)

expdata_norm <- as.matrix(expdata) %*% diag(deseq$sizeFactor)
colnames(expdata_norm) <- colnames(expdata)
annot_counts <- merge(x = res, y = expdata_norm, by = 0, all = TRUE)

Pull out just significant genes (change from 0.05 to 0.01 if needed)

DE_genes <- subset(annot_counts, padj < 0.05, select=colnames(annot_counts))

Export as a csv table

write.csv(DE_genes, file=paste0("./results_outliers_removed/DE_genes_", paste(degroups, collapse = "Vs"), ".csv"), row.names = FALSE)

6b. Volcano plot

p <- EnhancedVolcano(res, lab = row.names(res), selectLab = row.names(res)[1:20], drawConnectors = TRUE, title = NULL, subtitle = NULL, x = 'log2FoldChange', y = 'pvalue')

p <- EnhancedVolcano(res, lab = rownames(res), pointSize = 3, drawConnectors = TRUE, title = NULL, subtitle = NULL, x = 'log2FoldChange', y = 'pvalue')

p

NOTE: the above plot shows labels for the top significantly DE (i.e. by lowest adjusted p value) genes.

Output as publication quality (300dpi) tiff and pdf.

Create a (300dpi) tiff

ggsave(file = paste0("./results_outliers_removed/volcano_", paste(degroups, collapse = "Vs"), ".tiff"), dpi = 300, compression = "lzw", device = "tiff", width = 10, height = 8, plot = p)

Create a pdf

ggsave(file = paste0("./results_outliers_removed/volcano_", paste(degroups, collapse = "Vs"), ".pdf"), device = "pdf", width = 10, height = 8, plot = p)

6c. DE genes heatmaps and dendrograms

sort by p-value

DE_genes <- DE_genes[order(DE_genes$padj), ]
row.names(DE_genes) <- DE_genes$Row.names

Pull out normalised counts only

siggc <- DE_genes[colnames(DE_genes) %in% colnames(expdata)]

Scale and center each row. This is important to visualise relative differences between groups and not have row-wise colouration dominated by high or low gene expression.

xts <- scale(t(siggc))
xtst <- t(xts)

Define annotation column

annot_columns <- data.frame(meta$group[meta$group %in% degroups])

Make the row names the sample IDs

row.names(annot_columns) <- meta$sample_ID[meta$group %in% degroups]
colnames(annot_columns) <- "Treatment groups"

Need to factorise it

annot_columns[[1]] <- factor(annot_columns[[1]])

Generate dendrogram and heatmap for ALL DE genes

pheatmap(xtst, color=colorRampPalette(c("#D55E00", "white", "#0072B2"))(100), annotation_col=annot_columns, annotation_names_col = F, fontsize_col = 12, fontsize_row = 7, labels_row = row.names(siggc), show_rownames = T, filename = paste0("./results_outliers_removed/All_DEG_Heatmap_", paste(plotgroups, collapse = "Vs"), ".tiff"))