Example of using curatedAdipoChIP
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 differentiation stages 0 and 1
for the transcription factor (CEBPB) and treated them
as independent groups. The goal of this example is to show how a typical
differential binding analysis can be applied in the dataset. The main
focus is to explain how the the data and metadata in
peak_counts fit in each main piece of the analysis. We
started by sub-setting the object to the samples of interest, filtering
the low quality samples and low count genes. Then we applied the
DESeq2 method with the default values.
Subsetting the object to samples and peaks of interest
First, we subset the peak_counts object to all samples
in the differentiation stage 0 or 1 for the transcription factor CEBPB.
Using the samples IDs we get chose only the features/peaks that was
found in at least one of samples. Notice that other criteria for
choosing peaks can be applied depending on the purpose of the analysis.
The total number of samples is 8 with 4 samples in each stage. The total
number of peaks is 49028.
# subset peaks_count object
# select samples for the factor CEBPB and stage 0 and 1
sample_ind <- (peak_counts$factor == 'CEBPB') & (peak_counts$stage %in% c(0, 1))
sample_ids <- colnames(peak_counts)[sample_ind]
# select peaks from the selected samples
peak_ind <- lapply(mcols(peak_counts)$peak, function(x) sum(sample_ids %in% x))
peak_ind <- unlist(peak_ind) > 2
# subset the object
se <- peak_counts[peak_ind, sample_ind]
# show the number of samples in each group
table(se$stage)
#>
#> 0 1
#> 4 4
# show the number of peak in each sample
table(unlist(mcols(se)$peak))[sample_ids]
#>
#> GSM2257705 GSM2257706 GSM544720 GSM544721 GSM686970 GSM686971 GSM686972
#> 749 677 1806 14251 14074 15381 7653
#> GSM686973
#> 3367Filtering low quality samples
Since 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. This kind of information is available in
colData.
# check the number of files in qc
qc <- se$qc
table(lengths(qc))
#>
#> 1 2
#> 6 2
# flattening qc list
qc <- unlist(qc, recursive = FALSE)
length(qc)
#> [1] 10The qc object of the metadata 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 %>%
dplyr::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.
# a quick look at quality scores
summary(per_base)
#> run Base Mean
#> Length:362 Min. : 1.0 Min. :28.24
#> Class :character 1st Qu.:10.0 1st Qu.:31.90
#> Mode :character Median :19.0 Median :34.80
#> Mean :18.6 Mean :33.96
#> 3rd Qu.:28.0 3rd Qu.:35.57
#> Max. :37.0 Max. :36.50To identify the low quality samples, we categorize the runs
run_average which correspond to the average of the per base
mean scores. The following figure should make it easier to see why this
cutoff were used in this case.
# find low quality runs
per_base <- per_base %>%
group_by(run) %>%
mutate(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() +
lims(y = c(0,40))
The run IDs of the “bad” samples is then used to remove them from the dataset.
# get run ids of low quality samples
bad_runs <- unique(per_base$run[per_base$run_average == FALSE])
bad_samples <- lapply(se$runs, function(x) sum(x$run %in% bad_runs))
# subset the counts object
se2 <- se[, unlist(bad_samples) == 0]
# show remaining samples by stage
table(se2$stage)
#>
#> 0 1
#> 3 3Filtering peaks with low counts
To identify the low count feature/peaks, we keep only the features with at least 10 reads in 2 or more samples. Then we subset the object to exclude the rest.
Check sample reproducibility
One expects that samples from the same group should reflect similar biology. In fact, this is also expected for the differential binding analysis to be reliable. In the context of ChIP-Seq, several measures were proposed to quantify the similarities and discrepancies between the samples, check this article for more. Here, we use scatter plots of the pairs of samples in each group to show that there is a general agreement in terms of the counts in the peakset.
# plot scatters of samples from each group
par(mfrow = c(2,3))
lapply(split(colnames(se3), se3$stage), function(x) {
# get counts of three samples
y1 <- assay(se3)[, x[1]]
y2 <- assay(se3)[, x[2]]
y3 <- assay(se3)[, x[3]]
# plot scatters of pairs of samples
plot(log10(y1 + 1), log10(y2 + 1), xlab = x[1], ylab = x[2])
plot(log10(y1 + 1), log10(y3 + 1), xlab = x[1], ylab = x[3])
plot(log10(y2 + 1), log10(y3 + 1), xlab = x[2], ylab = x[3])
})
#> $`0`
#> NULL
#>
#> $`1`
#> NULLApplying differential binding using DESeq2
DESeq2
is a well documented and widely used R package for the differential
binding analysis of ChIP-Seq data. Here we use the default values of
DESeq to find the peaks which are deferentially bound in
stage 1 compared to 0.
# differential binding analysis
colData(se3) <- colData(se3)[, -2]
se3$stage <- factor(se3$stage)
dds <- DESeqDataSet(se3, ~stage)
#> 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
res <- results(dds)
table(res$padj < .1)
#>
#> FALSE TRUE
#> 13297 1240Next!
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 stage 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 stage 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 = 'stage')
#> 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.