PiGx-RNAseq - DESeq2 Report
BIMSB Bioinformatics Platform
2018-07-20 11:39:00
1 Description
PiGx RNAseq performs differential expression analysis using DESeq2, and produces this report. The report includes tables and figures summarizing similarities and differences between comparison groups as specified in the settings file. In addition to the differential expression analysis, there are plots and statistics for quality control of the experiment in general with an emphasis on the reproducibility of the sequencing results among the biological replicates.
This report was generated with PiGx RNAseq version 0.0.4.
2 Input Settings
countDataFile <- params$countDataFile
colDataFile <- params$colDataFile
gtfFile <- params$gtfFile
caseSampleGroups <- params$caseSampleGroups
controlSampleGroups <- params$controlSampleGroups
covariates <- params$covariates
prefix <- params$prefix
workdir <- params$workdir
organism <- gsub('\\s*', '', params$organism)
#whether to do GO analysis or not
runGO <- TRUE
#if organism is not provided, it is not possible to do GO analysis
if(organism == '') {
runGO <- FALSE
}
#create a folder to save high quality images generated by the report script
imagesDir <- file.path(workdir, paste0(prefix, '_images'))
if(! dir.exists(imagesDir)) {
dir.create(path = imagesDir)
}
inputParameterDesc <- c('Count Data File',
'Experiment Data File',
'GTF File',
'Case sample groups',
'Control sample groups',
'Covariates to control for',
'Prefix for output files',
'Working directory',
'Analyzed organism'
)
inputParameterValues <- c(countDataFile,
colDataFile,
gtfFile,
caseSampleGroups,
controlSampleGroups,
covariates,
prefix,
workdir,
organism)
inputSettings <- data.frame(parameters = inputParameterDesc,
values = inputParameterValues,
stringsAsFactors = FALSE)
DT::datatable(data = inputSettings,
extensions = 'FixedColumns',
options = list(fixedColumns = TRUE,
scrollX = TRUE,
pageLength = 9,
dom = 't'))gtfData <- rtracklayer::import.gff(con = gtfFile, format = 'gtf')
caseSamples <- gsub(' ', '', unlist(strsplit(x = caseSampleGroups, split = ',')))
controlSamples <- gsub(' ', '', unlist(strsplit(x = controlSampleGroups, split = ',')))
covariates <- gsub(' ', '', unlist(strsplit(x = covariates, split = ',')))
#read colData and countData files
colData = read.table(colDataFile, header=T, row.names = 1, sep='\t', stringsAsFactors = T, check.names = FALSE)
countData = read.table(countDataFile, header=TRUE, row.names=1, sep='\t', stringsAsFactors = T, check.names = FALSE)
#subset colData and countData - only keep case and control samples
colData <- colData[colData$group %in% c(caseSamples, controlSamples),]
countData <- subset(countData, select = rownames(colData))
#split samples as case/control for deseq
colData$AnalysisGroup <- 'Control'
colData[colData$group %in% caseSamples,]$AnalysisGroup <- 'Case'mapIdsToNames <- function(ids, gtfData) {
#first figure out if the given ids are transcript or gene ids
transcripts <- gtfData[gtfData$type == 'transcript']
df <- unique(data.frame('transcript_id' = transcripts$transcript_id,
'gene_id' = transcripts$gene_id,
'gene_name' = transcripts$gene_name, stringsAsFactors = FALSE))
m <- apply(head(df[,1:2], 1000), 2, function(x) sum(x %in% ids))
#then map the ids to gene names
if(m['transcript_id'] > m['gene_id']){
return(df[match(ids, df$transcript_id),]$gene_name)
} else {
return(df[match(ids, df$gene_id),]$gene_name)
}
}
if(length(covariates) > 0){
designFormula <- paste("~", paste(covariates, collapse = ' + '), "+ AnalysisGroup")
} else {
designFormula <- "~ AnalysisGroup"
}
message("design formula:", designFormula)
dds <- DESeq2::DESeqDataSetFromMatrix(countData = countData, colData = colData, design = stats::as.formula(designFormula))
dds <- dds[ rowSums(counts(dds)) > 1, ]
dds <- DESeq2::DESeq(dds)
vsd <- DESeq2::varianceStabilizingTransformation(dds)
vsd.counts = SummarizedExperiment::assay(vsd)
DEtable = DESeq2::results(dds, contrast = c("AnalysisGroup", 'Case', 'Control'))
DEtable <- DEtable[order(DEtable$padj),]
DE <- as.data.frame(DEtable)
DE$geneName <- mapIdsToNames(rownames(DE), gtfData)
DEnormalizedCountsFile <- file.path(workdir, paste0(prefix, '.normalized_counts.tsv'))
write.table(x = vsd.counts,
file = DEnormalizedCountsFile,
quote = FALSE, sep = '\t')
DEresultsFile <- file.path(workdir, paste0(prefix, '.deseq_results.tsv'))
write.table(x = DE,
file = DEresultsFile,
quote = FALSE, sep = '\t')3 Differential Expression Analysis
Differential expression (DE) analysis was done using the DESeq2 R package. First, read counts are transformed using a variance stabilizing transformation, and then the expression values of each gene is compared between the control and sample groups using a negative binomial distribution as a model.
3.1 Differential Expression Results Table (top 1000)
This is the table of top 1000 differentially expressed genes comparing cases to controls (as specified in the input settings listed above). The table is first filtered by absolute log2foldChange > 1 and sorted by adjusted P value after multiple testing correction (padj). The baseMean refers to expression level in the controls, and the log fold change column denotes the expression in the cases, as compared to the control.
The full table of DESeq2 results and normalized counts tables can be found at:
- DESeq2 results table: /data/local/buyar/pigx/rnaseq/manuscript-use-case/results_resubmission/report/WT_diff_day6.star.deseq_results.tsv
- DESeq2 normalized counts table: /data/local/buyar/pigx/rnaseq/manuscript-use-case/results_resubmission/report/WT_diff_day6.star.normalized_counts.tsv
DEsubset <- DE[!is.na(DE$padj) & abs(DE$log2FoldChange) > 1,]
max <- 1000
if(nrow(DEsubset) < max) {
max <- nrow(DEsubset)
}
DEsubset <- DEsubset[1:max,]
DT::datatable(DEsubset,
extensions = c('Buttons', 'FixedColumns', 'Scroller'),
options = list(fixedColumns = TRUE,
scrollY = 400,
scrollX = TRUE,
scroller = TRUE,
dom = 'Bfrtip',
buttons = c('colvis', 'copy', 'print', 'csv','excel', 'pdf'),
columnDefs = list(
list(targets = c(3,4,5), visible = FALSE)
)),
filter = 'bottom'
)4 Diagnostic Plots
This section holds a number of plots meant for a quick diagnostic and/or sanity check of the analysis.
4.1 Number of reads assigned to genes
This plot shows the number of reads, in each sample, that are assigned to genes/transcripts. Outlier samples may be faulty and should be examined.
readCounts <- as.data.frame(colSums(countData))
readCounts$group <- colData[rownames(readCounts),]$group
readCounts$sample <- rownames(readCounts)
colnames(readCounts)[1] <- 'readCounts'
quantiles <- quantile(readCounts$readCounts, c(1:20)/20)[c(1,5,15,19)]
p <- ggplot(readCounts, aes(x = sample, y = readCounts)) + geom_bar(aes(fill = group), stat = 'identity') +
geom_hline(yintercept = as.numeric(quantiles), color = 'red') +
geom_label_repel(data = data.frame(x = 0, y = as.numeric(quantiles)), aes(x = x, y = y, label = names(quantiles))) + theme(legend.position = 'bottom') + scale_y_continuous(labels = scales::comma) + coord_flip()
print(p)
#save image to folder
pdf(file = file.path(imagesDir, 'readcounts.pdf'))
print(p)
invisible(dev.off())4.2 p-value histogram
The P value distribution from the DE analysis. The expected shape depends on the expected difference / similarity between the controls and samples.
p <- ggplot(data = DE, aes(x = pvalue)) + geom_histogram(bins = 100)
print(p)
#save image to folder
pdf(file = file.path(imagesDir, 'pvalue_histogram.pdf'))
print(p)
invisible(dev.off())4.3 MA plot
The MA plot gives an overview of the comparison between the two groups in the experiment. The log2 fold change for each gene is plotted on the y axis, against the average expression of that gene on the x axis.
DESeq2::plotMA(DEtable, main=paste("MA plot"))
#save image to folder
pdf(file = file.path(imagesDir, 'MA_plot.pdf'))
DESeq2::plotMA(DEtable, main=paste("MA plot"))
invisible(dev.off())plotGroups <- c(covariates, 'AnalysisGroup', 'group')
pcaPlots <- lapply(plotGroups, function(g) {
pca <- stats::prcomp(t(vsd.counts), center = TRUE)
pcaSummary <- summary(pca)
df <- merge(as.data.frame(pca$x), colData, by = 'row.names')
ggplot(df, aes(x = PC1, y = PC2)) +
geom_point(aes_string(color = g)) +
geom_label_repel(aes(label = Row.names), size = 3) +
labs(x = paste0('PC1 (',round(pcaSummary$importance[2, 'PC1'] * 100, 1),'%)'),
y = paste0('PC2 (',round(pcaSummary$importance[2, 'PC2'] * 100, 1),'%)')) +
theme_bw()
})
#save image to folder
pdf(file = file.path(imagesDir, 'pcaPlots.pdf'))
for(p in pcaPlots){
print(p)
}
invisible(dev.off())4.4 PCA plots
The 2-dimensional principal component analysis plot shows which samples group together when plotted in a reduced dimension. The 2D PCA reduced dimension conserves as much of the variance in the dataset as is possible for any 2D embedding of the data. It provides a useful birds-eye view of the data and an intuition as to which factors may drive the differences between samples or groups.
4.4.1 AnalysisGroup
4.4.2 group
4.5 Correlation Plot
The pairwise correlation plot provides a more detailed view of which samples are more similar or different.
M <- stats::cor(vsd.counts)
corrplot::corrplot(corr = M, order = 'hclust', method = 'square', type = 'lower', tl.srt = 45, addCoef.col = 'white')
#save image to folder
pdf(file = file.path(imagesDir, 'correlationPlot.pdf'))
corrplot::corrplot(corr = M, order = 'hclust', method = 'square', type = 'lower', tl.srt = 45, addCoef.col = 'white')
invisible(dev.off())4.6 Heatmaps
4.6.1 Top 100 most highly variable genes
The heatmap below summarizes the experiment, and the apparent relationship between samples, based on the 100 highest variance genes. Each column is a sample, and each row is a gene. Both rows and columns are clustered using euclidean distance and complete linkage.
select <- na.omit(names(sort(apply(X = vsd.counts, MARGIN = 1, FUN = var),decreasing = T))[1:100])
df <- as.data.frame(colData[,c("group","AnalysisGroup")])
pheatmap::pheatmap(vsd.counts[select,],
cluster_rows=TRUE,
show_rownames=FALSE,
cluster_cols=TRUE,
annotation_col=df,
main = 'Heatmap of the Normalized Expression Values (VST) of \n Top 100 Genes with highest variance across samples')
#save image to folder
pdf(file = file.path(imagesDir, 'heatmap.pdf'))
pheatmap::pheatmap(vsd.counts[select,],
cluster_rows=TRUE,
show_rownames=FALSE,
cluster_cols=TRUE,
annotation_col=df,
main = 'Heatmap of the Normalized Expression Values (VST) of \n Top 100 Genes with highest variance across samples')
invisible(dev.off())5 Exploratory Plots and Tables
5.1 Summary of up/down regulated genes - volcano plot
This volcano plot summarizes the differential expression landscape in the comparison between the two groups.
p <- ggplot(DE, aes(x = log2FoldChange, y = -log10(pvalue))) + geom_point(aes(color = padj < 0.1))
print(p)
#save image to folder
pdf(file = file.path(imagesDir, 'volcanoPlot.pdf'))
print(p)
invisible(dev.off())5.2 Summary of up/down regulated genes - bar plots
These bar plots summarizes the number of significantly upregulated/downregulated number of genes based on different adjusted p-value (selected adjusted p-values are 0.001, 0.01, 0.05, and 0.1 - see facet headers) and log2 fold change thresholds (on the x-axis) used to define the significance levels.
filterUP <- function(df, log2fc = 1, p = 0.1) {nrow(df[df$log2FoldChange >= log2fc & !is.na(df$padj) & df$padj <= p,])}
filterDOWN <- function(df, log2fc = 1, p = 0.1) {nrow(df[df$log2FoldChange < -log2fc & !is.na(df$padj) & df$padj <= p,])}
pVals <- c(0.001, 0.01, 0.05, 0.1)
fcVals <- c(0:(max(DE$log2FoldChange)+1))
summary <- do.call(rbind, lapply(pVals, function(p) {
do.call(rbind, lapply(fcVals, function(f){
up <- filterUP(DE, f, p)
down <- filterDOWN(DE, f, p)
return(data.frame("log2FoldChange" = f, "padj" = p,
"upRegulated" = up, "downRegulated" = down))
}))
}))
mdata <- melt(summary, id.vars = c('log2FoldChange', 'padj'))
p <- ggplot(mdata, aes(x = log2FoldChange, y = value)) + geom_bar(aes(fill = variable), stat = 'identity', position = 'dodge') + facet_grid(~ padj) + theme(legend.position = 'bottom', legend.title = element_blank()) + labs(title = 'Number of differentially up/down regulated genes', subtitle = 'based on different p-value and log2foldChange cut-off values')
print(p)
#save image to folder
pdf(file = file.path(imagesDir, 'up_down_regulated_genes_summary.pdf'))
print(p)
invisible(dev.off())5.3 Interactive box plots of genes with significant differential expression
This interactive plot lets you see genes’ position in the volcano plot, as well as their expression levels in the cases and the controls (in the box plot on the left side). Use the search box to find genes of interest. Notice that only top 1000 genes that have an adjusted p-value less than 0.1 and absolute log2 fold change value of greater than 1 are plotted.
select <- rownames(DEsubset)
if(length(select) > 1) {
expressionLevels <- reshape2::melt(vsd.counts[select,])
colnames(expressionLevels) <- c('geneId', 'sampleName', 'expressionLevel')
expressionLevels$group <- colData[expressionLevels$sampleName,]$group
expressionLevels$AnalysisGroup <- colData[expressionLevels$sampleName,]$AnalysisGroup
matchIds <- match(expressionLevels$geneId, rownames(DE))
expressionLevels$padj <- DE[matchIds,]$padj
expressionLevels$log2FoldChange <- DE[matchIds,]$log2FoldChange
sd <- SharedData$new(expressionLevels, ~geneId)
lineplot <- plot_ly(sd, x = ~sampleName, y = ~expressionLevel) %>%
group_by(geneId) %>%
add_lines(text = ~geneId, hoverinfo = "text", color = ~AnalysisGroup)
volcanoPlot <- plot_ly(sd, x = ~log2FoldChange, y = ~-log10(padj)) %>%
add_markers(text = ~geneId, hoverinfo = "text")
subplot(
plot_ly(sd, y = ~expressionLevel, color = ~AnalysisGroup) %>%
add_boxplot(),
volcanoPlot
) %>% highlight(on = 'plotly_click', off = 'plotly_doubleclick', selectize = TRUE)
} else {
cat("Couldn't detect at least two genes satisfying the p-value and fold change thresholds\n")
}if(runGO == FALSE) {
cat("Warning:Skipping GO analysis because `organism` option is not set in settings.yaml file\n")
} else if (curl::has_internet() == FALSE){
#gProfiler tool needs internet access to work. So, go analysis module is conditional
runGO <- FALSE
cat("Warning:Skipping GO analysis as there is no internet connection to query https://biit.cs.ut.ee/gprofiler/\n")
}cat("# GO Term Enrichment Analysis\n")6 GO Term Enrichment Analysis
cat("\n The following table list GO terms for differentially expressed genes. GO term analysis was carried out using [g:Profiler](https://cran.r-project.org/web/packages/gProfileR/) R package. Differentially expressed genes are defined as those genes with adjusted p-value of less than 0.1.")The following table list GO terms for differentially expressed genes. GO term analysis was carried out using g:Profiler R package. Differentially expressed genes are defined as those genes with adjusted p-value of less than 0.1.
#filter genes for differential expression
DEgenes <- rownames(DE[!is.na(DE$padj) & DE$padj < 0.1,])
#search for enriched GO terms within the GO and pathway domains
goResults <- gProfileR::gprofiler(query = DEgenes,
organism = organism,
hier_filtering = "none",
significant = TRUE,
src_filter = c('GO', 'KEGG', 'REAC', 'CORUM'))
#order by p-value
goResults <- goResults[order(goResults$p.value),]
#save full GO term table to disk
goResultsFile <- file.path(workdir, paste0(prefix, '.GOterms.tsv'))
write.table(x = goResults,
file = goResultsFile,
quote = FALSE, sep = '\t')
#only display top GO terms in the HTML report.
max <- ifelse(nrow(goResults) > 1000, 1000, nrow(goResults))
DT::datatable(goResults[1:max,],
extensions = c('Buttons', 'FixedColumns', 'Scroller'),
options = list(fixedColumns = TRUE,
scrollY = 400,
scrollX = TRUE,
scroller = TRUE,
dom = 'Bfrtip',
buttons = c('colvis', 'copy', 'print', 'csv','excel', 'pdf'),
columnDefs = list(
list(targets = c(0,1,2,5,9,11,13,14), visible = FALSE)
)
),
filter = 'bottom'
)