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.