Commit d219419f authored by Erwan DELAGE's avatar Erwan DELAGE
Browse files

Add plots to visualize DESeq2 significant ASVs along with their abundance and taxonomy

parent 0e2c7b24
......@@ -12,4 +12,5 @@ dependencies:
- r-svglite
- r-ggrepel
- r-ggiraph
- r-rcolorbrewer
\ No newline at end of file
- r-rcolorbrewer
- r-gridextra
......@@ -25,6 +25,8 @@ library(hrbrthemes, quietly = TRUE)
library("ggrepel", quietly=TRUE)
library("ggiraph", quietly=TRUE)
library("RColorBrewer", quietly=TRUE)
library(gridExtra)
library(grid)
#######################
# FILTERING FUNCTIONS #
......@@ -209,6 +211,139 @@ volcano <- function(df, cond1, cond2, diffThreshold=0, adjPValThreshold=0.05, sh
return(p)
}
getLastTaxo <- function(deseq2_res){
#' Get lower taxonomic assignatino (to plot on the graph)
#'
#' @description Get lower taxonomic assignatino (to plot on the graph)
#'
#' @param deseq2_res DataFrame of deseq2 results
#' @return vector of annotation
lastTaxo = c()
taxoRanks = c("Species","Genus")
for (i in 1:dim(deseq2_res[1])){
taxoVal = ""
for (taxoRank in taxoRanks){
if (!grepl("unidentified", deseq2_res[i,taxoRank], fixed=TRUE) && !grepl("uncultured", deseq2_res[i,taxoRank], fixed=TRUE) && deseq2_res[i,taxoRank] != ""){
taxoVal = deseq2_res[i,taxoRank]
break
}
}
lastTaxo <- append(lastTaxo, taxoVal)
}
return(lastTaxo)
}
heatmapColorScale <- function(v){
#' Define a discrete color scale and categorize ASVs according to abundance
#'
#' @description Define a discrete color scale and categorize ASVs according to abundance. Temporary solution, waiting for
#' the proper way to define a multi-continuous scale
#'
#' @param v vector of TSS abundance
#' @return factor of level of abundance
y = c()
for(i in 1:length(v)){
if (v[i] == 0){
y <- append(y, "0")
}else if (v[i] > 20){
y <- append(y, ">20")
}else if (v[i] > 10) {
y <- append(y, "10-20")
}else if (v[i] > 5) {
y <- append(y, "5-10")
}else if (v[i] > 1) {
y <- append(y, "1-5")
}else if (v[i] > 0.1) {
y <- append(y, "1e-1")
}else if (v[i] > 0.01) {
y <- append(y, "1e-2")
}else if (v[i] > 0.001) {
y <- append(y, "1e-3")
}else if (v[i] > 0.0001) {
y <- append(y, "1e-4")
}else if (v[i] > 0.00001) {
y <- append(y, "1e-5")
}else {
y <- append(y, "<1e-5")
}
}
return(fct_rev(factor(y,levels=c("0","<1e-5", "1e-5", "1e-4", "1e-3", "1e-2", "1e-1", "1-5", "5-10","10-20", ">20"))))
}
plot_deseq2 <- function(deseq2_res, tss_abundance, context, context_condition, cond1, cond2) {
#' Build a plot showing differential ASVs, their fold change, their taxonomic assignation and their ordre of abundance
#'
#'
#' @param deseq2_res DataFrame of deseq2 results
#' @param tss_abundance DataFrame of TSS normalized abundance
#' @param context DataFrame of contextual metadata
#' @param context_condition Contextual variable deseq2 was tested on
#' @param cond1 First condition tested
#' @param cond2 Second condition tested
#' @return ggplot
# Add mean abundance according to condition
for (feature_id in rownames(deseq2_res)) {
deseq2_res[feature_id,cond1] <- cond1_mean <- rowMeans(tss_abundance[feature_id,context[context_condition] == cond1])
deseq2_res[feature_id,cond2] <- cond2_mean <- rowMeans(tss_abundance[feature_id,context[context_condition] == cond2])
}
# Format data frame
deseq2_res <- deseq2_res %>%
arrange(log2FoldChange) %>% # Sort rows according to log fold change
mutate(id=as_factor(rownames(deseq2_res))) %>% # Add an "id" column from rownames, with factors
mutate(pvalPlot=-log10(padj)) %>% # Transform adjusted p-value to reflect it as the node size
mutate(posLabel = case_when(
log2FoldChange < 0 ~ -0.2,
log2FoldChange > 0 ~ 1.4
)) # Create a column for node label position, right side of the node when logFoldChange is negative, left side otherwise
# If info is available at the genus or specie level, we add it on the plot
deseq2_res <- deseq2_res %>%
mutate(lastTaxo = getLastTaxo(deseq2_res)) # Get last taxonomic assignation
# Fold change plot
foldChangePlot <- deseq2_res %>%
ggplot(aes(log2FoldChange, id, colour = Family, size=pvalPlot)) +
geom_point() +
geom_vline(xintercept=0, linetype="dotted", size=0.5) + # Add vertical bar at null fold change
scale_y_discrete(position = "right") + # Put features labels on the right side of the plot
geom_text(aes(label=lastTaxo),hjust=deseq2_res$posLabel, vjust=0.3, size=4) +
theme(
legend.position = c(.15, .75) # Set legend to the uper left part of the figure
) +
labs(size="-log10 padj", y="Feature ID") # Rename legend title for node size, and y axis
# Wide to long dataframe
abundance_df <- deseq2_res %>%
select(eval(cond1),eval(cond2),id) %>% # Select only interesting values
gather(!!context_condition,"TSS", -id) %>% # Turn from wide to long
mutate(Proportion=heatmapColorScale(TSS)) # Assign TSS values to a discrete category
# Add factors to order heatmap
abundance_df <- abundance_df %>%
mutate(!!context_condition := factor(abundance_df[[context_condition]], levels = c(cond2,cond1)))
# Heatmap
heatmap <- abundance_df %>% ggplot(aes_string(context_condition, "id", fill="Proportion")) +
geom_tile() +
scale_fill_viridis_d() +
theme(axis.title.y=element_blank(),
axis.text.y=element_blank())
# Combine the two plots into one
p <- grid.arrange(foldChangePlot,heatmap,ncol=2, widths=c(10,2), top = textGrob(paste("DESeq2 :",cond2,"vs",cond1) ,gp=gpar(fontsize=20,font=3)))
return(p)
}
boxplot_significant_features <- function(feature_id, abundance_table, context_table, context_condition, taxon, outfile="") {
......@@ -392,6 +527,12 @@ for(i in 1:(length(category)-1)){
h <- girafe(ggobj = p, pointsize = 16, width_svg = 12, height_svg = 10)
htmlwidgets::saveWidget(h, paste0(outdir,"/volcano_plot.html"), selfcontained = TRUE)
# Log fold change + heatmap plot
if (taxoProvided){
p <- plot_deseq2(resSig, tss_abundance, context, deseq2_params[["condition"]], category[i], category[j])
ggsave(file=paste0(outdir,"/deseq2_results.svg"), plot = p, width = 20, height = 13)
}
# Put boxplot in subdirectory
outdir = paste0(outdir, "/boxplot_tss_abundance")
dir.create(outdir, showWarnings = FALSE)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment