{"type":"doc","content":[{"type":"heading","attrs":{"level":2},"content":[{"text":"Public small RNA-seq data","type":"text"}]},{"type":"table","attrs":{"layout":"default","width":1800.0,"localId":"75a260a8-078b-496f-a0a8-7a1740110402"},"content":[{"type":"tableRow","content":[{"type":"tableHeader","attrs":{"colspan":1,"rowspan":1,"colwidth":[109.0]},"content":[{"type":"paragraph","content":[{"text":"Species","type":"text","marks":[{"type":"strong"}]}]}]},{"type":"tableHeader","attrs":{"colspan":1,"rowspan":1,"colwidth":[173.0]},"content":[{"type":"paragraph","content":[{"text":"ENA link","type":"text","marks":[{"type":"strong"}]}]}]},{"type":"tableHeader","attrs":{"colspan":1,"rowspan":1,"colwidth":[1518.0]},"content":[{"type":"paragraph","content":[{"text":"Description","type":"text","marks":[{"type":"strong"}]}]}]}]},{"type":"tableRow","content":[{"type":"tableCell","attrs":{"colspan":1,"rowspan":1,"colwidth":[109.0]},"content":[{"type":"paragraph","content":[{"text":"Human","type":"text"}]}]},{"type":"tableCell","attrs":{"colspan":1,"rowspan":1,"colwidth":[173.0]},"content":[{"type":"paragraph","content":[{"type":"inlineCard","attrs":{"url":"https://www.ebi.ac.uk/ena/browser/view/PRJEB5212?show=publications"}}]}]},{"type":"tableCell","attrs":{"colspan":1,"rowspan":1,"colwidth":[1518.0]},"content":[{"type":"paragraph","content":[{"text":"RNA-seq of micro RNAs (miRNAs) in Human prefrontal cortex to identify differentially expressed miRNAs between Huntington's Disease and control brain samples ","type":"text"}]}]}]}]},{"type":"heading","attrs":{"level":3},"content":[{"text":"1. Connect to an rVDI virtual desktop machine","type":"text"}]},{"type":"paragraph","content":[{"text":"To access and run an rVDI virtual desktop:","type":"text","marks":[{"type":"strong"}]}]},{"type":"paragraph","content":[{"text":"Go to ","type":"text"},{"text":"https://rvdi.qut.edu.au/","type":"text","marks":[{"type":"link","attrs":{"href":"https://rvdi.qut.edu.au/"}}]}]},{"type":"paragraph","content":[{"text":"Click on ‘","type":"text"},{"text":"VMware Horizon HTML Access","type":"text","marks":[{"type":"strong"}]},{"text":"’","type":"text"}]},{"type":"paragraph","content":[{"text":"Log on with your QUT username and password","type":"text"}]},{"type":"paragraph","content":[{"text":"*","type":"text"},{"text":"NOTE","type":"text","marks":[{"type":"strong"}]},{"text":": you need to be connected to the QUT network first, either being on campus or connecting remotely via VPN.","type":"text"}]},{"type":"heading","attrs":{"level":3},"content":[{"text":"2. Open PuTTY terminal","type":"text"}]},{"type":"bulletList","content":[{"type":"listItem","content":[{"type":"paragraph","content":[{"text":"Click on the PuTTY icon","type":"text"}]}]},{"type":"listItem","content":[{"type":"paragraph","content":[{"text":"Double-click on “Lyra”","type":"text"}]}]},{"type":"listItem","content":[{"type":"paragraph","content":[{"text":"Fill your password and connect to the HPC","type":"text"}]}]}]},{"type":"heading","attrs":{"level":3},"content":[{"text":"Precomputed results from session 6:","type":"text"}]},{"type":"paragraph","content":[{"text":"We ran the small RNA seq samples against the mirBase database and the results can be found at:","type":"text"}]},{"type":"codeBlock","content":[{"text":"/work/training/2024/smallRNAseq/runs/run2_human/results/mirna_quant/edger_qc/mature_counts.csv\n/work/training/2024/smallRNAseq/data/human_disease/metadata_microRNA.txt","type":"text"}]},{"type":"paragraph","content":[{"text":"The results of the miRNA profiling can be found in the folder called “","type":"text"},{"text":"mirna_quant/edger_qc","type":"text","marks":[{"type":"textColor","attrs":{"color":"#ff991f"}},{"type":"strong"}]},{"text":"”:","type":"text"}]},{"type":"codeBlock","content":[{"text":"├── results\n│ ├── bowtie_index\n│ ├── fastp\n│ ├── fastqc\n│ ├── mirna_quant\n│ ├── mirtrace\n| ├── multiqc\n| └── pipeline_info","type":"text"}]},{"type":"paragraph","content":[{"text":"inside the “","type":"text"},{"text":"mirna_quant/edger_qc","type":"text","marks":[{"type":"textColor","attrs":{"color":"#ff991f"}},{"type":"strong"}]},{"text":"” folder find the “","type":"text"},{"text":"mature_counts.csv","type":"text","marks":[{"type":"textColor","attrs":{"color":"#ff991f"}},{"type":"strong"}]},{"text":"” file:","type":"text"}]},{"type":"codeBlock","content":[{"text":"hairpin_counts.csv\nhairpin_CPM_heatmap.pdf\nhairpin_edgeR_MDS_distance_matrix.txt\nhairpin_edgeR_MDS_plot_coordinates.txt\nhairpin_edgeR_MDS_plot.pdf\nhairpin_log2CPM_sample_distances_dendrogram.pdf\nhairpin_log2CPM_sample_distances_heatmap.pdf\nhairpin_log2CPM_sample_distances.txt\nhairpin_logtpm.csv\nhairpin_logtpm.txt\nhairpin_normalized_CPM.txt\nhairpin_unmapped_read_counts.txt\nmature_counts.csv <-- we will use this file for the statistical analysis in the next section\nmature_counts.txt\nmature_CPM_heatmap.pdf\nmature_edgeR_MDS_distance_matrix.txt\nmature_edgeR_MDS_plot_coordinates.txt\nmature_edgeR_MDS_plot.pdf\nmature_log2CPM_sample_distances_dendrogram.pdf\nmature_log2CPM_sample_distances_heatmap.pdf\nmature_log2CPM_sample_distances.txt\nmature_logtpm.csv\nmature_logtpm.txt\nmature_normalized_CPM.txt\nmature_unmapped_read_counts.txt","type":"text"}]},{"type":"paragraph","content":[{"text":"Note: the “","type":"text","marks":[{"type":"strong"}]},{"text":"mature_counts.csv","type":"text","marks":[{"type":"textColor","attrs":{"color":"#ff991f"}},{"type":"strong"}]},{"text":"” needs to be transposed prior to running the statistical analysis. This can be done either user the R script or using a script called “","type":"text","marks":[{"type":"strong"}]},{"text":"transpose_csv.py","type":"text","marks":[{"type":"textColor","attrs":{"color":"#ff991f"}},{"type":"strong"}]},{"text":"”. ","type":"text","marks":[{"type":"strong"}]}]},{"type":"paragraph","content":[{"text":"Let’s initially create a “","type":"text"},{"text":"DESeq2","type":"text","marks":[{"type":"textColor","attrs":{"color":"#ff991f"}},{"type":"strong"}]},{"text":"” folder and ","type":"text"},{"text":"copy the files needed for the statistical analysis","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#d3f1a7"}}]},{"text":":","type":"text"}]},{"type":"codeBlock","content":[{"text":"mkdir -p $HOME/workshop/2024-2/session6_smallRNAseq/runs/run1_human_miRBase/DESeq2\nmkdir -p $HOME/workshop/2024-2/session6_smallRNAseq/runs/run2_human_MirGeneDB/DESeq2\nmkdir -p $HOME/workshop/2024-2/session6_smallRNAseq/runs/run3_drosha_miRBase/DESeq2\nmkdir -p $HOME/workshop/2024-2/session6_smallRNAseq/scripts\nmkdir -p $HOME/workshop/2024-2/session6_smallRNAseq/data\n\ncp /work/training/2024/smallRNAseq/scripts/* $HOME/workshop/2024-2/session6_smallRNAseq/scripts\ncp /work/training/2024/smallRNAseq/data/human_disease/* $HOME/workshop/2024-2/session6_smallRNAseq/data\n\ncp $HOME/workshop/2024-2/session6_smallRNAseq/scripts/transpose_csv.py $HOME/workshop/2024-2/session6_smallRNAseq/runs/run1_human_miRBase/DESeq2\ncp $HOME/workshop/2024-2/session6_smallRNAseq/data/metadata_microRNA.txt $HOME/workshop/2024-2/session6_smallRNAseq/runs/run1_human_miRBase/DESeq2\ncp /work/training/2024/smallRNAseq/runs/run2_human/results/mirna_quant/edger_qc/mature_counts.csv $HOME/workshop/2024-2/session6_smallRNAseq/runs/run1_human_miRBase/DESeq2\ncd $HOME/workshop/2024-2/session6_smallRNAseq/runs/run1_human_miRBase/DESeq2","type":"text"}]},{"type":"paragraph","content":[{"text":"To transpose the initial “","type":"text"},{"text":"mature_counts.csv","type":"text","marks":[{"type":"textColor","attrs":{"color":"#ff991f"}},{"type":"strong"}]},{"text":"” file do the following:","type":"text"}]},{"type":"codeBlock","content":[{"text":"python transpose_csv.py --input mature_counts.csv --out mature_counts.txt","type":"text"}]},{"type":"heading","attrs":{"level":2},"content":[{"text":"Differential expression analysis using RStudio","type":"text","marks":[{"type":"strong"}]}]},{"type":"paragraph","content":[{"text":"Differential expression analysis for smRNA-Seq is similar to regular RNA-Seq.","type":"text"}]},{"type":"paragraph","content":[{"text":"As with the previous RNA-Seq, we will also be running this smRNA-Seq differential expression analysis in RStudio on an rVDI virtual machine. The reason is the same as before - to save time as the required R packages are pre-installed on these virtual machines. And, as before, you can also copy and paste this script to RStudio on your local computer and adapt it to your own dataset.","type":"text"}]},{"type":"heading","attrs":{"level":3},"content":[{"text":"3. Run analysis script in RStudio","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"a. Open RStudio","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#d3f1a7"}}]}]},{"type":"paragraph","content":[{"text":"b. Create a new R script ('File'->'New File'-> ‘R script’)","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#d3f1a7"}}]}]},{"type":"paragraph","content":[{"text":"c. Hit the save button and save this file in the working directory you created above (H:\\workshop\\2024-2\\session6_smallRNAseq\\runs\\run1_human_miRBase\\DESeq2). Name the R script ‘DESeq2.R’.","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#d3f1a7"}}]}]},{"type":"paragraph","content":[{"text":"d. Copy and paste the entire script from the code window below into your R script.","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#d3f1a7"}}]}]},{"type":"paragraph","content":[{"text":"e. Run the entire script ('Code'-> ‘Run region’ → ‘Run all’)","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#d3f1a7"}}]}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"LOAD PACKAGES","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#d3f1a7"}}]}]},{"type":"codeBlock","content":[{"text":"#### 4. Loading required packages ####\n\n# This section needs to be run every time\n# Load packages\nbioconductor_packages <- c(\"DESeq2\", \"EnhancedVolcano\")\ncran_packages <- c(\"ggrepel\", \"ggplot2\", \"plyr\", \"reshape2\", \"FactoMineR\", \"factoextra\", \"pheatmap\")\nlapply(cran_packages, require, character.only = TRUE)\nlapply(bioconductor_packages, require, character.only = TRUE)\n","type":"text"}]},{"type":"paragraph","content":[{"text":"IMPORT DATA","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#d3f1a7"}}]}]},{"type":"codeBlock","content":[{"text":"#### 5. Import your count data ####\n\n# Set working directory. \n# Change this to your working directory\n# Set your home working directory\n# NOTE: \n# Working directory - \nsetwd(\"H:/workshop/2024-2/session6_smallRNAseq/runs/run1_human_miRBase/DESeq2\")\n\n# Import your count data. make sure you've created a 'data' subdirectory and put the count table file there.\nmetacounts <- read.csv(\"H:/workshop/2024-2/session6_smallRNAseq/runs/run1_human_miRBase/DESeq2/mature_counts.csv\", header = TRUE, row.names = 1)\n# Count table needs to be transposed and converted to a data frame\nmetacounts <- as.data.frame(t(metacounts))\n\n# Import metadata. Again, need a metadata_microRNA.txt file in the data subdirectory.\nmeta <- read.table(\"H:/workshop/2024-2/session6_smallRNAseq/runs/run1_human_miRBase/DESeq2/metadata_microRNA.txt\", header = TRUE)\n\n# Rename sample names to new sample IDs\ncounts <- metacounts[as.character(meta$sample_name)]\ncolnames(counts) <- meta$sample_ID\n","type":"text"}]},{"type":"paragraph","content":[{"text":"OUTLIERS AND BATCH EFFECTS - TRANSFORM DATA","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#d3f1a7"}}]}]},{"type":"paragraph","content":[{"text":"Q1. What line would you change if you wanted to change the threshold for low-coverage transcripts that you are filtering out?","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#dfd8fd"}}]}]},{"type":"codeBlock","content":[{"text":"#### 6. Outliers and batch effects ####\n\n# This section normalises and transforms the count data so that it can be plotted on a PCA plot and a heatmap\n\n## USER INPUT\n# 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\n# To see what groups are present, run the following:\nunique(meta$group)\n# 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.\nplotgroups <- c(\"normal\", \"Huntingtons_disease\")\n\n# Pull out only the counts from the above groups\ngroupcounts <- counts[meta$group %in% plotgroups]\n\n# 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.\n# Set up the initial DeSeq2 experimental parameters.\ncondition <- factor(1:length(groupcounts))\n# Set up the column data. A data frame of sample ID's and conditions\ncoldata <- data.frame(row.names=colnames(groupcounts), condition)\n# Set up the DeSeq2 data set structure\nf <- DESeqDataSetFromMatrix(countData = groupcounts, colData = coldata, design= ~ condition)\n# Estimate the size factors. See DeSeq2 manual for details\nf <- estimateSizeFactors(f)\n# Size factors can be viewed by: sizeFactors(f)\n\n# Multiply each row (sample) by the corresponding size factor\nsubcount_norm <- as.matrix(groupcounts) %*% diag(sizeFactors(f))\n# Re-add column names\ncolnames(subcount_norm) <- colnames(groupcounts)\n\n## Remove low coverage transcripts (mean count < 10) ##\n\n# Find the mean of each row (and output as a data frame object)\nmeans <- as.data.frame(rowMeans(subcount_norm))\n# Then join the means data with the counts\nmeans <- cbind(means, subcount_norm)\n# Then subset out only genes with mean > 10\ndata <- subset(means, means[ , 1] > 10)\n# Remove the means column\ndata <- data[,-1]\n\n# Transform data\ndata_log <- vst(round(as.matrix(data)), nsub = nrow(data)-20)\n# Transformation can create some infinite values. Can't generate PCA data on these. Can see how many by: sum(sapply(data_log, is.infinite))\n# To remove infinite rows, use 'is.finte' or '!is.infinite'\ndata_log <- data_log[is.finite(rowSums(data_log)),]\ncolnames(data_log) <- colnames(groupcounts)\n\n### Set up the PCA plot base data ###\n\n# We're using the FactoMineR package to generate PCA plots (http://factominer.free.fr/index.html)\n\n# Need to transpose the data first\ndata_log_t <- t(data_log)\n# Add the group data\ndata_log_t_vars <- data.frame(meta$group[meta$group %in% plotgroups], data_log_t)\n# Generate the PCA data using FactoMineR package\nres.pca <- PCA(data_log_t_vars, quali.sup = 1, graph=FALSE)\n\n## Set up the dendogram/heatmaps base data ##\n\n# Calculate the distance matrix:\ndistance_matrix <- as.matrix(dist(t(data_log)))\n","type":"text"}]},{"type":"paragraph","content":[{"text":"OUTLIERS AND BATCH EFFECTS - VISUALISE DATA (PCA)","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#d3f1a7"}}]}]},{"type":"paragraph","content":[{"text":"Q2. What line would you change if you wanted to show the 99% confidence level ellipses on the PCA?","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#dfd8fd"}}]}]},{"type":"codeBlock","content":[{"text":"#### 6a. PCA plot ####\n\n# 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.\n# 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.\np <- fviz_pca_ind(res.pca,\n geom.ind = c(\"point\", \"text\"), # show points only (but not \"text\")\n col.ind = meta$group[meta$group %in% plotgroups], # color by groups\n pointsize = 5, label = \"all\", title = \"\", legend.title = \"Treatment groups\", palette = \"Dark2\",\n 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))\n\np # If this doesn't show you a plot in the plot tab of the bottom right pane, then run this next line and try it again:\n#dev.off()\n\n# Output as publication quality (300dpi) tiff and pdf.\n# This will name your output files with the treatment groups you selected.\n\n# Create a 'results_outliers_included' subdirectory where all results_outliers_included will be output\ndir.create(\"results_outliers_included\", showWarnings = FALSE)\n# Create a (300dpi) tiff\nggsave(file = paste0(\"./results_outliers_included/PCA_\", paste(plotgroups, collapse = \"_Vs_\"), \".tiff\"), dpi = 300, compression = \"lzw\", device = \"tiff\", width = 10, height = 8, plot = p)\n# Create a pdf\nggsave(file = paste0(\"./results_outliers_included/PCA_\", paste(plotgroups, collapse = \"_Vs_\"), \".pdf\"), device = \"pdf\", width = 10, height = 8, plot = p)\n","type":"text"}]},{"type":"paragraph","content":[{"text":"OUTLIERS AND BATCH EFFECTS - VISUALISE DATA (HEATMAP)","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#d3f1a7"}}]}]},{"type":"paragraph","content":[{"text":"Q3. What line would you change if you wanted to change the colours of the heatmap purple to green?","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#dfd8fd"}}]}]},{"type":"codeBlock","content":[{"text":"#### 6b. Samples heatmap and dendrogram ####\n\n# This section plots a heatmap and dendrogram of pairwise relationships between samples. In this way you can see if samples cluster by treatment group.\n\n# See here: https://davetang.org/muse/2018/05/15/making-a-heatmap-in-r-with-the-pheatmap-package/\n\n# Define annotation column\nannot_columns <- data.frame(meta$group[meta$group %in% plotgroups])\n# Make the row names the sample IDs\nrow.names(annot_columns) <- meta$sample_ID[meta$group %in% plotgroups]\ncolnames(annot_columns) <- \"Treatment groups\"\n# Need to factorise it\nannot_columns[[1]] <- factor(annot_columns[[1]])\n\n# Generate dendrogram and heatmap\nprint(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\")))\n\n# Notes about heatmap colours.\n# You can change the colours used in the heatmap itself by changing the colour names (color=colorRampPalette....)\n# If you want to change the annotation colours, see here: https://zhiganglu.com/post/pheatmap_change_annotation_colors/\n","type":"text"}]},{"type":"paragraph","content":[{"text":"DIFFERENTIAL EXPRESSION ANALYSIS","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#d3f1a7"}}]}]},{"type":"paragraph","content":[{"text":"Q4. What line(s) would you change if you wanted to adjust the p-value of significantly differentially expresssed genes?","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#dfd8fd"}}]}]},{"type":"codeBlock","content":[{"text":"#### 7. Differential expression analysis ####\n\n# In this section we use the Deseq2 package to identify differentially expressed genes.\n\n## USER INPUT\n# Choose the treatment groups you want to compare. \n# To see what groups are present, run the following:\nunique(meta$group)\n# Enter which groups you want to compare (two groups only). BASELINE OR CONTROL GROUP SHOULD BE LISTED FIRST.\ndegroups <- c(\"normal\", \"Huntingtons_disease\")\n\n# From the count table, pull out only the counts from the above groups\nexpdata <- as.matrix(counts[,meta$group %in% degroups])\n\n# Set up the experimental condition\n# 'factor' sets up the reference level, i.e. which is the baseline group (otherwise the default baseline level is in alphabetic order)\ncondition <- factor(meta$group[meta$group %in% degroups], levels = degroups)\n# Type 'condition' in the console to see is the levels are set correctly\n\n# Set up column data (treatment groups and sample ID)\ncoldata <- data.frame(row.names=colnames(expdata), condition)\n\n# Create the DESeq2 dataset (dds)\ndds <- DESeq2::DESeqDataSetFromMatrix(countData=expdata, colData=coldata, design=~condition)\ndds$condition <- factor(dds$condition, levels = degroups)\n\n# Run DESeq2 to identify differentially expressed genes\ndeseq <- DESeq(dds)\n\n# Extract a results table from the DESeq analysis\nres <- results(deseq)\n\n# Reorder results by adjusted p vales, so that the most signififcantly DE genes are at the top\nres <- res[order(res$padj), ]\n\n# 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\nsummary(res, alpha=0.05)\n\n# Convert from DESeq object to a data frame.\nres <- data.frame(res)\n\n# Look at the top 6 DE genes\nhead(res)\n\n\n# Add normalised counts to the output table. This is so you can later plot expression trends for individual genes in R, Excel, etc.\n# Need to normalise the counts first, using the size factors calculated by DESeq2 (in the 'deseq' object)\nexpdata_norm <- as.matrix(expdata) %*% diag(deseq$sizeFactor)\ncolnames(expdata_norm) <- colnames(expdata)\nannot_counts <- merge(x = res, y = expdata_norm, by = 0, all = TRUE)\n\n\n# Pull out just significant genes (change from 0.05 to 0.01 if needed)\nDE_genes <- subset(annot_counts, padj < 0.05, select=colnames(annot_counts))\n\n# Export as a csv table\nwrite.csv(DE_genes, file=paste0(\"./results_outliers_included/DE_genes_\", paste(degroups, collapse = \"_Vs_\"), \".csv\"), row.names = FALSE)\n","type":"text"}]},{"type":"paragraph","content":[{"text":"DIFFERENTIAL EXPRESSION ANALYSIS - VISUALISATION (VOLCANO PLOT)","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#d3f1a7"}}]}]},{"type":"paragraph","content":[{"text":"Q5. What line would you change if you wanted to change the dot point size on the Volcano plot?","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#dfd8fd"}}]}]},{"type":"paragraph","content":[{"text":"Q6. What line would you change if you wanted to show different cutoff lines - vertically and horizontally? ","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#c6edfb"}}]},{"text":"https://bioconductor.org/packages/release/bioc/manuals/EnhancedVolcano/man/EnhancedVolcano.pdf","type":"text","marks":[{"type":"link","attrs":{"href":"https://bioconductor.org/packages/release/bioc/manuals/EnhancedVolcano/man/EnhancedVolcano.pdf"}}]}]},{"type":"codeBlock","content":[{"text":"#### 7b. Volcano plot ####\n\np <- EnhancedVolcano(res, lab = row.names(res), selectLab = row.names(res)[1:20], drawConnectors = TRUE, title = NULL, subtitle = NULL, x = 'log2FoldChange', y = 'pvalue')\n\np <- EnhancedVolcano(res, lab = rownames(res), pointSize = 3, drawConnectors = TRUE, title = NULL, subtitle = NULL, x = 'log2FoldChange', y = 'pvalue')\n\np\ndev.off()\n# NOTE: the above plot shows labels for the top significantly DE (i.e. by lowest adjusted p value) genes.\n\n# Output as publication quality (300dpi) tiff and pdf.\n# Create a (300dpi) tiff\nggsave(file = paste0(\"./results_outliers_included/volcano_\", paste(degroups, collapse = \"_Vs_\"), \".tiff\"), dpi = 300, compression = \"lzw\", device = \"tiff\", width = 10, height = 8, plot = p)\n# Create a pdf\nggsave(file = paste0(\"./results_outliers_included/volcano_\", paste(degroups, collapse = \"_Vs_\"), \".pdf\"), device = \"pdf\", width = 10, height = 8, plot = p)\n","type":"text"}]},{"type":"paragraph","content":[{"text":"DIFFERENTIAL EXPRESSION ANALYSIS - VISUALISATION (HEATMAP)","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#d3f1a7"}}]}]},{"type":"paragraph","content":[{"text":"Q7. What happens when you change the annotation_names_col to T?","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#dfd8fd"}}]}]},{"type":"codeBlock","content":[{"text":"#### 7c. DE genes heatmaps and dendrograms ####\n\n# sort by p-value\nDE_genes <- DE_genes[order(DE_genes$padj), ]\nrow.names(DE_genes) <- DE_genes$Row.names\n\n# Pull out normalised counts only\nsiggc <- DE_genes[colnames(DE_genes) %in% colnames(expdata)]\n\n# 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.\nxts <- scale(t(siggc))\nxtst <- t(xts)\n\n# Define annotation column\nannot_columns <- data.frame(meta$group[meta$group %in% degroups])\n# Make the row names the sample IDs\nrow.names(annot_columns) <- meta$sample_ID[meta$group %in% degroups]\ncolnames(annot_columns) <- \"Treatment groups\"\n# Need to factorise it\nannot_columns[[1]] <- factor(annot_columns[[1]])\n\n# Generate dendrogram and heatmap for ALL DE genes\nprint(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\")))\n\n","type":"text"}]},{"type":"paragraph","content":[{"text":"DIFFERENTIAL EXPRESSION ANALYSIS - OUTLIER REMOVAL AND REPEATING OF CODE ABOVE (you can run all this at once but it will overwrite the plots in the RStudio plot pane)","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#d3f1a7"}}]}]},{"type":"paragraph","content":[{"text":"Q8. Are there many outliers removed? How would you tell?","type":"text","marks":[{"type":"backgroundColor","attrs":{"color":"#dfd8fd"}}]}]},{"type":"codeBlock","content":[{"text":"\n#### OUTLIER REMOVAL ####\n\n# This section repeats the above, but removes outliers first\n\n# REMOVE OUTLIERS FROM METADATA TABLE AN COUNT TABLE.\nmeta <- meta[- grep(c(\"WT4\"), meta$sample_ID),]\ncounts <- counts[- grep(c(\"WT4\"), colnames(counts))]\n\n\n\n\n\n\n\n#### 8. Outliers and batch effects ####\n\n# This section normalises and transforms the count data so that it can be plotted on a PCA plot and a heatmap\n\n## USER INPUT\n# 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\n# To see what groups are present, run the following:\nunique(meta$group)\n# 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.\nplotgroups <- c(\"normal\", \"Huntingtons_disease\")\n\n# Pull out only the counts from the above groups\ngroupcounts <- counts[meta$group %in% plotgroups]\n\n# 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.\n# Set up the initial DeSeq2 experimental parameters.\ncondition <- factor(1:length(groupcounts))\n# Set up the column data. A data frame of sample ID's and conditions\ncoldata <- data.frame(row.names=colnames(groupcounts), condition)\n# Set up the DeSeq2 data set structure\nf <- DESeqDataSetFromMatrix(countData = groupcounts, colData = coldata, design= ~ condition)\n# Estimate the size factors. See DeSeq2 manual for details\nf <- estimateSizeFactors(f)\n# Size factors can be viewed by: sizeFactors(f)\n\n# Multiply each row (sample) by the corresponding size factor\nsubcount_norm <- as.matrix(groupcounts) %*% diag(sizeFactors(f))\n# Re-add column names\ncolnames(subcount_norm) <- colnames(groupcounts)\n\n## Remove low coverage transcripts (mean count < 10) ##\n\n# Find the mean of each row (and output as a data frame object)\nmeans <- as.data.frame(rowMeans(subcount_norm))\n# Then join the means data with the counts\nmeans <- cbind(means, subcount_norm)\n# Then subset out only genes with mean > 10\ndata <- subset(means, means[ , 1] > 10)\n# Remove the means column\ndata <- data[,-1]\n\n# Transform data\ndata_log <- vst(round(as.matrix(data)), nsub = nrow(data)-20)\n# Transformation can create some infinite values. Can't generate PCA data on these. Can see how many by: sum(sapply(data_log, is.infinite))\n# To remove infinite rows, use 'is.finte' or '!is.infinite'\ndata_log <- data_log[is.finite(rowSums(data_log)),]\ncolnames(data_log) <- colnames(groupcounts)\n\n### Set up the PCA plot base data ###\n\n# We're using the FactoMineR package to generate PCA plots (http://factominer.free.fr/index.html)\n\n# Need to transpose the data first\ndata_log_t <- t(data_log)\n# Add the group data\ndata_log_t_vars <- data.frame(meta$group[meta$group %in% plotgroups], data_log_t)\n# Generate the PCA data using FactoMineR package\nres.pca <- PCA(data_log_t_vars, quali.sup = 1, graph=FALSE)\n\n## Set up the dendogram/heatmaps base data ##\n\n# Calculate the distance matrix:\ndistance_matrix <- as.matrix(dist(t(data_log)))\n\n\n\n#### 8a. PCA plot ####\n\n# 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.\n# 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.\np <- fviz_pca_ind(res.pca,\n geom.ind = c(\"point\", \"text\"), # show points only (but not \"text\")\n col.ind = meta$group[meta$group %in% plotgroups], # color by groups\n pointsize = 5, label = \"all\", title = \"\", legend.title = \"Treatment groups\", palette = \"Dark2\",\n 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))\n\np\ndev.off()\n\n# Output as publication quality (300dpi) tiff and pdf.\n# This will name your output files with the treatment groups you selected.\n\n# Create a 'results_outliers_removed' subdirectory where all results_outliers_removed will be output\ndir.create(\"results_outliers_removed\", showWarnings = FALSE)\n# Create a (300dpi) tiff\nggsave(file = paste0(\"./results_outliers_removed/PCA_\", paste(plotgroups, collapse = \"_Vs_\"), \".tiff\"), dpi = 300, compression = \"lzw\", device = \"tiff\", width = 10, height = 8, plot = p)\n# Create a pdf\nggsave(file = paste0(\"./results_outliers_removed/PCA_\", paste(plotgroups, collapse = \"_Vs_\"), \".pdf\"), device = \"pdf\", width = 10, height = 8, plot = p)\n\n\n\n\n#### 8b. Samples heatmap and dendrogram ####\n\n# This section plots a heatmap and dendrogram of pairwise relationships between samples. In this way you can see if samples cluster by treatment group.\n\n# See here: https://davetang.org/muse/2018/05/15/making-a-heatmap-in-r-with-the-pheatmap-package/\n\n# Define annotation column\nannot_columns <- data.frame(meta$group[meta$group %in% plotgroups])\n# Make the row names the sample IDs\nrow.names(annot_columns) <- meta$sample_ID[meta$group %in% plotgroups]\ncolnames(annot_columns) <- \"Treatment groups\"\n# Need to factorise it\nannot_columns[[1]] <- factor(annot_columns[[1]])\n\n# Generate dendrogram and heatmap\nprint(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\")))\n\n# Notes about heatmap colours.\n# You can change the colours used in the heatmap itself by changing the colour names (color=colorRampPalette....)\n# If you want to change the annotation colours, see here: https://zhiganglu.com/post/pheatmap_change_annotation_colors/\n\n\n\n#### 9. Differential expression analysis ####\n\n# In this section we use the Deseq2 package to identify differentially expressed genes.\n\n## USER INPUT\n# Choose the treatment groups you want to compare. \n# To see what groups are present, run the following:\nunique(meta$group)\n# Enter which groups you want to compare (two groups only). BASELINE OR CONTROL GROUP SHOULD BE LISTED FIRST.\ndegroups <- c(\"normal\", \"Huntingtons_disease\")\n\n# From the count table, pull out only the counts from the above groups\nexpdata <- as.matrix(counts[,meta$group %in% degroups])\n\n# Set up the experimental condition\n# 'factor' sets up the reference level, i.e. which is the baseline group (otherwise the default baseline level is in alphabetic order)\ncondition <- factor(meta$group[meta$group %in% degroups], levels = degroups)\n# Type 'condition' in the console to see is the levels are set correctly\n\n# Set up column data (treatment groups and sample ID)\ncoldata <- data.frame(row.names=colnames(expdata), condition)\n\n# Create the DESeq2 dataset (dds)\ndds <- DESeq2::DESeqDataSetFromMatrix(countData=expdata, colData=coldata, design=~condition)\ndds$condition <- factor(dds$condition, levels = degroups)\n\n# Run DESeq2 to identify differentially expressed genes\ndeseq <- DESeq(dds)\n\n# Extract a results table from the DESeq analysis\nres <- results(deseq)\n\n# Reorder results by adjusted p vales, so that the most signififcantly DE genes are at the top\nres <- res[order(res$padj), ]\n\n# 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\nsummary(res, alpha=0.05)\n\n# Convert from DESeq object to a data frame.\nres <- data.frame(res)\n\n# Look at the top 6 DE genes\nhead(res)\n\n\n# Add normalised counts to the output table. This is so you can later plot expression trends for individual genes in R, Excel, etc.\n# Need to normalise the counts first, using the size factors calculated by DESeq2 (in the 'deseq' object)\nexpdata_norm <- as.matrix(expdata) %*% diag(deseq$sizeFactor)\ncolnames(expdata_norm) <- colnames(expdata)\nannot_counts <- merge(x = res, y = expdata_norm, by = 0, all = TRUE)\n\n\n# Pull out just significant genes (change from 0.05 to 0.01 if needed)\nDE_genes <- subset(annot_counts, padj < 0.05, select=colnames(annot_counts))\n\n# Export as a csv table\nwrite.csv(DE_genes, file=paste0(\"./results_outliers_removed/DE_genes_\", paste(degroups, collapse = \"_Vs_\"), \".csv\"), row.names = FALSE)\n\n\n#### 9b. Volcano plot ####\n\np <- EnhancedVolcano(res, lab = row.names(res), selectLab = row.names(res)[1:20], drawConnectors = TRUE, title = NULL, subtitle = NULL, x = 'log2FoldChange', y = 'pvalue')\n\np <- EnhancedVolcano(res, lab = rownames(res), pointSize = 3, drawConnectors = TRUE, title = NULL, subtitle = NULL, x = 'log2FoldChange', y = 'pvalue')\n\n\np\ndev.off()\n\n# NOTE: the above plot shows labels for the top significantly DE (i.e. by lowest adjusted p value) genes.\n\n# Output as publication quality (300dpi) tiff and pdf.\n# Create a (300dpi) tiff\nggsave(file = paste0(\"./results_outliers_removed/volcano_\", paste(degroups, collapse = \"_Vs_\"), \".tiff\"), dpi = 300, compression = \"lzw\", device = \"tiff\", width = 10, height = 8, plot = p)\n# Create a pdf\nggsave(file = paste0(\"./results_outliers_removed/volcano_\", paste(degroups, collapse = \"_Vs_\"), \".pdf\"), device = \"pdf\", width = 10, height = 8, plot = p)\n\n\n\n\n#### 9c. DE genes heatmaps and dendrograms ####\n\n# sort by p-value\nDE_genes <- DE_genes[order(DE_genes$padj), ]\nrow.names(DE_genes) <- DE_genes$Row.names\n\n# Pull out normalised counts only\nsiggc <- DE_genes[colnames(DE_genes) %in% colnames(expdata)]\n\n# 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.\nxts <- scale(t(siggc))\nxtst <- t(xts)\n\n# Define annotation column\nannot_columns <- data.frame(meta$group[meta$group %in% degroups])\n# Make the row names the sample IDs\nrow.names(annot_columns) <- meta$sample_ID[meta$group %in% degroups]\ncolnames(annot_columns) <- \"Treatment groups\"\n# Need to factorise it\nannot_columns[[1]] <- factor(annot_columns[[1]])\n\n# Generate dendrogram and heatmap for ALL DE genes\nprint(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\")))\n\n\n","type":"text"}]},{"type":"table","attrs":{"layout":"align-start","width":1800.0,"localId":"88992cd5-dc4e-45ee-8317-07e82785c5ce"},"content":[{"type":"tableRow","content":[{"type":"tableHeader","attrs":{"colspan":1,"rowspan":1,"colwidth":[883.0]},"content":[{"type":"mediaSingle","attrs":{"layout":"center","width":610,"widthType":"pixel"},"content":[{"type":"media","attrs":{"width":2099,"alt":"All_DEG_Heatmap_normal_Vs_Huntingtons_disease_outliers.tiff","id":"5bc7dc4c-0ba0-4733-9efa-12cc57fadefe","collection":"contentId-2406940692","type":"file","height":2099}},{"type":"caption","content":[{"text":"Figure 1. All DEG with outliers included","type":"text"}]}]}]},{"type":"tableHeader","attrs":{"colspan":1,"rowspan":1,"colwidth":[909.0]},"content":[{"type":"mediaSingle","attrs":{"layout":"center","width":614,"widthType":"pixel"},"content":[{"type":"media","attrs":{"width":2099,"alt":"All_DEG_Heatmap_normal_Vs_Huntingtons_disease.tiff","id":"9c4ca1bd-e52f-421b-88ec-a9f90f09d688","collection":"contentId-2406940692","type":"file","height":2099}},{"type":"caption","content":[{"text":"Figure 2. All DEG without outliers (i.e. removed)","type":"text"}]}]},{"type":"paragraph"}]}]}]},{"type":"paragraph","content":[{"text":" ","type":"text"}]},{"type":"table","attrs":{"layout":"align-start","width":1800.0,"localId":"678abc96-bf7d-4d28-9495-4c5a42e17dc3"},"content":[{"type":"tableRow","content":[{"type":"tableHeader","attrs":{"colspan":1,"rowspan":1,"colwidth":[881.0]},"content":[{"type":"mediaSingle","attrs":{"layout":"center","width":622,"widthType":"pixel"},"content":[{"type":"media","attrs":{"width":2100,"alt":"Pairwise_sample_heatmap_normal_Vs_Huntingtons_disease_outliers.tiff","id":"098e73bf-70b2-4699-894a-9a2e42e0f6e5","collection":"contentId-2406940692","type":"file","height":2100}},{"type":"caption","content":[{"text":"Figure 3. Heatmap with outliers included","type":"text"}]}]},{"type":"paragraph"}]},{"type":"tableHeader","attrs":{"colspan":1,"rowspan":1,"colwidth":[908.0]},"content":[{"type":"mediaSingle","attrs":{"layout":"center","width":614,"widthType":"pixel"},"content":[{"type":"media","attrs":{"width":2100,"alt":"Pairwise_sample_heatmap_normal_Vs_Huntingtons_disease.tiff","id":"7d9044da-e7a3-4a22-ba18-f09dd2ecdbe4","collection":"contentId-2406940692","type":"file","height":2100}},{"type":"caption","content":[{"text":"Figure 4. Heatmap without outliers (i.e. removed)","type":"text"}]}]},{"type":"paragraph"}]}]}]},{"type":"table","attrs":{"layout":"align-start","width":1800.0,"localId":"5f131fdc-0ebb-4078-b416-a29c490a6781"},"content":[{"type":"tableRow","content":[{"type":"tableHeader","attrs":{"colspan":1,"rowspan":1},"content":[{"type":"mediaSingle","attrs":{"layout":"center","width":616,"widthType":"pixel"},"content":[{"type":"media","attrs":{"width":3000,"alt":"PCA_normal_Vs_Huntingtons_disease_outliers.tiff","id":"628fa6ac-0e26-4d2c-917e-614682f03ac2","collection":"contentId-2406940692","type":"file","height":2400}},{"type":"caption","content":[{"text":"Figure 5. PCA with outliers included","type":"text"}]}]}]},{"type":"tableHeader","attrs":{"colspan":1,"rowspan":1},"content":[{"type":"mediaSingle","attrs":{"layout":"center","width":616,"widthType":"pixel"},"content":[{"type":"media","attrs":{"width":3000,"alt":"PCA_normal_Vs_Huntingtons_disease.tiff","id":"f8c96bfa-8abd-4f61-8b3f-a60ba4c0a8c2","collection":"contentId-2406940692","type":"file","height":2400}},{"type":"caption","content":[{"text":"Figure 6. PCA without outliers (i.e. removed)","type":"text"}]}]},{"type":"paragraph"}]}]}]},{"type":"table","attrs":{"layout":"align-start","width":1800.0,"localId":"018c56bd-3669-4248-957a-ceba9e44784a"},"content":[{"type":"tableRow","content":[{"type":"tableHeader","attrs":{"colspan":1,"rowspan":1},"content":[{"type":"mediaSingle","attrs":{"layout":"center","width":616,"widthType":"pixel"},"content":[{"type":"media","attrs":{"width":3000,"alt":"volcano_normal_Vs_Huntingtons_disease_outliers.tiff","id":"2d67f9fd-367b-4319-b582-0770d85b4a8b","collection":"contentId-2406940692","type":"file","height":2400}},{"type":"caption","content":[{"text":"Figure 7. Volcano with outliers included","type":"text"}]}]},{"type":"paragraph"}]},{"type":"tableHeader","attrs":{"colspan":1,"rowspan":1},"content":[{"type":"mediaSingle","attrs":{"layout":"center","width":619,"widthType":"pixel"},"content":[{"type":"media","attrs":{"width":3000,"alt":"volcano_normal_Vs_Huntingtons_disease.tiff","id":"70bee215-9ea6-4e78-8ce0-fb571479ee03","collection":"contentId-2406940692","type":"file","height":2400}},{"type":"caption","content":[{"text":"Figure 8. Volcano plot without outliers (i.e. removed)","type":"text"}]}]},{"type":"paragraph"}]}]}]},{"type":"layoutSection","content":[{"type":"layoutColumn","attrs":{"width":50.0},"content":[{"type":"paragraph","content":[{"text":"Introduction","type":"text","marks":[{"type":"link","attrs":{"href":"https://eresearchqut.atlassian.net/wiki/spaces/EG/pages/2364702806"}}]}]}]},{"type":"layoutColumn","attrs":{"width":50.0},"content":[{"type":"paragraph","marks":[{"type":"alignment","attrs":{"align":"end"}}],"content":[{"text":"Exercises","type":"text","marks":[{"type":"link","attrs":{"href":"https://eresearchqut.atlassian.net/wiki/spaces/EG/pages/2406711323"}}]}]}]}]}],"version":1}

Browser not supported