Browse code

fixed RMD check for branch

Ludt authored on 18/02/2021 09:49:03
Showing 5 changed files

... ...
@@ -23,6 +23,7 @@ Imports:
23 23
     AnnotationDbi,
24 24
     bs4Dash,
25 25
     colorspace,
26
+	colourpicker,
26 27
     ComplexHeatmap,
27 28
     dendextend,
28 29
     DESeq2,
... ...
@@ -18,6 +18,7 @@ export(geneinfo_2_html)
18 18
 export(get_aggrscores)
19 19
 export(get_expression_values)
20 20
 export(ggs_graph)
21
+export(ggs_volcano)
21 22
 export(go_2_html)
22 23
 export(gs_alluvial)
23 24
 export(gs_dendro)
... ...
@@ -92,6 +93,7 @@ importFrom(bs4Dash,bs4ValueBoxOutput)
92 93
 importFrom(bs4Dash,renderbs4InfoBox)
93 94
 importFrom(bs4Dash,renderbs4ValueBox)
94 95
 importFrom(colorspace,rainbow_hcl)
96
+importFrom(colourpicker,colourInput)
95 97
 importFrom(dendextend,branches_attr_by_clusters)
96 98
 importFrom(dendextend,set)
97 99
 importFrom(dplyr,"%>%")
... ...
@@ -13,6 +13,7 @@
13 13
 #' bs4ValueBoxOutput renderbs4InfoBox renderbs4ValueBox
14 14
 #' bs4TabPanel bs4TabCard
15 15
 #' @importFrom colorspace rainbow_hcl
16
+#' @importFrom colourpicker colourInput
16 17
 #' @importFrom ComplexHeatmap Heatmap HeatmapAnnotation draw
17 18
 #' @importFrom dendextend branches_attr_by_clusters set
18 19
 #' @importFrom DESeq2 vst counts estimateSizeFactors normalizationFactors sizeFactors
... ...
@@ -17,7 +17,7 @@
17 17
 #' @param FDR Numeric value, specifying the significance level for thresholding
18 18
 #' adjusted p-values. Defaults to 0.05.
19 19
 #' @param color Character string to specify color of filtered points in the plot.
20
-#'  Default color is red.
20
+#'  Default color is #0092AC.
21 21
 #' @param volcano_labels Integer, maximum number of labels for the gene sets to be
22 22
 #' plotted as labels on the volcano scatter plot.
23 23
 #' @param plot_title Character string, to specify the title of the plot,
... ...
@@ -39,7 +39,6 @@
39 39
 #' rownames(dds_macrophage) <- substr(rownames(dds_macrophage), 1, 15)
40 40
 #' dds_macrophage <- estimateSizeFactors(dds_macrophage)
41 41
 #'
42
-#' vst_macrophage <- vst(dds_macrophage)
43 42
 #'
44 43
 #' # annotation object
45 44
 #' anno_df <- data.frame(
... ...
@@ -61,8 +60,7 @@
61 60
 #' res_enrich <- shake_topGOtableResult(topgoDE_macrophage_IFNg_vs_naive)
62 61
 #' res_enrich <- get_aggrscores(res_enrich, res_de, anno_df)
63 62
 #'
64
-#' ggs_volcano(vst_macrophage,
65
-#'            res_de,
63
+#' ggs_volcano(res_de,
66 64
 #'            res_enrich,
67 65
 #'            anno_df,
68 66
 #'            geneset_id = res_enrich$gs_id[1]
