In step 1, we use the read_tsv() function in the readr package to load the tab delimited text file of counts into a dataframe called count_dataframe. Then, from that, we extract the 'gene' column to a new variable, genes, and erase it from count_dataframe, by assigning NULL. This is all done so we can easily convert into the count_matrix matrix with the base as.matrix() function and add the gene information back as rownames. Finally, we load the phenotype data we'll need from file using the readr read_table2() function.
Step 2 is concerned with working out which columns in count_matrix we want to use. We define a variable, experiments_of_interest, which holds the column names we want and then use the %in% operator and which() functions to create a binary vector that matches the number of columns. If, say, the third column of the columns_of_interest vector is TRUE it indicates the name was in the experiments_of interest variable.
Step 3 begins with loading the magrittr package to get the %>% operator, which will allow piping. We then use R indexing with the binary columns_of_interest factor to select the names of columns we want and send it to the forcats as_factor() function to get a factor object for our grouping variable. Sample grouping information is basically a factor that tells us which samples are replications of the same thing and it's important for the experimental design description. We need to create a grouping vector, each index of which refers to a column in the counts table. So, in the following example, the first three columns in the data would be replicates of one sample, the second three columns in the counts table would be replicates of a different replicate, and so on. We can use any symbols in the grouping vector to represent the groups. The more complicated the grouping vector, the more complicated the experiment design can be. In the recipe here, we'll use a simple test/control design:
numeric_groups <- c(1,1,1,2,2,2)
letter_groups <- c("A","A","A", "B","B","B")
A simple vector like this will do, but you can also use a factor object. The factor is R's categorical data type and is implemented as a vector of integers that have associated name labels, called levels. When a factor is displayed, the name labels are taken instead of the integers. The factor object has a memory of sorts, and even when a subset of levels is used, all of the levels that could have been used are retained so that when, for example, the levels are used as categories, empty levels can still be displayed.
In Step 4, we use indexing to extract the columns of data we want to actually analyze.
By Step 5, our preparatory work is done and we can build the DGEList object we need to do differential analysis. To start, we load the edgeR library and use the DGEList() function on counts_of_interest and our grouping object.
In Step 6, with DGEList, we can go through the edgeR process. First, we create the experimental design descriptor design object with the base model.matrix() function. A model design is required to tell the functions how to compare samples; this is a common thing in R and so has a base function. We use the grouping variable we created. We must estimate the dispersions of each gene with the estimateDisp() function, then we can use that measure of variability in tests. Finally, a generalized linear model is fit and the quasi-likelihood F-test is applied with the two uses of glmQLFTest(), first with the dispersal estimates, eset_dge, then with the resulting fit object.
We can use the topTags() function to see the details of differentially expressed genes. We get the following output:
## Coefficient: groupingL2Larvae
## logFC logCPM F PValue FDR
## FBgn0027527 6.318665 11.14876 42854.72 1.132951e-41 1.684584e-37
## [ reached 'max' / getOption("max.print") -- omitted 9 rows ]
The columns show the gene name, the logFC value of the gene, the F value, the P value and the False Detection Rate (FDR). Usually, the column we want to make statistical conclusions from is FDR.
