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"