Differential Expression Exercise - Revision history

Rf at 13:27, 10 May 2017

2017-05-10T13:27:51Z

Rf at 12:17, 10 May 2017

2017-05-10T12:17:55Z

Rf at 23:15, 9 May 2017

2017-05-09T23:15:08Z

Rf at 23:02, 9 May 2017

2017-05-09T23:02:24Z

Show changes

Rf at 22:16, 9 May 2017

2017-05-09T22:16:18Z

Show changes

Rf at 14:12, 4 May 2017

2017-05-04T14:12:06Z

Rf at 13:48, 4 May 2017

2017-05-04T13:48:14Z

Rf: Rf moved page I2rda dex to Differential Expression Exercise over redirect

2017-05-04T13:12:00Z

Rf moved page I2rda dex to Differential Expression Exercise over redirect

Rf: Rf moved page Differential Expression Exercise to I2rda dex

2017-05-04T12:55:50Z

Rf moved page Differential Expression Exercise to I2rda dex

Rf: s2

2017-04-08T22:08:44Z

@@ Line 376: / Line 376: @@
 Now quit R:
-  quit()
+  q()
 Have a look at the file pv_glm.txt:
@@ Line 384: / Line 384: @@
 Get rid of the first column, the first line and the quotation marks, replace the whitespace between the columns by a tab, and write the results to a file named <code>pv_glm.rnk</code>:
-This exercise has been done for you so you can skip it if you wish.
+This exercise has been done for you so you can skip it if you wish. We use the following old Unix tools: <code>cut</code> and <code>sed</code>. ALthough it is not necessary to understand is, you should be able to recognise the <code>s///g</code> subsitution idiom which <code>vim</code> and the Perl programming language also use.
   cut -f2,3 -d" " pv_glm.txt | sed '1d' | sed 's/\"//g' | sed 's/\s/\t/g' > pv_glm.rnk
@@ Line 400: / Line 400: @@
 == edgeR.R script ==
