Example of using curatedAdipoRNA
Motivation
All the samples in this dataset come from the 3T3-L1 cell line. The MDI induction media, were used to induce adipocyte differentiation. The two important variables in the dataset are time and stage, which correspond to the time point and stage of differentiation when the sample were captured. Ideally, this dataset should be treated as a time course. However, for the purposes of this example, we only used samples from two time points 0 and 24 hours and treated them as independent groups. The goal of this example is to show how a typical differential expression analysis can be applied in the dataset. The main focus is to explain how the the data and metadata in adipo_counts fit in each main piece of the analysis. We started by filtering the low quality samples and low count genes. Then we applied the DESeq2 method with the default values.
Filtering low quality samples
First, we subset the adipo_counts object to all samples that has time points 0 or 24. The total number of samples is 30; 22 at 0 hour and 8 samples at 24 hours. The total number of features/genes in the set is 23916.
# subsetting counts to 0 and 24 hours
se <- adipo_counts[, adipo_counts$time %in% c(0, 24)]
# showing the numbers of features, samples and time groups
dim(se)
#> [1] 23916 30
table(se$time)
#>
#> 0 24
#> 22 8Since the quality metrics are reported per run file, we need to get the SSR* id for each of the samples. Notice that, some samples would have more than one file. In this case because some of the samples are paired-end, so each of them would have two files SRR\*_1 and SRR\*_2.
# filtering low quality samples
# chek the library layout
table(se$library_layout)
#>
#> PAIRED SINGLE
#> 12 18
# check the number of files in qc
qc <- se$qc
table(lengths(qc))
#>
#> 1 2
#> 18 12
# flattening qc list
qc <- unlist(qc, recursive = FALSE)
length(qc)
#> [1] 42The qc object of the colData contains the output of fastqc in a qc_read class. More information on this object can be accessed by calling ?fastqcr::qc_read. Here, we only use the per_base_sequence_quality to filter out low quality samples. This is by no means enough quality control but it should drive the point home.
# extracting per_base_sequence_quality
per_base <- lapply(qc, function(x) {
df <- x[['per_base_sequence_quality']]
df %>%
select(Base, Mean) %>%
transform(Base = strsplit(as.character(Base), '-')) %>%
unnest(Base) %>%
mutate(Base = as.numeric(Base))
}) %>%
bind_rows(.id = 'run')After tidying the data, we get a data.frame with three columns; run, Mean and Base for the run ID, the mean quality score and the base number in each read. fastqc provide thorough documentation of this quality control module and others. Notice that read length varies significantly between the runs and that the average of the mean score is suitable.
# a quick look at quality scores
summary(per_base)
#> run Base Mean
#> Length:3408 Min. : 1.00 Min. :10.74
#> Class :character 1st Qu.: 21.00 1st Qu.:34.45
#> Mode :character Median : 42.00 Median :36.44
#> Mean : 50.24 Mean :35.28
#> 3rd Qu.: 70.00 3rd Qu.:37.90
#> Max. :368.00 Max. :40.17To identify the low quality samples, we categorize the runs by length and run_average which are the read length and the average of the per base mean scores. The following figure should make it easier to see why these cutoff were used in this case.
# find low quality runs
per_base <- per_base %>%
group_by(run) %>%
mutate(length = max(Base) > 150,
run_average = mean(Mean) > 34)
# plot average per base quality
per_base %>%
ggplot(aes(x = Base, y = Mean, group = run, color = run_average)) +
geom_line() +
facet_wrap(~length, scales = 'free_x')
The run IDs of the “bad” samples is then used to remove them from the dataset.
# get run ids of low quality samples
bad_samples <- data.frame(samples = unique(per_base$run[per_base$run_average == FALSE]))
bad_samples <- separate(bad_samples, col = samples, into = c('id', 'run'), sep = '\\.')
# subset the counts object
se2 <- se[, !se$id %in% bad_samples$id]
table(se2$time)
#>
#> 0 24
#> 19 6Filtering low count genes
To identify the low count feature/genes (possibly not expressed), we keep only the features with at least 10 reads in 2 or more samples. Then we subset the object to exclude these genes.
Applying differential expression using DESeq2
DESeq2 is a well documented and widely used R package for the differential expression analysis. Here we use the default values of DESeq to find the genes which are deferentially expressed between the samples at 24 hours and 0 hours.
# differential expression analysis
se3$time <- factor(se3$time)
dds <- DESeqDataSet(se3, ~time)
#> renaming the first element in assays to 'counts'
#> converting counts to integer mode
dds <- DESeq(dds)
#> estimating size factors
#> estimating dispersions
#> gene-wise dispersion estimates
#> mean-dispersion relationship
#> final dispersion estimates
#> fitting model and testing
#> -- replacing outliers and refitting for 53 genes
#> -- DESeq argument 'minReplicatesForReplace' = 7
#> -- original counts are preserved in counts(dds)
#> estimating dispersions
#> fitting model and testing
res <- results(dds)
table(res$padj < .1)
#>
#> FALSE TRUE
#> 6194 7462Next!
In this example, we didn’t attempt to correct for the between study factors that might confound the results. To show how is this possible, we use the PCA plots with a few of these factors in the following graphs. The first uses the time factor which is the factor of interest in this case. We see that the DESeq transformation did a good job separating the samples to their expected groups. However, it also seems that the time is not the only factor in play. For example, we show in the second and the third graphs two other factors library_layout and instrument_model which might explain some of the variance between the samples. This is expected because the data were collected from different studies using slightly different protocols and different sequencing machines. Therefore, it is necessary to account for these differences to obtain reliable results. There are multiple methods to do that such as Removing Unwanted Variation (RUV) and Surrogate Variable Analysis (SVA).
# explaining variabce
plotPCA(rlog(dds), intgroup = 'time')
#> using ntop=500 top features by variance
Citing the studies in this subset of the data
Speaking of studies, as mentioned earlier the studies object contains full information of the references of the original studies in which the data were published. Please cite them when using this dataset.
# keys of the studies in this subset of the data
unique(se3$bibtexkey)
#> [1] "Duteil2014" "zhao_fto-dependent_2014"
#> [3] "siersbaek_molecular_2014" "Lim2015"
#> [5] "brunmeir_comparative_2016" "Brier2017"
#> [7] "park_distinct_2017" "chen_diabetes_2018"
#> [9] "siersbaek_dynamic_2017" "ryu_metabolic_2018"