... ...
@@ -76,7 +74,7 @@ ggs_volcano <- function(res_de,
76 74
                        genelist = NULL,
77 75
                        FDR = 0.05,
78 76
                        color = "red",
79
-                       volcano_labels = 10,
77
+                       volcano_labels = 25,
80 78
                        plot_title = NULL
81 79
                        ) {
82 80
 
... ...
@@ -96,27 +94,48 @@ ggs_volcano <- function(res_de,
96 94
       thisset_members_ids <- annotation_obj$gene_id[match(thisset_members, annotation_obj$gene_name)]
97 95
     }
98 96
   }else {
99
-    # overridable via a list
100
-    if (!all(genelist %in% rownames(se))) {
101
-      not_there <- genelist[!(genelist %in% rownames(se))]
97
+    # overwritable via a list
98
+    if (!all(genelist %in% res_de@rownames)) {
99
+      not_there <- genelist[!(genelist %in% res_de@rownames)]
102 100
       warning("Some of the provided gene ids were not found in the SummarizedExperiment",
103 101
               "\nNot found: ", not_there)
104 102
     }
105
-    thisset_members_ids <- intersect(rownames(se), genelist)
103
+    thisset_members_ids <- intersect(res_de@rownames, genelist)
106 104
     thisset_name <- "Custom list"
107 105
   }
108 106
 
109 107
   thisset_members_ids
110 108
   thisset_members
111 109
   
110
+  complete_genes_ids <- res_de@rownames
111
+  complete_genes <- annotation_obj$gene_name[match(complete_genes_ids, annotation_obj$gene_id)]
112
+  padj_complete <- res_de[complete_genes_ids, "padj"]
113
+  filter_info_complete <- sapply(padj_complete, function(x) x <= FDR)
114
+  padj_complete <- sapply(padj_complete, function(x) -log10(x))
115
+  log2FoldChange_complete <- res_de[complete_genes_ids, "log2FoldChange"]
116
+  gene_set_belong <- complete_genes_ids %in% thisset_members_ids
117
+  
118
+  filter_info_complete <- filter_info_complete & gene_set_belong
119
+
120
+  
121
+  thisset_complete_data <- data.frame(complete_genes_ids, padj_complete, log2FoldChange_complete, filter_info_complete, gene_set_belong)
122
+  colnames(thisset_complete_data) <- c("genes", "logTransformedpvalue", "log2FoldChange", "significant", "belonging")
123
+  
124
+  
125
+  # Prepare plotting
126
+  volcano_df_complete <- thisset_complete_data
127
+  volcano_df_complete$genes_name <- complete_genes
128
+  max_x <- max(abs(range(thisset_complete_data["log2FoldChange"])))
129
+  limit_x <- max_x * c(-1, 1)
130
+  
112 131
   # retrieve adjusted pvalues and log2Fold Change for genes
113
-  padj <- res_de[thisset_members_ids, "padj"]
114
-  padjlog <- sapply(padj, function(x) -log10(x))
115
-  log2FoldChange <- res_de[thisset_members_ids, "log2FoldChange"]
116
-  filter_info <- sapply(padj, function(x) x <= FDR)
132
+  #padj <- res_de[thisset_members_ids, "padj"]
133
+  #padjlog <- sapply(padj, function(x) -log10(x))
134
+  #log2FoldChange <- res_de[thisset_members_ids, "log2FoldChange"]
135
+  #filter_info <- sapply(padj, function(x) x <= FDR)
117 136
   
118 137
   # Prepare data to plot
119
-  thisset_data <- data.frame(thisset_members_ids, padj, log2FoldChange, padjlog, filter_info)
138
+  #thisset_data <- data.frame(thisset_members_ids, padj, log2FoldChange, padjlog, filter_info)
120 139
 
121 140
   
122 141
   if (is.null(plot_title)) {
... ...
@@ -125,30 +144,55 @@ ggs_volcano <- function(res_de,
125 144
     title <- plot_title
126 145
   }
127 146
   
128
-  colnames(thisset_data) <- c("genes", "pvalue", "log2FoldChange", "logTransformedpvalue", "significant")
147
+ # colnames(thisset_data) <- c("genes", "pvalue", "log2FoldChange", "logTransformedpvalue", "significant")
129 148
   
130 149
   
131 150
   # Prepare plotting
132
-  volcano_df <- thisset_data
133
-  volcano_df$genes_name <- thisset_members
134
-  max_x <- max(abs(range(thisset_data["log2FoldChange"])))
135
-  limit_x <- max_x * c(-1, 1)
151
+ # volcano_df <- thisset_data
152
+  #volcano_df$genes_name <- thisset_members
153
+  #max_x <- max(abs(range(thisset_data["log2FoldChange"])))
154
+  #limit_x <- max_x * c(-1, 1)
136 155
 
137 156
 
138 157
   # Plot data
139
-   p <- ggplot(volcano_df, aes(x = log2FoldChange, y = logTransformedpvalue, label = genes_name)) +
140
-    geom_point(aes(color = significant)) +
158
+  #p <- ggplot(volcano_df, aes(x = log2FoldChange, y = logTransformedpvalue, label = genes_name)) +
159
+  #geom_point(aes(color = significant)) +
160
+  #  labs(x = "log2Fold Change",
161
+  #       y = "-log10 p-value", color = "p-value <= FDR") + 
162
+  #  scale_x_continuous(limits = limit_x) +
163
+  #  scale_color_manual(labels = c("significant", "not significant"), 
164
+  #                     breaks = c("TRUE", "FALSE"), values = c(color, "grey25")) +
165
+  #  theme(legend.title = element_text(size = 11, face = "bold"), legend.text = element_text(size = 10)) +
166
+  #  geom_text_repel()
167
+  
168
+  p <- ggplot(volcano_df_complete, aes_string(x = "log2FoldChange", y = "logTransformedpvalue")) +
169
+    geom_point(aes_string(color = "significant", alpha = "belonging")) +
141 170
     labs(x = "log2Fold Change",
142
-         y = "-log10 p-value", color = "p-value <= FDR") + 
171
+         y = "-log10 p-value",
172
+         color = "p-value <= FDR") + 
143 173
     scale_x_continuous(limits = limit_x) +
144 174
     scale_color_manual(labels = c("significant", "not significant"), 
145
-                       breaks = c("TRUE", "FALSE"), values = c(color, "grey25")) +
146
-    theme(legend.title = element_text(size = 11, face = "bold"), legend.text = element_text(size = 10)) +
147
-    geom_text_repel()
175
+                       breaks = c("TRUE", "FALSE"),
176
+                       values = c(color, "grey25")) +
177
+    scale_alpha_manual(breaks = c("TRUE", "FALSE"),
178
+                       values = c(1, 1/10)) +
179
+    theme_bw() +
180
+    theme(legend.title = element_text(size = 11, face = "bold"),
181
+          legend.text = element_text(size = 10))  
182
+    
148 183
 