-So far we have been typing (single commands in a certain sequence.
+So far we have been typing (or copying pasting) single commands in a certain sequence, which is good and instructive for the first time. However, for later occasions, we would prefer to run that straight through, in other words turn the procedure into a script which can then be run as a single command, but this time, non-interactively. So the entire analysis can be run from the command line using the script <code>edgeR.R</code>:
   R --vanilla < edgeR.R
-== Extra tasks (optional) ==
+We can also run it from inside R with
 If you have time left, you can:
 # repeat the analysis increasing the minimum number of reads a gene must have in order to be included in the analysis.
-# repeat the analysis excluding patient A13 (the outlier in the MDS plot).
+# repeat the analysis excluding patient <code>A13</code> (the outlier in the MDS plot).
-The quickest way to do this is to change the appropriate command in the edgeR.R script and run it.
+The quickest way to do this is to change the appropriate command in the <code>edgeR.R</code> script and run it.
-To open the script with the nano text editor <code>nano edgeR.R</code> or with vim i.e. vim <code>edgeR.R</code>.
+Edit the with vim i.e. <code>vim edgeR.R</code>.
 You can now edit the script:
@@ Line 420: / Line 424: @@
 * For (2), change dgList2 in <code>> dgList2 <- dgList[, which(dgList$samples$patient != "A13")]</code> into dgList.
-When you're done, save via <code>Ctrl+X</code>, followed by <code>Y</code> and <code>Enter</code>.
+When you're done, save via <code>:sav! edgeR2.R</code> so that your edited script has a new name.
 <ins>Questions</ins>:
 * What is the effect of increasing the minimum number of reads on the number of differentially expressed genes?
 * Do you see any effect on the speed of the analysis?
-* What is the effect of excluding patient A13 on the number of differentially expressed genes?
+* What is the effect of excluding patient <code>A13</code> on the number of differentially expressed genes?

@@ Line 167: / Line 167: @@
 For closely located points, it is difficult to see the labels. However, we assigned the output of the <code>plotMDS</code> function to a variable called <code>mdp</code> and we might be able to see which are these samples.
-  > which(mdp$x>1.7 & mdp<1.85)
+  > which(mdp$x>1.7 & mdp$x<1.85)
 gives us the indices of the samples, so we can go back to dgList data frame and find out their names.
@@ Line 400: / Line 400: @@
 == edgeR.R script ==
 The entire analysis can be run from the command line using the script <code>edgeR.R</code>:
-  R –-vanilla < edgeR.R
+  R --vanilla < edgeR.R
 == Extra tasks (optional) ==

@@ Line 326: / Line 326: @@
 == Perform the differential expression calling ==
-Run a generalised linear model (GLM) likelihood ratio test on the data and write the results to a <code>DGEGLM</code> object named lrt:
+Run a generalised linear model (GLM) likelihood ratio test on the data and write the results to a <code>DGEGLM</code> object named <code>lrt</code>:
   > fit <- glmFit(dgGlm, design)
   > lrt <- glmLRT(fit)
-This tests by default for the last coefficient in the design matrix, i.e. the condition Cancer.
+This tests by default for the last coefficient in the design matrix, i.e. the condition "Cancer".
-Write the genes with a FDR < 0.05 to de_table_lrt:
+Write the genes with a FDR < 0.05 to <code>de_table_lrt</code>:
   > num <- length(lrt$table[,1])
@@ Line 340: / Line 340: @@
 <ins>Questions</ins>:
-* How many genes are found to be differentially expressed at 5% FDR?
+* How many genes are found to be differentially expressed at 5% FDR? (you need to remember the command).
 * Which method finds a larger number of differentially expressed genes, the group wise comparison or the matched pair comparison?

@@ Line 181: / Line 181: @@
 == Perform a simple two group comparison ==
 === Calculate the dispersion factors ===
@@ Line 232: / Line 234: @@
 * <code>logCPM</code>: average log2-counts per million
 * <code>PValue</code>: p-value
 Determine the number of genes in et and write the value to the variable num:
   > num <- length(et$table[,1])
-Rank the genes by p-value and write the results to a TopTags object named top:
+This uses the first column to calculate the total genes, though we could have chosen columns 2, 3, or 4 as well.
   > top <- topTags(et, n = num)
 <ins>top contains</ins>:
@@ Line 251: / Line 251: @@
 * <code>$table</code>: a data frame containing the logFC, logCPM, PValue and FDR (false discovery rate) for all genes
-Write the genes with a FDR < 0.05 to de_table:
+Let's have a look at the first ten genes that appear in top?
-  > de_table <- top[top$table[,4] < 0.05,]
+  > head(top$table, n=10)
-Look at the number of differentially expressed genes at 5% FDR:
+Say we're interested in genes with an FDR of less than 0.05. We can create a new TopTags data structure with just those genes in the following manner:
-  > nrow(de_table)
+ > delowfdr <- top[top$table[,4] < 0.05,]
 === Visualise the results ===
-  > detags <- rownames(de_table)
+A smear plot is almost identical to an MA plot, except some visual smoothing occurs when points appear together. It acts on the exactTest table, but will color certain points red if the genes names are supplied to it. Create a smear plot showing the fold change of genes against the log count per million for each gene:
   > plotSmear(et, de.tags = detags)
   > abline(h = c(-1,1), col = "blue")
@@ Line 272: / Line 277: @@
 * <code>col</code>: colour of lines
-What can you conclude from this plot?
+What can you say about this plot?
 === Write the results to a file ===

@@ Line 1: / Line 1: @@
 = Aims =
 In this part you will learn to:
 * perform differential expression analysis
@@ Line 5: / Line 6: @@
 You will use the following tools, which have been pre-installed on our bioinformatics training server at the University of St Andrews:
 * edgeR: http://www.bioconductor.org/packages/release/bioc/html/edgeR.html
 The data set you will investigate is from the study described in "RNA-Seq Analyses Generate Comprehensive Transcriptomic
-Landscape and Reveal Complex Transcript Patterns in Hepatocellular Carcinoma Data" by Huang et al. (http://journals.plos.org/plo
+The data set you will investigate is from the study described in "RNA-Seq Analyses Generate Comprehensive Transcriptomic Landscape and Reveal Complex Transcript Patterns in Hepatocellular Carcinoma Data" by Huang et al. (http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0026168). The raw data is available from ENA (http://www.ebi.ac.uk/ena/data/view/PRJNA134241).
-sone/article?id=10.1371/journal.pone.0026168). The raw data is available from ENA (http://www.ebi.ac.uk/ena/data/view/PRJNA13
-). The data set is from 26 samples of both cancerous and non-cancerous tissue from 10 patients. For two patients (A13 and
+The data set is from 26 samples of both cancerous and non-cancerous tissue from 10 patients. For two patients (A13 and
 A39) there is more than one cancer and non-cancer sample.
 You will use the following files:
@@ Line 36: / Line 37: @@
   less groups.tsv
-The file contains three tab-separated columns, that contain the name of the count file for a sample, the group to which the sample belongs (Canc
+The file contains three tab-separated columns, that contain the name of the count file for a sample, the group to which the sample belongs (Canc er or Non cancer) and the ID of the patient from which the sample was taken.
 = Differential expression analysis with edgeR =
@@ Line 56: / Line 56: @@
 == Read in the data ==
 Read in the data from the file groups.tsv and store them in a data frame named targets:
@@ Line 67: / Line 68: @@
 From the targets data frame, create a DGEList object named <code>dgList</code>:
-> dgList <- readDGE(targets, comment.char = "_")
+ > dgList <- readDGE(targets, comment.char = "_")
 You may get the following warning message, which you can ignore:
@@ Line 75: / Line 76: @@
   Name partially matched in data frame
-The function readDGE reads the count files, calculates the sizes of the count libraries (i.e. the total number of reads across all the genes in the
+The function readDGE reads the count files, calculates the sizes of the count libraries (i.e. the total number of reads across all the genes in the sample) and produces a DGEList object.
 sample) and produces a DGEList object.
 The argument comment.char = "_" is useful if you are using the output straight from htseq-count, as this tells readDGE to ignore the
 summary lines at the end of the file that start with __.
@@ Line 85: / Line 86: @@
 dgList contains:
-$samples: the meta data for the samples, which we provided in the file groups.tsv (files, groups and patient), plus some data
+* $samples: the meta data for the samples, which we provided in the file groups.tsv (files, groups and patient), plus some data calculated by edgeR, i.e the total number of reads across all the genes in the sample (lib.size), and the factor by which the library will be normalised (norm.factors), which is set at 1 at the moment.
-calculated by edgeR, i.e the total number of reads across all the genes in the sample (lib.size), and the factor by which the library will be
+* $counts: the read counts mapped to each gene for each of the 26 samples.
 By default, the levels of group (Cancer and Non cancer) are ordered in alphabetical order, with the first one considered the reference:
@@ Line 109: / Line 108: @@
 == Filter the data ==
 Keep only those genes that have at least 20 reads across all samples:
@@ Line 135: / Line 135: @@
 == Normalise the data ==
 Normalise the read counts based on library size using the Trimmed Mean of M-values (TMM) method:
@@ Line 147: / Line 148: @@
 Compare the library sizes and normalisation factors. Is there anything that strikes you?
-== Explore the data ==
+== Exploration of our Dataset ==
 Create an MDS (Multi Dimensional Scaling) plot for dgList:
-  > plotMDS(dgList, col = as.integer(dgList$samples$group), labels = dgList$samples$patient)
+  > mdp <- plotMDS(dgList, col = as.integer(dgList$samples$group), labels = dgList$samples$patient)
 <ins>where</ins>:
-dgList: data object
+* dgList: data object
-col: numeric or character vector of colors for the plotting characters
+* col: numeric or character vector of colors for the plotting characters
-as.integer: convert data type to integer (the group levels Non cancer and Cancer are of the data type factor, to base the colour of the data
+* as.integer: convert data type to integer (the group levels Non cancer and Cancer are of the data type factor, to base the colour of the data points on these they have to be converted to the data type integer)
-points on these they have to be converted to the data type integer)
+* labels: character vector of sample names or labels
 <ins>Questions</ins>:
 * Do multiple samples for the same condition from the same patient (A13 and A39) cluster together?
 * Do the cancer and non cancer samples separate nicely? Are there any outliers?
 Outliers can be removed from the analysis by creating a new DGEList object that doesn't contain the sample(s) in question. To create a new DGEList object named dgList2, that doesn't contain patient A13:
@@ Line 172: / Line 183: @@
 === Calculate the dispersion factors ===
 Estimate a single dispersion factor across all genes:
@@ Line 192: / Line 204: @@
 Have a look at dgList again:
-> dgList
+ > dgList
 The following information has been added to dgList:
@@ Line 204: / Line 216: @@
 This plots the biological coefficient of variation (square root of tagwise dispersions) against abundance (log2 counts per million) for each gene.
 What can you conclude form this plot?
 === Perform the differential expression calling ===
 Run an exact test on the data and write the results to a DGEExact object named et:
@@ Line 381: / Line 395: @@
 == edgeR.R script ==
 The entire analysis can be run from the command line using the script <code>edgeR.R</code>:
@@ Line 386: / Line 401: @@
 == Extra tasks (optional) ==
 If you have time left, you can:

← Older revision	Revision as of 13:12, 4 May 2017
(No difference)

← Older revision	Revision as of 12:55, 4 May 2017
(No difference)