149 184
     
185
+  # adding labels to the significant points of the geneset
186
+  p <- p + geom_text_repel(
187
+    data = subset(volcano_df_complete, filter_info_complete),
188
+    aes_string(label = "genes_name"),
189
+    size = 4,
190
+    max.overlaps = volcano_labels)
191
+
192
+
150 193
   
151 194
   # handling the title
152 195
   p <- p + ggtitle(title)
196
+  p <- p + guides(alpha = FALSE)
153 197
   return(p)
154 198
 }
155 199
new file mode 100644
... ...
@@ -0,0 +1,97 @@
1
+% Generated by roxygen2: do not edit by hand
2
+% Please edit documentation in R/ggs_volcano.R
3
+\name{ggs_volcano}
4
+\alias{ggs_volcano}
5
+\title{Plot a volcano plot of the gene signature on the data}
6
+\usage{
7
+ggs_volcano(
8
+  res_de,
9
+  res_enrich,
10
+  annotation_obj = NULL,
11
+  gtl = NULL,
12
+  geneset_id = NULL,
13
+  genelist = NULL,
14
+  FDR = 0.05,
15
+  color = "red",
16
+  volcano_labels = 25,
17
+  plot_title = NULL
18
+)
19
+}
20
+\arguments{
21
+\item{res_de}{A \code{DESeqResults} object.}
22
+
23
+\item{res_enrich}{A \code{data.frame} object, storing the result of the functional
24
+enrichment analysis. See more in the main function, \code{\link[=GeneTonic]{GeneTonic()}}, to check the
25
+formatting requirements (a minimal set of columns should be present).}
26
+
27
+\item{annotation_obj}{A \code{data.frame} object with the feature annotation
28
+information, with at least two columns, \code{gene_id} and \code{gene_name}.}
29
+
30
+\item{gtl}{A \code{GeneTonic}-list object, containing in its slots the arguments
31
+specified above: \code{dds}, \code{res_de}, \code{res_enrich}, and \code{annotation_obj} - the names
32
+of the list \emph{must} be specified following the content they are expecting}
33
+
34
+\item{geneset_id}{Character specifying the gene set identifier to be plotted}
35
+
36
+\item{genelist}{A vector of character strings, specifying the identifiers
37
+contained in the row names of the \code{se} input object.}
38
+
39
+\item{FDR}{Numeric value, specifying the significance level for thresholding
40
+adjusted p-values. Defaults to 0.05.}
41
+
42
+\item{color}{Character string to specify color of filtered points in the plot.
43
+Default color is #0092AC.}
44
+
45
+\item{volcano_labels}{Integer, maximum number of labels for the gene sets to be
46
+plotted as labels on the volcano scatter plot.}
47
+
48
+\item{plot_title}{Character string, to specify the title of the plot,
49
+displayed over the heatmap. If left to \code{NULL} as by default, it tries to use
50
+the information on the geneset identifier provided}
51
+}
52
+\value{
53
+A plot returned by the \code{\link[=ggplot]{ggplot()}} function
54
+}
55
+\description{
56
+Plot a volcano plot for the selected gene signature on the provided data, with the possibility to compactly display also DE only genes
57
+}
58
+\examples{
59
+library("macrophage")
60
+library("DESeq2")
61
+library("org.Hs.eg.db")
62
+library("AnnotationDbi")
63
+
64
+# dds object
65
+data("gse", package = "macrophage")
66
+dds_macrophage <- DESeqDataSet(gse, design = ~line + condition)
67
+rownames(dds_macrophage) <- substr(rownames(dds_macrophage), 1, 15)
68
+dds_macrophage <- estimateSizeFactors(dds_macrophage)
69
+
70
+
71
+# annotation object
72
+anno_df <- data.frame(
73
+  gene_id = rownames(dds_macrophage),
74
+  gene_name = mapIds(org.Hs.eg.db,
75
+                     keys = rownames(dds_macrophage),
76
+                     column = "SYMBOL",
77
+                     keytype = "ENSEMBL"),
78
+  stringsAsFactors = FALSE,
79
+  row.names = rownames(dds_macrophage)
80
+)
81
+
82
+# res object
83
+data(res_de_macrophage, package = "GeneTonic")
84
+res_de <- res_macrophage_IFNg_vs_naive
85
+
86
+# res_enrich object
87
+data(res_enrich_macrophage, package = "GeneTonic")
88
+res_enrich <- shake_topGOtableResult(topgoDE_macrophage_IFNg_vs_naive)
89
+res_enrich <- get_aggrscores(res_enrich, res_de, anno_df)
90
+
91
+ggs_volcano(res_de,
92
+           res_enrich,
93
+           anno_df,
94
+           geneset_id = res_enrich$gs_id[1]
95
+           )
96
+
97
+}