|
0
|
1 #!/usr/bin/env Rscript
|
|
|
2 Sys.setlocale("LC_CTYPE", "en_US.UTF-8") # this is necessary for handling unicode characters (data.tree package)
|
|
|
3 suppressPackageStartupMessages(library(data.tree))
|
|
|
4 suppressPackageStartupMessages(library(stringr))
|
|
|
5 suppressPackageStartupMessages(library(R2HTML))
|
|
|
6 suppressPackageStartupMessages(library(hwriter))
|
|
|
7 suppressPackageStartupMessages(library(DT))
|
|
|
8 suppressPackageStartupMessages(library(tools))
|
|
|
9 suppressPackageStartupMessages(library(scales))
|
|
|
10 suppressPackageStartupMessages(library(igraph))
|
|
|
11 suppressPackageStartupMessages(library(plotrix))
|
|
|
12 suppressPackageStartupMessages(library(png))
|
|
|
13
|
|
|
14 source("htmlheader.R")
|
|
|
15 source("config.R") # load tandem ranks info
|
|
|
16 source("utils.R")
|
|
|
17 DT_OPTIONS = options = list(pageLength = 1000, lengthMenu = c(10,50,100,1000,5000,10000))
|
|
|
18 WD = getwd() # to get script directory when run from Rserve
|
|
|
19 HTMLHEADER = htmlheader ## header (character) loaded from htmlheader.R
|
|
|
20 htmlheader = gsub("Superclusters summary","TAREAN summary", htmlheader)
|
|
|
21
|
|
|
22 evaluate_LTR_detection = function(f){
|
|
|
23 NO_LTR=NULL
|
|
|
24 if (length(readLines(f)) == 11 ){
|
|
|
25 return(NO_LTR)
|
|
|
26 }
|
|
|
27 df=read.table(f,as.is=TRUE,sep="\t", skip=11,fill=TRUE)
|
|
|
28 if (ncol(df) != 23){
|
|
|
29 #df is smaller if no pbs is detected!
|
|
|
30 return(NO_LTR)
|
|
|
31 }
|
|
|
32 df=df[!df$V13 == "",,drop=FALSE]
|
|
|
33 if (nrow(df)==0){
|
|
|
34 return(NO_LTR)
|
|
|
35 }
|
|
|
36 # criteria:
|
|
|
37 df_part=df[df$V15 >=12 & df$V20 == 23 & df$V21<df$V20,,drop=FALSE]
|
|
|
38 if (nrow(df_part) == 0){
|
|
|
39 return(NO_LTR)
|
|
|
40 }
|
|
|
41 PBS_type = gsub("_","",(str_extract_all(df_part$V13, pattern="_([A-Z][a-z]{2})", simplify=TRUE))) %>%
|
|
|
42 paste(collapse=" ")
|
|
|
43 return(PBS_type)
|
|
|
44 }
|
|
|
45
|
|
|
46
|
|
|
47
|
|
|
48 ## annotate superclusters
|
|
|
49 select_reads_id = function(index, search_type = c("cluster","supercluster")) {
|
|
|
50 ## select read if base on the supecluster index need database connection
|
|
|
51 ## HITSORTDB!
|
|
|
52 search_type = match.arg(search_type)
|
|
|
53 x = dbGetQuery(HITSORTDB,
|
|
|
54 paste0("SELECT vertexname FROM vertices WHERE vertexindex IN ",
|
|
|
55 "(SELECT vertexindex FROM communities ",
|
|
|
56 "WHERE ", search_type,"=\"", index,
|
|
|
57 "\")"))
|
|
|
58 return(x$vertexname)
|
|
|
59 }
|
|
|
60
|
|
|
61
|
|
|
62 get_reads_annotation = function(reads_id) {
|
|
|
63 ## select annotation from tables in SEQDB which has name in format *_database
|
|
|
64 annot_types = grep("_database", dbListTables(SEQDB), value = TRUE)
|
|
|
65 annot = list()
|
|
|
66 for (i in annot_types) {
|
|
|
67 query = paste0("SELECT * FROM ", i, " WHERE name IN (", paste0("\"", reads_id,
|
|
|
68 "\"", collapse = ", "), ")")
|
|
|
69 annot[[i]] = dbGetQuery(SEQDB, query)
|
|
|
70 }
|
|
|
71 return(annot)
|
|
|
72 }
|
|
|
73
|
|
|
74 supercluster_size = function(supercluster) {
|
|
|
75 x = dbGetQuery(HITSORTDB, paste0("SELECT count(*) FROM vertices WHERE vertexindex IN ",
|
|
|
76 "(SELECT vertexindex FROM communities ", "WHERE supercluster=\"", supercluster,
|
|
|
77 "\")"))
|
|
|
78 return(x$"count(*)")
|
|
|
79 }
|
|
|
80
|
|
|
81
|
|
|
82 cluster_annotation = function(cluster, search_type = c("cluster", "supercluster")){
|
|
|
83 ## searcheither for cluster or supercluster annotation
|
|
|
84 ## read annotation from sqlite databases database is access though SEQDB (sequence
|
|
|
85 ## annotation) and HITSORTDB - clustering information
|
|
|
86 search_type = match.arg(search_type)
|
|
|
87 reads_id = select_reads_id(cluster, search_type)
|
|
|
88 annot = get_reads_annotation(reads_id)
|
|
|
89 return(annot)
|
|
|
90 }
|
|
|
91
|
|
|
92 get_tarean_info = function(cluster, search_type = c("cluster", "supercluster")){
|
|
|
93 search_type = match.arg(search_type)
|
|
|
94 if (search_type == "cluster") {
|
|
|
95 search_type = "[index]"
|
|
|
96 }
|
|
|
97 tarean_info = dbGetQuery(HITSORTDB,
|
|
|
98 paste0(
|
|
|
99 "SELECT [index], supercluster, satellite_probability, tandem_rank, size_real, satellite FROM cluster_info WHERE ",
|
|
|
100 search_type, " = ", cluster))
|
|
|
101 nhits = sum(tarean_info$size_real[tarean_info$tandem_rank %in% 1:2])
|
|
|
102 proportion = nhits/sum(tarean_info$size_real)
|
|
|
103 return( list(
|
|
|
104 nhits = nhits,
|
|
|
105 proportion = proportion,
|
|
|
106 tarean_info = tarean_info)
|
|
|
107 )
|
|
|
108
|
|
|
109 }
|
|
|
110
|
|
|
111 get_ltr_info = function(supercluster){
|
|
|
112 ltr_info = dbGetQuery(HITSORTDB,
|
|
|
113 paste0(
|
|
|
114 "SELECT [index], supercluster, ltr_detection, size_real FROM cluster_info WHERE ",
|
|
|
115 "supercluster", " = ", supercluster))
|
|
|
116 return(ltr_info)
|
|
|
117 }
|
|
|
118
|
|
|
119 summarize_annotation = function(annot, size = NULL){
|
|
|
120 db_id = do.call(rbind, annot)$db_id
|
|
|
121 cl_string = gsub("^.+/","",gsub(":.+$", "", gsub("^.+#", "", db_id)))
|
|
|
122 weight = do.call(rbind, annot)$bitscore
|
|
|
123 domain = str_replace(str_extract(str_replace(db_id, "^.+#", ""), ":.+"), ":",
|
|
|
124 "")
|
|
|
125 domain[is.na(domain)] = ""
|
|
|
126 domain_table = table(cl_string, domain)
|
|
|
127 dt = as.data.frame(domain_table)
|
|
|
128 ## remove no count
|
|
|
129 dt = dt[dt$Freq != 0, ,drop = FALSE]
|
|
|
130 rownames(dt) = paste(dt$cl_string,dt$domain)
|
|
|
131 if (!is.null(size)){
|
|
|
132 dt$proportion = dt$Freq / size
|
|
|
133
|
|
|
134 }
|
|
|
135 return(dt)
|
|
|
136 }
|
|
|
137
|
|
|
138 annot2colors = function(annot, ids, base_color = "#00000050"){
|
|
|
139 annot_table = do.call(rbind, annot)
|
|
|
140 domains = str_replace(str_extract(str_replace(annot_table$db_id, "^.+#", ""), ":.+"), ":","")
|
|
|
141 other_annot = str_replace(annot_table$db_id, "^.+#", "")
|
|
|
142 complete_annot = ifelse(is.na(domains), other_annot, domains)
|
|
|
143 names(complete_annot) = annot_table$name
|
|
|
144 unique_complete_annot = names (table(complete_annot) %>% sort(decreasing = TRUE))
|
|
|
145 color_key = unique_colors[seq_along(unique_complete_annot)]
|
|
|
146 names(color_key) = unique_complete_annot
|
|
|
147 color_table = rep(base_color, length(ids))
|
|
|
148 names(color_table) = ids
|
|
|
149 color_table[names(complete_annot)] = color_key[complete_annot]
|
|
|
150 return(list( color_table = color_table, legend = color_key))
|
|
|
151 }
|
|
|
152
|
|
|
153
|
|
|
154 read_annotation_to_tree = function(supercluster, TREE_TEMPLATE = CLS_TREE) {
|
|
|
155 ## CLS_TREE is datatree containing 'taxonomy' information of repeats, attribute of
|
|
|
156 ## each node/ leave if filled up from annotation annotation is added to leaves
|
|
|
157 ## only!! Inner node are NA or NULL
|
|
|
158 annot0 = cluster_annotation(supercluster, search_type = "supercluster")
|
|
|
159 ## keep only built in annotation
|
|
|
160 annot = list(protein_databse = annot0$protein_database, dna_database = annot0$dna_database)
|
|
|
161 annot$dna_database$bitscore = annot$dna_database$bitscore / 2 #DNA bitscore corection, blastx - more weight
|
|
|
162 annot = do.call(rbind, annot)
|
|
|
163 ## remove duplicated hits, keep best
|
|
|
164 annot_clean = annot[order(annot$bitscore, decreasing = TRUE), ]
|
|
|
165 annot_clean = annot_clean[!duplicated(annot_clean$name), ]
|
|
|
166 tarean_info = get_tarean_info(supercluster, search_type = "supercluster")
|
|
|
167 ltr_info = get_ltr_info(supercluster)
|
|
|
168 db_id = annot_clean$db_id
|
|
|
169 cl_string = gsub(":.+$", "", gsub("^.+#", "", db_id))
|
|
|
170 weight = annot_clean$bitscore
|
|
|
171 domain = str_replace(str_extract(str_replace(db_id, "^.+#", ""), ":.+"), ":",
|
|
|
172 "")
|
|
|
173 domain[is.na(domain)] = "NA"
|
|
|
174 mean_weight_table = by(weight, INDICES = list(cl_string, domain), FUN = function(x) signif(mean(x)))
|
|
|
175 mean_weight_table[is.na(mean_weight_table)] = 0
|
|
|
176 total_weight_table = by(weight, INDICES = list(cl_string, domain), FUN = sum)
|
|
|
177 total_weight_table[is.na(total_weight_table)] = 0
|
|
|
178 cl_table = table(cl_string)
|
|
|
179 domain_table = table(cl_string, domain)
|
|
|
180 cls_tree = Clone(TREE_TEMPLATE)
|
|
|
181 cls_tree$size = supercluster_size(supercluster)
|
|
|
182 cls_tree$ltr_info = ltr_info
|
|
|
183 for (i in cls_tree$leaves) {
|
|
|
184 if (i$full_name %in% names(cl_table)) {
|
|
|
185 i$nhits = cl_table[[i$full_name]]
|
|
|
186 i$domains = domain_table[i$full_name, ]
|
|
|
187 names(i$domains) = colnames(domain_table)
|
|
|
188 i$mean_weight = mean_weight_table[i$full_name, ]
|
|
|
189 i$total_weight = total_weight_table[i$full_name, ]
|
|
|
190 i$proportion = signif(i$nhits/cls_tree$size, 3)
|
|
|
191 } else {
|
|
|
192 if (i$name == "satellite"){
|
|
|
193 i$nhits = tarean_info$nhits
|
|
|
194 i$tandem_rank = tarean_info$tarean_info$tandem_rank
|
|
|
195 i$proportion = tarean_info$proportion
|
|
|
196 i$tarean_info = tarean_info$tarean_info
|
|
|
197 }else{
|
|
|
198 i$proportion = 0
|
|
|
199 }
|
|
|
200 }
|
|
|
201 }
|
|
|
202 ## create domain string for printing:
|
|
|
203 for (i in Traverse(cls_tree)) {
|
|
|
204 if (is.null(i$domains)) {
|
|
|
205 i$features = ""
|
|
|
206 } else if (sum(i$domains) == 0) {
|
|
|
207 i$features = ""
|
|
|
208 } else {
|
|
|
209 dom = i$domains[!names(i$domains) == "NA"]
|
|
|
210 i$features = gsub(" *\\(\\)", "", pasteDomains(dom))
|
|
|
211 }
|
|
|
212 }
|
|
|
213 ## add LTR info:
|
|
|
214 if (any(!cls_tree$ltr_info$ltr_detection %in% "None")){
|
|
|
215 tr = FindNode(cls_tree, "LTR")
|
|
|
216 tr$features = "LTR/PBS"
|
|
|
217 }
|
|
|
218 return(cls_tree)
|
|
|
219 }
|
|
|
220
|
|
|
221 pasteDomains = function(dom){
|
|
|
222 d = dom[dom != 0]
|
|
|
223 if (length(d) == 0){
|
|
|
224 return("")
|
|
|
225 }else{
|
|
|
226 paste0(d, " (", names(d), ")", sep = ", ", collapse="")
|
|
|
227 }
|
|
|
228
|
|
|
229 }
|
|
|
230
|
|
|
231 rescale = function(x, newrange) {
|
|
|
232 # taken from plotrix package
|
|
|
233 if (missing(x) | missing(newrange)) {
|
|
|
234 usage.string <- paste("Usage: rescale(x,newrange)\n", "\twhere x is a numeric object and newrange is the new min and max\n",
|
|
|
235 sep = "", collapse = "")
|
|
|
236 stop(usage.string)
|
|
|
237 }
|
|
|
238 if (is.numeric(x) && is.numeric(newrange)) {
|
|
|
239 xna <- is.na(x)
|
|
|
240 if (all(xna))
|
|
|
241 return(x)
|
|
|
242 if (any(xna))
|
|
|
243 xrange <- range(x[!xna]) else xrange <- range(x)
|
|
|
244 if (xrange[1] == xrange[2])
|
|
|
245 return(x)
|
|
|
246 mfac <- (newrange[2] - newrange[1])/(xrange[2] - xrange[1])
|
|
|
247 return(newrange[1] + (x - xrange[1]) * mfac)
|
|
|
248 } else {
|
|
|
249 warning("Only numeric objects can be rescaled")
|
|
|
250 return(x)
|
|
|
251 }
|
|
|
252 }
|
|
|
253
|
|
|
254 add_value_to_nodes = function(tr, value_names = c("nhits", "domains", "total_weight",
|
|
|
255 "proportion")) {
|
|
|
256 ## propagate values from leaves to nodes, assume that nodes are not defined
|
|
|
257 ## (either NA or NULL) leaves coud be be either numeric of list, if list values
|
|
|
258 ## are summed up based on the names
|
|
|
259 for (vn in value_names) {
|
|
|
260 for (i in Traverse(tr)) {
|
|
|
261 w = add_leaves_value(i$leaves, vn)
|
|
|
262 i[[vn]] = w
|
|
|
263 }
|
|
|
264 }
|
|
|
265 }
|
|
|
266
|
|
|
267 add_leaves_value = function(x, name) {
|
|
|
268 ## check if annotation is multivalue or single value, propagates values from
|
|
|
269 ## leaves to root
|
|
|
270 n = unique(unlist(lapply(lapply(x, "[[", name), names)))
|
|
|
271 if (is.null(n)) {
|
|
|
272 ## no names given, we can sum everything to one value
|
|
|
273 total = sum(unlist(sapply(x, "[[", name)), na.rm = TRUE)
|
|
|
274 return(total)
|
|
|
275 } else {
|
|
|
276 total = numeric(length(n))
|
|
|
277 names(total) = n
|
|
|
278 xvals = lapply(x, "[[", name)
|
|
|
279 N = 0
|
|
|
280 for (i in xvals) {
|
|
|
281 for (j in names(i)) {
|
|
|
282 total[j] = total[j] + i[j]
|
|
|
283 N = N + 1
|
|
|
284 }
|
|
|
285 }
|
|
|
286 return(total)
|
|
|
287 }
|
|
|
288 }
|
|
|
289
|
|
|
290
|
|
|
291 trmap = function(tr, xl = 0, yb = 0, xr = 1, yt = 1, first = TRUE, vertical = TRUE) {
|
|
|
292 # treemap for data.tree plotting - experimental
|
|
|
293 if (first) {
|
|
|
294 plot(xl:xr, yb:yt, type = "n", axes = FALSE, xlab = "", ylab = "")
|
|
|
295 }
|
|
|
296 if (tr$nhits > 0) {
|
|
|
297 width = 2 * (6 - tr$level)
|
|
|
298 rect(xl, yb, xr, yt, col = "#00008805", lwd = width, border = "#00000080")
|
|
|
299 if (tr$level != 1) {
|
|
|
300 text(x = (xl + xr)/2, y = (yb + yt)/2, labels = tr$name, cex = width/4)
|
|
|
301 }
|
|
|
302 } else {
|
|
|
303 return(NULL)
|
|
|
304 }
|
|
|
305 Nchildren = length(tr$children)
|
|
|
306 cat("children:", tr$name, tr$level, "\n")
|
|
|
307 if (Nchildren == 0) {
|
|
|
308 return(NULL)
|
|
|
309 }
|
|
|
310 sizes = sapply(tr$children, "[[", "nhits")
|
|
|
311 rng = c(sizes, 0) / sum(sizes)
|
|
|
312 if (vertical) {
|
|
|
313 ## x does not change
|
|
|
314 xl2 = rep(xl, Nchildren)
|
|
|
315 xr2 = rep(xr, Nchildren)
|
|
|
316 yb2 = rescale(rng, c(yb, yt))[1:(Nchildren + 1)]
|
|
|
317 yt2 = rescale(rng, c(yb, yt))[-1]
|
|
|
318 } else {
|
|
|
319 yb2 = rep(yb, Nchildren)
|
|
|
320 yt2 = rep(yt, Nchildren)
|
|
|
321 xr2 = rescale(rng, c(xl, xr))[1:(Nchildren + 1)]
|
|
|
322 xl2 = rescale(rng, c(xl, xr))[-1]
|
|
|
323 }
|
|
|
324 for (i in 1:Nchildren) {
|
|
|
325 trmap(tr$children[[i]], xl2[i], yb2[i], xr2[i], yt2[i], first = FALSE, vertical = !vertical)
|
|
|
326 }
|
|
|
327
|
|
|
328 }
|
|
|
329
|
|
|
330 filter_tree = function(tr) {
|
|
|
331 ## keep only nodes with positive nhits values
|
|
|
332 ## must be used on trees where leaves were already propagated
|
|
|
333 ## using add_values_to_nodes
|
|
|
334 tr_filtered = Clone(tr)
|
|
|
335 Prune(tr_filtered, function(x) x$nhits > 0 | containLTR(x))
|
|
|
336 return(tr_filtered)
|
|
|
337 }
|
|
|
338
|
|
|
339 containLTR = function(tr){
|
|
|
340 for (i in Traverse(tr)){
|
|
|
341 if (i$features == "LTR/PBS"){
|
|
|
342 return(TRUE)
|
|
|
343 }
|
|
|
344 }
|
|
|
345 return(FALSE)
|
|
|
346 }
|
|
|
347
|
|
|
348 filter_tree2 = function(tr) {
|
|
|
349 tr_filtered = Clone(tr)
|
|
|
350 Prune(tr_filtered, function(x) x$parent$nhits > 0)
|
|
|
351 return(tr_filtered)
|
|
|
352 }
|
|
|
353
|
|
|
354 ## this is for testing purposesm in futurem each node will have its function
|
|
|
355 ## named find_best_hit, it will work like decisin tree.
|
|
|
356 ## functions will be set manually based on training data
|
|
|
357 formatWidth = function(x, w){
|
|
|
358 if (length (dim(x)) > 1){
|
|
|
359 for (i in seq_along(x)){
|
|
|
360 delta = nchar(x[,i], type="bytes") - nchar(x[,i], type="char")
|
|
|
361 ## delta correction is neccessary for utf-8 correct formating
|
|
|
362 x[,i] = sprintf(paste0("%-",w[i] + delta,"s"), x[,i])
|
|
|
363 }
|
|
|
364 return(x)
|
|
|
365 }else{
|
|
|
366 delta = nchar(x, type="bytes") - nchar(x, type="char")
|
|
|
367 sprintf(paste0("%",w + delta,"s"),
|
|
|
368 x) %>% return
|
|
|
369 }
|
|
|
370 }
|
|
|
371
|
|
|
372
|
|
|
373 find_best_hit_repeat = function(cls_tree){
|
|
|
374 ## three children:
|
|
|
375 ## repeat
|
|
|
376 ## |--rDNA
|
|
|
377 ## |--satellite this need special handling, rank 1 or 2
|
|
|
378 ## |--mobile_element
|
|
|
379 ##
|
|
|
380 ## in this case we require that satellite has proportion 0.95
|
|
|
381 ## params of good hit
|
|
|
382 min_prop_of_positive = 0.7
|
|
|
383 min_prop_second_best = 0.9
|
|
|
384 min_proportion_x_nhits = 2.5 # based on 0.05 x 50
|
|
|
385 min_satellite_proportion = 0.95
|
|
|
386 if (isLeaf(cls_tree)) {
|
|
|
387 return(cls_tree)
|
|
|
388 }
|
|
|
389 cond1 = cls_tree$nhits * cls_tree$proportion < min_proportion_x_nhits
|
|
|
390 if (cond1){
|
|
|
391 return(cls_tree)
|
|
|
392 }
|
|
|
393 nhits = sapply(cls_tree$children, "[[", "nhits")
|
|
|
394 all_hits = cls_tree$root$nhits
|
|
|
395 name_of_best_hit = cls_tree$children[[which.max(nhits)]]$name
|
|
|
396 best_hit_child = cls_tree$children[[which.max(nhits)]]
|
|
|
397 second_max_hits = ifelse(length(nhits) == 1, 0, max(nhits[-which.max(nhits)]))
|
|
|
398 cond2 = max(nhits) / sum(nhits) > min_prop_of_positive
|
|
|
399 cond3 = max(nhits) / (second_max_hits + max(nhits)) > min_prop_second_best
|
|
|
400 cond_satellite = best_hit_child$proportion > min_satellite_proportion & name_of_best_hit == "satellite"
|
|
|
401 cond_other = ! name_of_best_hit == "satellite"
|
|
|
402 cat(cond2, cond3, cond_satellite, cond_other, "\n")
|
|
|
403 if (cond2 & cond3 & (cond_satellite | cond_other)) {
|
|
|
404 ## clear case satellite or other
|
|
|
405 if ("find_best_hit" %in% names(best_hit_child)){
|
|
|
406 ## custom rules for node is defined
|
|
|
407 best_hit = best_hit_child$find_best_hit(cls_tree$children[[which.max(nhits)]])
|
|
|
408 }else{
|
|
|
409 # use generic rules
|
|
|
410 best_hit = find_best_hit(cls_tree$children[[which.max(nhits)]])
|
|
|
411 }
|
|
|
412 }else{
|
|
|
413 ## do more specific tests for rDNA, or mobile element
|
|
|
414 ## rDNA o mobile_elements must be above threshold!
|
|
|
415 cond_sat_rank = any(cls_tree$satellite$tandem_rank == 2) & cls_tree$satellite$nhits != 0
|
|
|
416 cond_rdna = cls_tree$rDNA$nhits * cls_tree$rDNA$proportion > min_proportion_x_nhits
|
|
|
417 cond_mobile = cls_tree$mobile_element$nhits * cls_tree$mobile_element$proportion > min_proportion_x_nhits
|
|
|
418 cat(cond_sat_rank, "\n")
|
|
|
419 cat(cond_rdna,"\n")
|
|
|
420 cat(cond_mobile, "\n")
|
|
|
421 if (cond_sat_rank & (cond_rdna | cond_mobile) ){
|
|
|
422 cls_tree$satellite$nhits = 0
|
|
|
423 best_hit = find_best_hit(cls_tree)
|
|
|
424 }else{
|
|
|
425 return(cls_tree)
|
|
|
426 }
|
|
|
427 }
|
|
|
428 return(best_hit)
|
|
|
429 }
|
|
|
430
|
|
|
431 find_best_hit = function(cls_tree){
|
|
|
432 ## general params of good hit
|
|
|
433 min_prop_of_positive = 0.7
|
|
|
434 min_prop_second_best = 0.8
|
|
|
435 min_proportion_x_nhits = 2.5 # based on 0.05 x 50
|
|
|
436 if (isLeaf(cls_tree)) {
|
|
|
437 return(cls_tree)
|
|
|
438 }
|
|
|
439
|
|
|
440
|
|
|
441 cond1 = cls_tree$nhits * cls_tree$proportion < min_proportion_x_nhits
|
|
|
442
|
|
|
443 if (cond1){
|
|
|
444 ## return if proportions is lower than threshold
|
|
|
445 return(cls_tree)
|
|
|
446 }
|
|
|
447
|
|
|
448 nhits = sapply(cls_tree$children, "[[", "nhits")
|
|
|
449 best_hit_child = cls_tree$children[[which.max(nhits)]]
|
|
|
450 ## more special cases:
|
|
|
451 second_max_hits = ifelse(length(nhits) == 1, 0, max(nhits[-which.max(nhits)]))
|
|
|
452 cond2 = max(nhits) / sum(nhits) > min_prop_of_positive
|
|
|
453 cond3 = max(nhits) / (second_max_hits + max(nhits)) > min_prop_second_best
|
|
|
454 if (cond2 & cond3) {
|
|
|
455 if ("find_best_hit" %in% names(best_hit_child)){
|
|
|
456 ## custom rules for node is defined
|
|
|
457 best_hit = best_hit_child$find_best_hit(cls_tree$children[[which.max(nhits)]])
|
|
|
458 }else{
|
|
|
459 # use generic rules
|
|
|
460 best_hit = find_best_hit(cls_tree$children[[which.max(nhits)]])
|
|
|
461 }
|
|
|
462 } else {
|
|
|
463 return(cls_tree)
|
|
|
464 }
|
|
|
465 return(best_hit)
|
|
|
466 }
|
|
|
467
|
|
|
468
|
|
|
469 get_annotation_groups = function(tr){
|
|
|
470 tr_all = Clone(tr)
|
|
|
471 Prune(tr_all, pruneFun=function(...)FALSE) ## prune everithing -> keep root
|
|
|
472 tr_all$name = "Unclassified repeat (No evidence)"
|
|
|
473 tr_contamination = Clone(tr)$contamination
|
|
|
474 tr_organelle = Clone(tr)$organelle
|
|
|
475 tr_repeat = Clone(tr)$"repeat"
|
|
|
476 tr_repeat$name = "Unclassified_repeat (conflicting evidences)"
|
|
|
477 return(
|
|
|
478 list(
|
|
|
479 tr_repeat = tr_repeat,
|
|
|
480 tr_organelle = tr_organelle,
|
|
|
481 tr_all = tr_all,
|
|
|
482 tr_contamination = tr_contamination
|
|
|
483 )
|
|
|
484 )
|
|
|
485 }
|
|
|
486
|
|
|
487
|
|
|
488
|
|
|
489 format_tree = function(cls_tree, ...) {
|
|
|
490 df = ToDataFrameTree(cls_tree, ...)
|
|
|
491 ## try to get rid off unicode characters
|
|
|
492 df[,1] = gsub("\U00B0", "'", df[,1], fixed = TRUE, useBytes = TRUE)
|
|
|
493 df[,1] = gsub("\U00A6", "|", df[,1], fixed = TRUE, useBytes = TRUE)
|
|
|
494 colnames(df)[1] = " " ## originally levelName
|
|
|
495 if ("proportion" %in% colnames(df)){
|
|
|
496 df$proportion = signif(df$proportion,2)
|
|
|
497 }
|
|
|
498 if ("Proportion[%]" %in% colnames(df)){
|
|
|
499 df$"Proportion[%]" = round(df$"Proportion[%]", 2)
|
|
|
500 }
|
|
|
501 ## format header
|
|
|
502 w1 = apply(df, 2, function(x) max(nchar(x)))
|
|
|
503 w2 = nchar(colnames(df))
|
|
|
504 # use whatever with is bigger for formating
|
|
|
505 w = ifelse(w1 > w2, w1, w2)
|
|
|
506
|
|
|
507 out = character()
|
|
|
508 header = mapply(FUN = function(x, y) formatWidth(x, w = y), colnames(df), w) %>%
|
|
|
509 paste0(collapse=" | ")
|
|
|
510 hr_line = gsub(".","-",header)
|
|
|
511 # create output
|
|
|
512 # classification lines
|
|
|
513 class_lines = formatWidth(df, w) %>%
|
|
|
514 apply(1, FUN = function(x) paste0(x,collapse = " | ")) %>%
|
|
|
515 paste(collapse = "\n")
|
|
|
516 paste(
|
|
|
517 header,
|
|
|
518 hr_line,
|
|
|
519 class_lines,
|
|
|
520 sep="\n") %>% return
|
|
|
521 }
|
|
|
522
|
|
|
523
|
|
|
524 make_final_annotation_template = function(
|
|
|
525 annot_attributes = list(
|
|
|
526 Nreads = 0,
|
|
|
527 Nclusters = 0,
|
|
|
528 Nsuperclusters = 0,
|
|
|
529 "Proportion[%]" = 0,
|
|
|
530 cluster_list = c())
|
|
|
531 )
|
|
|
532 {
|
|
|
533 ct = Clone(CLS_TREE)
|
|
|
534 for (i in Traverse(ct)){
|
|
|
535 for (j in names(annot_attributes)){
|
|
|
536 i[[j]] = annot_attributes[[j]]
|
|
|
537 }
|
|
|
538 }
|
|
|
539 return(ct)
|
|
|
540 }
|
|
|
541
|
|
|
542
|
|
|
543 common_ancestor = function(tr1, tr2){
|
|
|
544 a1 = tr1$Get('name', traversal = "ancestor")
|
|
|
545 a2 = tr2$Get('name', traversal = "ancestor")
|
|
|
546 ancestor = intersect(a1,a2)[1]
|
|
|
547 return(ancestor)
|
|
|
548 }
|
|
|
549
|
|
|
550 create_all_superclusters_report = function(max_supercluster = 100,
|
|
|
551 paths,
|
|
|
552 libdir,
|
|
|
553 superclusters_dir,
|
|
|
554 seqdb, hitsortdb,
|
|
|
555 classification_hierarchy_file,
|
|
|
556 HTML_LINKS)
|
|
|
557 {
|
|
|
558 ## connect to sqlite databases
|
|
|
559 ## seqdb and hitsortdb are path to sqlite files
|
|
|
560 HTML_LINKS = nested2named_list(HTML_LINKS)
|
|
|
561 paths = nested2named_list(paths)
|
|
|
562 report_file = paths[['supercluster_report_html']]
|
|
|
563 ## create SEQDB, HTISORTDB and CLS_TREE
|
|
|
564 connect_to_databases(seqdb, hitsortdb, classification_hierarchy_file)
|
|
|
565
|
|
|
566 ## append specific classication rules
|
|
|
567 CLS_TREE$"repeat"$find_best_hit = find_best_hit_repeat
|
|
|
568
|
|
|
569 ###
|
|
|
570 cat("Superclusters summary\n", file = report_file)
|
|
|
571 cat("---------------------\n\n", file = report_file, append = TRUE)
|
|
|
572 empty = rep(NA,max_supercluster)
|
|
|
573 SC_table = data.frame(
|
|
|
574 Supercluster = 1 : max_supercluster, "Number of reads" = empty, Automatic_annotation = empty,
|
|
|
575 Similarity_hits = empty, TAREAN_annotation = empty,
|
|
|
576 Clusters = empty, stringsAsFactors = FALSE, check.names = FALSE
|
|
|
577 )
|
|
|
578 SC_csv = SC_table # table as spreadsheet - no html formatings
|
|
|
579 final_cluster_annotation = make_final_annotation_template()
|
|
|
580
|
|
|
581 for (sc in 1 : max_supercluster) {
|
|
|
582 cat("supercluster: ",sc,"\n")
|
|
|
583 cls_tree = read_annotation_to_tree(sc)
|
|
|
584 sc_summary = get_supercluster_summary(sc)
|
|
|
585 add_value_to_nodes(cls_tree)
|
|
|
586 cls_tree_filtered = filter_tree(cls_tree)
|
|
|
587 best_hit = find_best_hit(cls_tree)
|
|
|
588 ## exception to decision tree put here.
|
|
|
589 ## All -> class I ?
|
|
|
590 if (best_hit$name == "All"){
|
|
|
591 ## check LTR
|
|
|
592 if (any(!(unique(cls_tree$ltr_info$ltr_detection) %in% "None"))){
|
|
|
593 ## if LTR found - move classification to Class I
|
|
|
594
|
|
|
595 LTR = FindNode(best_hit, "LTR")
|
|
|
596 nhits_without_class_I = best_hit$nhits - LTR$nhits
|
|
|
597 prop_without_class_I = best_hit$proportion - LTR$proportion
|
|
|
598 cond3 = nhits_without_class_I * prop_without_class_I < 1
|
|
|
599 cond2 = best_hit$nhits * best_hit$proportion < 0.5 # other hits weak
|
|
|
600 cond1 = LTR$nhits >= (0.95 * best_hit$nhits) # hits are pre in sub class_I
|
|
|
601
|
|
|
602 if (cond1 | cond2 | cond3){
|
|
|
603 best_hit = LTR
|
|
|
604 }
|
|
|
605 }
|
|
|
606 }else{
|
|
|
607 ## check is conflict os best hit with LTR/PBS exists
|
|
|
608 if (any(!(unique(cls_tree$ltr_info$ltr_detection) %in% "None"))){
|
|
|
609 LTR = FindNode(cls_tree, "LTR")
|
|
|
610 ca_name = common_ancestor(LTR, best_hit)
|
|
|
611 if (ca_name != "LTR"){
|
|
|
612 ## reclassify
|
|
|
613 best_hit = FindNode(cls_tree, ca_name)
|
|
|
614 }
|
|
|
615 }
|
|
|
616 }
|
|
|
617
|
|
|
618 best_hit_name = best_hit$name
|
|
|
619 best_hit_path = paste(best_hit$path, collapse = "/")
|
|
|
620 ## add best hit do database:
|
|
|
621 sqlcmd = paste0(
|
|
|
622 "UPDATE cluster_info SET supercluster_best_hit = \"", best_hit_path, "\" ",
|
|
|
623 " WHERE supercluster = ", sc
|
|
|
624 )
|
|
|
625 dbExecute(HITSORTDB, sqlcmd)
|
|
|
626 # add annotation annotation summary stored in final_cluster_annotation - summing up
|
|
|
627 best_hit_node = FindNode(final_cluster_annotation, best_hit_name)
|
|
|
628 best_hit_node$Nsuperclusters = best_hit_node$Nsuperclusters + 1
|
|
|
629 best_hit_node$Nclusters = best_hit_node$Nclusters + sc_summary$Nclusters
|
|
|
630 best_hit_node$Nreads = best_hit_node$Nreads + sc_summary$Nreads
|
|
|
631 best_hit_node$"Proportion[%]" = best_hit_node$"Proportion[%]" + sc_summary$proportion*100
|
|
|
632 best_hit_node$cluster_list = append(best_hit_node$cluster_list, sc_summary$cluster_list)
|
|
|
633 best_hit_label = ifelse(best_hit_name == "All", "NA", best_hit_name)
|
|
|
634 SC_csv$Automatic_annotation[sc] = best_hit_label
|
|
|
635 SC_csv$"Number of reads"[sc] = SC_table$"Number of reads"[sc] = cls_tree$size
|
|
|
636 SC_csv$Similarity_hits[sc] = format_tree(cls_tree_filtered,
|
|
|
637 "nhits", "proportion", "features")
|
|
|
638 SC_table$Similarity_hits[sc] = format_tree(cls_tree_filtered,
|
|
|
639 "nhits", "proportion", "features") %>%
|
|
|
640 preformatted
|
|
|
641
|
|
|
642 SC_table$Automatic_annotation[sc] = best_hit_label
|
|
|
643
|
|
|
644 G = get_supercluster_graph(
|
|
|
645 sc=sc, seqdb = seqdb,hitsortdb = hitsortdb,
|
|
|
646 classification_hierarchy_file = classification_hierarchy_file
|
|
|
647 )
|
|
|
648 ## append tarean annotation
|
|
|
649 clist = paste(V(G)$name, collapse =",")
|
|
|
650 cl_rank = dbGetQuery(HITSORTDB,
|
|
|
651 paste("SELECT [index], tandem_rank FROM cluster_info WHERE [index] IN (",clist,
|
|
|
652 ") AND tandem_rank IN (1,2)"))
|
|
|
653 ## add information about TR clusters is any
|
|
|
654 if (nrow(cl_rank) > 0){
|
|
|
655 tarean_annot = paste(
|
|
|
656 cl_rank$index, "-",
|
|
|
657 RANKS_TANDEM[cl_rank$tandem_rank]
|
|
|
658 )
|
|
|
659 SC_csv$TAREAN_annotation[sc] = paste(tarean_annot,collapse="\n")
|
|
|
660 SC_table$TAREAN_annotation[sc] = mapply(
|
|
|
661 hwrite, tarean_annot,
|
|
|
662 link = sprintf(HTML_LINKS$ROOT_TO_TAREAN,cl_rank$index)) %>% paste(collapse="<br>")
|
|
|
663 }
|
|
|
664 ## add links to clusters
|
|
|
665 SC_csv$Clusters[sc] = paste((V(G)$name), collapse = ",")
|
|
|
666 links = mapply(hwrite, V(G)$name, link = sprintf(HTML_LINKS$ROOT_TO_CLUSTER, as.integer(V(G)$name)))
|
|
|
667 SC_table$Clusters[sc] = paste(links, collapse =", ")
|
|
|
668 create_single_supercluster_report(G, superclusters_dir)
|
|
|
669 }
|
|
|
670
|
|
|
671
|
|
|
672
|
|
|
673 ## add html links
|
|
|
674 SC_table$Supercluster = mapply(hwrite, SC_table$Supercluster,
|
|
|
675 link = sprintf(HTML_LINKS$ROOT_TO_SUPERCLUSTER,
|
|
|
676 SC_table$Supercluster))
|
|
|
677
|
|
|
678
|
|
|
679
|
|
|
680 write.table(SC_csv, file = paths[["superclusters_csv_summary"]],
|
|
|
681 sep = "\t", row.names = FALSE)
|
|
|
682
|
|
|
683
|
|
|
684 DT_instance = datatable(SC_table, escape = FALSE, options = DT_OPTIONS)
|
|
|
685 saveWidget(DT_instance, file = normalizePath(report_file),
|
|
|
686 normalizePath(libdir), selfcontained = FALSE)
|
|
|
687 add_preamble(normalizePath(report_file),
|
|
|
688 preamble='<h2>Supercluster annotation</h2> <p><a href="documentation.html#superclust"> For table legend see documentation. <a> </p>')
|
|
|
689
|
|
|
690 annot_groups = get_annotation_groups(final_cluster_annotation)
|
|
|
691 saveRDS(object=annot_groups, file = paths[['repeat_annotation_summary_rds']])
|
|
|
692 final_cluster_annotation_formated = paste(
|
|
|
693 sapply(annot_groups, format_tree, 'Proportion[%]', 'Nsuperclusters', 'Nclusters', "Nreads"),
|
|
|
694 collapse = "\n\n\n"
|
|
|
695 )
|
|
|
696 ## TODO add singleton counts, total counts and extra text - csv output
|
|
|
697 ## export annotation summary
|
|
|
698 html_summary = start_html(paths[['summarized_annotation_html']], gsub("PAGE_TITLE", "Repeat Annotation Summary", HTMLHEADER))
|
|
|
699 html_summary("Repeat annotation summary", HTML.title, HR=2)
|
|
|
700 html_summary("This table summarizes the automatic annotations of superclusters that should be verified and manually corrected if necessary. Thus, the table should not be used as the final output of the analysis without critical evaluation.", HTML.title, HR=3)
|
|
|
701 html_summary(preformatted(final_cluster_annotation_formated), cat)
|
|
|
702
|
|
|
703
|
|
|
704
|
|
|
705 return()
|
|
|
706 }
|
|
|
707 ## for testing
|
|
|
708 if (FALSE){
|
|
|
709 create_supercluster_report(1:100, report_file, seqdb, hitsortdb, class_file)
|
|
|
710 }
|
|
|
711
|
|
|
712 create_single_supercluster_report = function(G, superclusters_dir){
|
|
|
713 sc_dir = paste0(superclusters_dir, "/dir_SC",sprintf("%04d", G$name))
|
|
|
714 htmlfile = paste0(sc_dir, "/index.html")
|
|
|
715 dir.create(sc_dir)
|
|
|
716 title = paste("Supercluster no. ",G$name)
|
|
|
717 htmlheader = gsub("PAGE_TITLE", title, HTMLHEADER)
|
|
|
718 cat(htmlheader, file = htmlfile)
|
|
|
719 file.copy(paste0(WD,"/style1.css"), sc_dir)
|
|
|
720 HTML.title(title, file = htmlfile)
|
|
|
721 HTML("Simmilarity hits (only hits with proportion above 0.1% in at least one cluster are shown) ", file = htmlfile)
|
|
|
722 img_file = paste0(sc_dir,"/SC",sprintf("%0d",G$name), ".png")
|
|
|
723 png(filename = img_file, width = 1400, height=1200)
|
|
|
724 coords = plot_supercluster(G = G)
|
|
|
725 dev.off()
|
|
|
726 ## basename - make link relative
|
|
|
727 href = paste0("../../clusters/dir_CL",sprintf("%04d", as.numeric(V(G)$name)),"/index.html")
|
|
|
728 html_insert_image(basename(img_file), htmlfile,coords=coords, href = href)
|
|
|
729
|
|
|
730 if (is_comparative()){
|
|
|
731 HTML("Comparative analysis", file = htmlfile)
|
|
|
732 img_file = paste0(sc_dir,"/SC",sprintf("%0d",G$name), "comparative.png")
|
|
|
733 png(filename = img_file, width = 1400, height=1200)
|
|
|
734 coords = plot_supercluster(G = G, "comparative")
|
|
|
735 mtext("comparative analysis")
|
|
|
736 dev.off()
|
|
|
737 ## basename - make link relative
|
|
|
738 href = paste0("../../clusters/dir_CL",sprintf("%04d", as.numeric(V(G)$name)),"/index.html")
|
|
|
739 html_insert_image(basename(img_file), htmlfile,coords=coords, href = href)
|
|
|
740
|
|
|
741 }
|
|
|
742 }
|
|
|
743
|
|
|
744 html_insert_image = function(img_file, htmlfile, coords=NULL, href=NULL) {
|
|
|
745 img_tag = sprintf('<p> <img src="%s" usemap="#clustermap" > \n</p>\n<br>\n', img_file)
|
|
|
746 cat(img_tag, file = htmlfile, append = TRUE)
|
|
|
747 if (!is.null(coords)){
|
|
|
748 formated_links = character()
|
|
|
749 for (i in seq_along(href)){
|
|
|
750 formated_links[i] = sprintf('<area shape="circle" coords="%f,%f,%f" href="%s" >\n',
|
|
|
751 coords[i,1],coords[i,2],coords[,3],href[i])
|
|
|
752 }
|
|
|
753 cat('<map name="clustermap" >\n',
|
|
|
754 formated_links,"</map>\n", file=htmlfile, append = TRUE)
|
|
|
755
|
|
|
756 }
|
|
|
757 }
|
|
|
758
|
|
|
759 html_insert_floating_image = function(img_file,
|
|
|
760 htmlfile,
|
|
|
761 width=NULL,
|
|
|
762 title= "",
|
|
|
763 footer = ""
|
|
|
764 ){
|
|
|
765 if (is.null(width)){
|
|
|
766 width=""
|
|
|
767 }
|
|
|
768 tag = paste(
|
|
|
769 '<div class="floating_img">\n',
|
|
|
770 ' <p>', title,'</p>\n',
|
|
|
771 ' <a href=',img_file, ">",
|
|
|
772 ' <img src="',img_file, '" alt="image" width="',width,'">\n',
|
|
|
773 ' </a>',
|
|
|
774 ' <p> ',footer,'</p>\n',
|
|
|
775 '</div>\n'
|
|
|
776 ,sep = "")
|
|
|
777 cat(tag, file = htmlfile, append = TRUE)
|
|
|
778 }
|
|
|
779
|
|
|
780 get_supercluster_graph = function(sc, seqdb, hitsortdb, classification_hierarchy_file){
|
|
|
781 ## sc - superclusted index
|
|
|
782 ## seqdb, hitsortdb - path to sqlite databases
|
|
|
783 ## classificationn_hierarchy_file - path data.tree rds file
|
|
|
784 ## SIZE of image
|
|
|
785 SIZE = 1000
|
|
|
786 connect_to_databases(seqdb, hitsortdb, classification_hierarchy_file)
|
|
|
787 clusters = dbGetQuery(HITSORTDB,
|
|
|
788 paste0("SELECT cluster, size FROM superclusters WHERE supercluster=",sc)
|
|
|
789 ) %>% as.data.frame
|
|
|
790 ## string for sql query:
|
|
|
791 cluster_list = paste0( "(", paste0(clusters$cluster, collapse = ","),")")
|
|
|
792 supercluster_ncol = dbGetQuery(HITSORTDB,
|
|
|
793 paste("SELECT c1,c2,w, k FROM cluster_mate_bond where c1 IN ",
|
|
|
794 cluster_list, " AND c2 IN ", cluster_list, "AND k > 0.1"
|
|
|
795 )
|
|
|
796 )
|
|
|
797 if (nrow(supercluster_ncol)>0){
|
|
|
798 ## at least two clusters
|
|
|
799 G = graph.data.frame(supercluster_ncol, directed=FALSE)
|
|
|
800 L = layout_with_kk(G)
|
|
|
801 ## layout is rescaled for easier linking html
|
|
|
802 L[,1] = scales::rescale(L[,1], to = c(1,SIZE) )
|
|
|
803 L[,2] = scales::rescale(L[,2], to = c(1,SIZE) )
|
|
|
804 G$L = L
|
|
|
805 }else{
|
|
|
806 ## only one cluster in supercluster
|
|
|
807 G = graph.full(n = 1)
|
|
|
808 V(G)$name = as.character(clusters$cluster)
|
|
|
809 G$L = matrix(c(SIZE/2, SIZE/2), ncol=2)
|
|
|
810 }
|
|
|
811 G = set_vertex_attr(G, "label", value = paste0("CL",names(V(G))))
|
|
|
812 G$name=sc ## names is supercluster
|
|
|
813 # create annotation of individual clusters, will be attached as node attribudes
|
|
|
814 annot = get_cluster_annotation_summary(clusters)
|
|
|
815 ## annotation of nodes(clusters)
|
|
|
816 G = set_vertex_attr(G, "annotation", value = annot$clusters[names(V(G))])
|
|
|
817 G = set_vertex_attr(G, "size",
|
|
|
818 value = clusters$size[match(names(V(G)), as.character(clusters$cluster))])
|
|
|
819 G$annotation = annot$supercluster
|
|
|
820 G$clusters = clusters
|
|
|
821 # TODO if comparative analysis - add proportion of species
|
|
|
822 if (is_comparative()){
|
|
|
823 comparative_counts = get_cluster_comparative_counts(cluster_list)
|
|
|
824 NACluster = comparative_counts$supercluster
|
|
|
825 # some clusters do not have comparative data - they are bellow threshold!
|
|
|
826 NACluster[,] = NA
|
|
|
827 counts_adjusted = comparative_counts$clusters[names(V(G))]
|
|
|
828 for (i in names(V(G))){
|
|
|
829 if(is.null(counts_adjusted[[i]])){
|
|
|
830 ## get count from database
|
|
|
831 seqid = dbGetQuery(HITSORTDB,
|
|
|
832 paste(
|
|
|
833 "SELECT vertexname FROM vertex_cluster where cluster = ",
|
|
|
834 i)
|
|
|
835 )$vertexname
|
|
|
836 codes = get_comparative_codes()$prefix
|
|
|
837 x = table(factor(substring(seqid, 1,nchar(codes[1])), levels = codes))
|
|
|
838 counts_adjusted[[i]] = data.frame(
|
|
|
839 Freq = as.numeric(x),
|
|
|
840 proportion = as.numeric(x/sum(x)),
|
|
|
841 row.names = names(x))
|
|
|
842
|
|
|
843 }
|
|
|
844 }
|
|
|
845
|
|
|
846 G = set_vertex_attr(G, "comparative_counts",
|
|
|
847 value = counts_adjusted[V(G)$name]
|
|
|
848 )
|
|
|
849 # for whole supercluster
|
|
|
850 G$comparative = comparative_counts$supercluster
|
|
|
851 }
|
|
|
852
|
|
|
853 return(G)
|
|
|
854 }
|
|
|
855
|
|
|
856 get_cluster_comparative_counts = function(cluster_list){
|
|
|
857 counts = dbGetQuery(
|
|
|
858 HITSORTDB,
|
|
|
859 paste0("SELECT * FROM comparative_counts WHERE clusterindex IN",
|
|
|
860 cluster_list
|
|
|
861 )
|
|
|
862 )
|
|
|
863 comparative_counts = apply(counts, 1, function(x)
|
|
|
864 y = data.frame(Freq = x[-1], proportion = x[-1]/sum(x[-1]))
|
|
|
865 )
|
|
|
866 names(comparative_counts) = counts$clusterindex
|
|
|
867 total_counts = data.frame(
|
|
|
868 Freq = colSums(counts[,-1]),
|
|
|
869 proportion = colSums(counts[,-1])/sum(counts[,-1]))
|
|
|
870 return(
|
|
|
871 list(
|
|
|
872 clusters = comparative_counts,
|
|
|
873 supercluster = total_counts
|
|
|
874 )
|
|
|
875 )
|
|
|
876
|
|
|
877 }
|
|
|
878
|
|
|
879 get_cluster_annotation_summary = function(clusters){
|
|
|
880 ## clusters - table of clusters, col: cluster, size
|
|
|
881 annot_list = apply(clusters,1,FUN = function(x)summarize_annotation(cluster_annotation(x[1]),x[2]))
|
|
|
882 ## if empty - not annotation at all
|
|
|
883 if (sum(sapply(annot_list, nrow)) == 0){
|
|
|
884 return(list (clusters = NULL,
|
|
|
885 superclusters = NULL
|
|
|
886 ))
|
|
|
887 }
|
|
|
888 names(annot_list) = as.character(clusters$cluster)
|
|
|
889 annot_all = do.call(rbind,annot_list)
|
|
|
890 total = by(annot_all$Freq,INDICES=with(annot_all, paste(cl_string,domain)), FUN=sum, simplify=TRUE)
|
|
|
891 total_df = data.frame(Freq = c(total), proportion = c(total) /sum(clusters$size))
|
|
|
892 ## for ploting it need to contain all annot categories as in s upercluster
|
|
|
893 annot_list_full = list()
|
|
|
894 for (i in seq_along(annot_list)){
|
|
|
895 annot_list_full[[i]] = total_df
|
|
|
896 for (j in rownames(total_df)){
|
|
|
897 if (!j %in% rownames(annot_list[[i]])){
|
|
|
898 annot_list_full[[i]][j,c('Freq','proportion')] = c(0,0)
|
|
|
899 }else{
|
|
|
900 annot_list_full[[i]][j,c('Freq','proportion')] = annot_list[[i]][j,c('Freq','proportion')]
|
|
|
901 }
|
|
|
902 }
|
|
|
903 }
|
|
|
904 names(annot_list_full) = names(annot_list)
|
|
|
905 return(list (clusters = annot_list_full,
|
|
|
906 supercluster = total_df
|
|
|
907 )
|
|
|
908 )
|
|
|
909 }
|
|
|
910
|
|
|
911 plot_supercluster = function(G,color_scheme=c("annotation","comparative")){
|
|
|
912 ## create plot in coords y: 1-1200, x: 1 - 1400
|
|
|
913 ## this is fixed for href to work in exact image areas!
|
|
|
914 color_scheme = match.arg(color_scheme)
|
|
|
915 LIMS0 = matrix(c (1,1000,1,1000), ncol = 2)
|
|
|
916 OFFSET_H = 100; OFFSET_W = 200
|
|
|
917 lims = LIMS0 + c(0,OFFSET_W * 2, 0, OFFSET_H * 2)
|
|
|
918 ## use full range
|
|
|
919 par(mar=c(0,0,0,0),xaxs="i", yaxs="i")
|
|
|
920 ## proportion of repeats
|
|
|
921 if (color_scheme == "annotation"){
|
|
|
922 prop = sapply (V(G)$annotation, FUN = function(x)x$proportion)
|
|
|
923 brv = c("#BBBBBB",unique_colors)
|
|
|
924 }else{
|
|
|
925 prop = sapply (V(G)$comparative_counts, FUN = function(x)x$proportion)
|
|
|
926 brv = unique_colors
|
|
|
927 }
|
|
|
928 ## handle special cases - single cluster with annot or single annotation group
|
|
|
929 if (is.null(dim(prop))){
|
|
|
930 prop=matrix(prop, nrow = 1)
|
|
|
931 }
|
|
|
932 if (length(prop)==0){
|
|
|
933 ## special case - not annotation at all
|
|
|
934 prop = matrix(rep(1,vcount(G)),ncol = vcount(G), dimnames=list("NAN", NULL))
|
|
|
935 }else{
|
|
|
936 if (color_scheme == "annotation"){
|
|
|
937 rownames(prop) = rownames(G$annotation)
|
|
|
938 }else{
|
|
|
939 rownames(prop) = rownames(G$comparative)
|
|
|
940 }
|
|
|
941 ## reorder by prop
|
|
|
942 if (color_scheme =="annotation"){
|
|
|
943 prop = prop[order(prop[,1], decreasing = TRUE),,drop=FALSE]
|
|
|
944 NAN = 1 - colSums(prop)
|
|
|
945 prop = rbind(NAN, prop)
|
|
|
946 brv = c("#BBBBBB",unique_colors)
|
|
|
947 }else{
|
|
|
948 }
|
|
|
949 }
|
|
|
950 L = G$L + c(OFFSET_H)
|
|
|
951 ## for href - necessary convert coordinater - image is flipped vertically
|
|
|
952 include = rowSums(prop>0.001, na.rm = TRUE)>=1
|
|
|
953 include[1] = TRUE
|
|
|
954 prop = prop[include, , drop=FALSE]
|
|
|
955 plot(0,0,type='n', , xlim = lims[,1], ylim = lims[,2])
|
|
|
956 plot_edges(G, L, lwd = 8)
|
|
|
957 RADIUS = radius_size(G)
|
|
|
958 coords = cbind(G$L[,1,drop=FALSE] + 100, 1100 - G$L[,2,drop=FALSE], RADIUS)
|
|
|
959 pieScatter(L[,1], L[,2], t(prop), col = brv,
|
|
|
960 radius = RADIUS, new.plot = FALSE, border=NA)
|
|
|
961 ## add ledend
|
|
|
962 legend("topright", legend = rownames(prop),
|
|
|
963 col = brv, pch = 15, cex= 1.5, pt.cex= 3)
|
|
|
964 ## add cluster labels
|
|
|
965 text(x = L[,1], y = L[,2], labels = paste(V(G)$label,"\n",V(G)$size))
|
|
|
966 return(coords)
|
|
|
967 }
|
|
|
968
|
|
|
969
|
|
|
970 radius_size = function(G){
|
|
|
971 if (vcount(G) == 1) RADIUS=120
|
|
|
972 if (vcount(G) == 2) RADIUS=50
|
|
|
973 if (vcount(G) %in% 3:8) RADIUS=40
|
|
|
974 if (vcount(G) > 8) RADIUS=30
|
|
|
975 return(RADIUS)
|
|
|
976 }
|
|
|
977
|
|
|
978
|
|
|
979 plot_edges = function(G,L,col="#33000040", lwd = 1){
|
|
|
980 e = get.edgelist(G, names = F)
|
|
|
981 X0 = L[e[, 1], 1]
|
|
|
982 Y0 = L[e[, 1], 2]
|
|
|
983 X1 = L[e[, 2], 1]
|
|
|
984 Y1 = L[e[, 2], 2]
|
|
|
985 segments(X0, Y0, X1, Y1, lwd = lwd,col = col)
|
|
|
986 }
|
|
|
987
|
|
|
988
|
|
|
989 pieScatter=function(x, y, prop,radius=1, col=NULL, edges=100, new.plot = TRUE, ...){
|
|
|
990 if (new.plot){
|
|
|
991 plot(x,y,type='n')
|
|
|
992 }
|
|
|
993 for (i in seq_along(x)){
|
|
|
994 p=prop[i,]
|
|
|
995 if (length(radius)==1){
|
|
|
996 r=radius
|
|
|
997 }else{
|
|
|
998 r=radius[i]
|
|
|
999 }
|
|
|
1000 include = p != 0
|
|
|
1001 floating.pie(x[i], y[i],p[include],
|
|
|
1002 col=col[include], radius=r,edges=50,...)
|
|
|
1003 }
|
|
|
1004 }
|
|
|
1005
|
|
|
1006 unique_colors = c("#00FF00", "#0000FF", "#FF0000",
|
|
|
1007 "#FFA6FE", "#FFDB66", "#006401", "#010067", "#95003A",
|
|
|
1008 "#007DB5", "#FF00F6", "#FFEEE8", "#774D00", "#90FB92", "#01FFFE",
|
|
|
1009 "#0076FF", "#D5FF00", "#FF937E", "#6A826C", "#FF029D",
|
|
|
1010 "#FE8900", "#7A4782", "#7E2DD2", "#85A900", "#FF0056",
|
|
|
1011 "#A42400", "#00AE7E", "#683D3B", "#BDC6FF", "#263400",
|
|
|
1012 "#BDD393", "#00B917", "#9E008E", "#001544", "#C28C9F",
|
|
|
1013 "#FF74A3", "#01D0FF", "#004754", "#E56FFE", "#788231",
|
|
|
1014 "#0E4CA1", "#91D0CB", "#BE9970", "#968AE8", "#BB8800",
|
|
|
1015 "#43002C", "#DEFF74", "#00FFC6", "#FFE502", "#620E00",
|
|
|
1016 "#008F9C", "#98FF52", "#7544B1", "#B500FF", "#00FF78",
|
|
|
1017 "#FF6E41", "#005F39", "#6B6882", "#5FAD4E", "#A75740",
|
|
|
1018 "#A5FFD2", "#FFB167", "#009BFF", "#E85EBE")
|
|
|
1019
|
|
|
1020 if (FALSE){
|
|
|
1021 ## For testing
|
|
|
1022 for ( i in 1:100){
|
|
|
1023 G = get_supercluster_graph(
|
|
|
1024 sc=i, seqdb = seqdb,
|
|
|
1025 hitsortdb = hitsortdb, classification_hierarchy_file = class_file
|
|
|
1026 )
|
|
|
1027 png(paste0("/mnt/raid/users/petr/tmp/test.figure_",i,".png"), width = 1400, height=1000)
|
|
|
1028 plot_supercluster(G)
|
|
|
1029 dev.off()
|
|
|
1030 print(i)
|
|
|
1031 }
|
|
|
1032
|
|
|
1033 }
|
|
|
1034
|
|
|
1035 plotg = function(GG,LL,wlim=NULL,...){
|
|
|
1036 ## quick ploting function
|
|
|
1037 e = get.edgelist(GG, names=F)
|
|
|
1038 w = E(GG)$weight
|
|
|
1039 if (!is.null(wlim)) {e = e[w > wlim,]; w = w[w > wlim]}
|
|
|
1040 X0 = LL[e[,1],1]
|
|
|
1041 Y0 = LL[e[,1],2]
|
|
|
1042 X1 = LL[e[,2],1]
|
|
|
1043 Y1 = LL[e[,2],2]
|
|
|
1044 plot(range(LL[,1]),range(LL[,2]),xlab="",ylab="",axes=FALSE,type="n",...)
|
|
|
1045 brv = 'grey'
|
|
|
1046 segments(X0,Y0,X1,Y1,lwd=.5,col=brv)
|
|
|
1047 points(LL,pch=18,cex=.8,...)
|
|
|
1048 }
|
|
|
1049
|
|
|
1050
|
|
|
1051
|
|
|
1052
|
|
|
1053 get_cluster_info = function(index){
|
|
|
1054 ## must be connected to hitsort database ( HITSORT)
|
|
|
1055 x = dbGetQuery(HITSORTDB,
|
|
|
1056 paste0("SELECT * FROM cluster_info WHERE [index]=",index))
|
|
|
1057 return(as.list(x))
|
|
|
1058 }
|
|
|
1059
|
|
|
1060 get_supercluster_info = function(sc){
|
|
|
1061 ## must be connected to hitsort database ( HITSORT)
|
|
|
1062 x = dbGetQuery(HITSORTDB,
|
|
|
1063 paste0("SELECT * FROM cluster_info WHERE supercluster=",sc))
|
|
|
1064 return(x)
|
|
|
1065 }
|
|
|
1066
|
|
|
1067 get_supercluster_summary = function(sc){
|
|
|
1068 info = get_supercluster_info(sc)
|
|
|
1069 Nclusters = nrow(info)
|
|
|
1070 Nreads = supercluster_size(sc)
|
|
|
1071 N = dbGetQuery(SEQDB, "SELECT count(*) from sequences")
|
|
|
1072 proportion = Nreads/N
|
|
|
1073 cluster_list = info$index
|
|
|
1074 return(list(
|
|
|
1075 Nreads = Nreads,
|
|
|
1076 Nclusters = Nclusters,
|
|
|
1077 proportion = proportion,
|
|
|
1078 cluster_list = cluster_list
|
|
|
1079 ))
|
|
|
1080 }
|
|
|
1081
|
|
|
1082
|
|
|
1083 get_cluster_connection_info = function(index, search_type=c("pair", "similarity")){
|
|
|
1084 search_type = match.arg(search_type)
|
|
|
1085 if (search_type == "pair"){
|
|
|
1086 x = dbGetQuery(HITSORTDB,
|
|
|
1087 sprintf(
|
|
|
1088 "SELECT * FROM cluster_mate_connections WHERE c1==%d OR c2=%d",
|
|
|
1089 index, index)
|
|
|
1090 )
|
|
|
1091 if (nrow(x) == 0 ){
|
|
|
1092 ## no connections
|
|
|
1093 return(NULL)
|
|
|
1094 }
|
|
|
1095 other_cl = ifelse(x$c1 == index, x$c2, x$c1)
|
|
|
1096 out = data.frame(cl = other_cl, N = x$N)
|
|
|
1097 out$ids = unlist(strsplit(x$ids, split = ", "))
|
|
|
1098 k = dbGetQuery(HITSORTDB,
|
|
|
1099 sprintf(
|
|
|
1100 "SELECT * FROM cluster_mate_bond WHERE c1==%d OR c2=%d",
|
|
|
1101 index, index)
|
|
|
1102 )
|
|
|
1103 if (k$c1 != x$c1 || k$c2 !=x$c2){
|
|
|
1104 ## tables not in same order
|
|
|
1105 stop("tables cluster_mate_bond and cluster_mate_connection are not in the same order")
|
|
|
1106 }
|
|
|
1107 out$k = k$k
|
|
|
1108 return(out)
|
|
|
1109 }else{
|
|
|
1110 x = dbGetQuery(HITSORTDB,
|
|
|
1111 sprintf(
|
|
|
1112 "SELECT * FROM cluster_connections WHERE c1==%d OR c2=%d",
|
|
|
1113 index, index)
|
|
|
1114 )
|
|
|
1115 if (nrow(x) == 0 ){
|
|
|
1116 ## no connections
|
|
|
1117 return(NULL)
|
|
|
1118 }
|
|
|
1119 other_cl = ifelse(x$c1 == index, x$c2, x$c1)
|
|
|
1120 out = data.frame(cl = other_cl, N = x$"count(*)")
|
|
|
1121 return(out)
|
|
|
1122 }
|
|
|
1123 }
|
|
|
1124
|
|
|
1125
|
|
|
1126 format_clinfo=function(x){
|
|
|
1127 include = c(
|
|
|
1128 "size",
|
|
|
1129 "size_real",
|
|
|
1130 "ecount",
|
|
|
1131 "supercluster",
|
|
|
1132 "annotations_summary",
|
|
|
1133 "ltr_detection",
|
|
|
1134 "pair_completeness",
|
|
|
1135 "TR_score",
|
|
|
1136 "TR_monomer_length",
|
|
|
1137 "loop_index",
|
|
|
1138 "satellite_probability",
|
|
|
1139 "TR_consensus",
|
|
|
1140 "tandem_rank",
|
|
|
1141 "orientation_score"
|
|
|
1142 )
|
|
|
1143 new_names = c(
|
|
|
1144 "TR_consensus" = "TAREAN consensus",
|
|
|
1145 "tandem_rank" = "TAREAN_annotation",
|
|
|
1146 "ltr_detection" = "PBS/LTR",
|
|
|
1147 'size' = 'number of reads in the graph',
|
|
|
1148 'size_real' = 'the total number of reads in the cluster',
|
|
|
1149 'ecount' = "number of edges of the graph:",
|
|
|
1150 "annotations_summary" = "similarity based annotation"
|
|
|
1151 )
|
|
|
1152 values = unlist(x[include]) %>% gsub("\n","<br>", .)
|
|
|
1153 include = replace(include, match(names(new_names),include), new_names)
|
|
|
1154 names(values) = include
|
|
|
1155 ## replace tandem rank with annotation description
|
|
|
1156 values["TAREAN_annotation"] = RANKS_TANDEM[values["TAREAN_annotation"]]
|
|
|
1157 ## create HTML table without header
|
|
|
1158 desc = df2html(data.frame(names(values), values))
|
|
|
1159 return(desc)
|
|
|
1160 }
|
|
|
1161
|
|
|
1162 create_cluster_report = function(index,
|
|
|
1163 seqdb, hitsortdb,
|
|
|
1164 classification_hierarchy_file,
|
|
|
1165 HTML_LINKS){
|
|
|
1166 ## index - index of cluster, integer
|
|
|
1167 ## seqdb, hitsortdb - sqlite databases
|
|
|
1168 ## classification_hierarchy_file - rds file with classification
|
|
|
1169 ## other information about cluster is collected from HITSORTDB
|
|
|
1170 connect_to_databases(seqdb, hitsortdb, classification_hierarchy_file)
|
|
|
1171 ## basic info about cluster
|
|
|
1172 clinfo = get_cluster_info(index)
|
|
|
1173 HTML_LINKS = nested2named_list(HTML_LINKS)
|
|
|
1174 PNGWIDTH = 1000 # this must be kept so than link in png work correctly
|
|
|
1175 PNGHEIGHT = 1000
|
|
|
1176 LEGENDWIDTH = 300
|
|
|
1177 PS = 30 ## pointsize for plotting
|
|
|
1178 MAX_N = 10 ## to show mates or similarities connections
|
|
|
1179
|
|
|
1180 #################################################################################
|
|
|
1181 ## create HTML title:
|
|
|
1182 title = paste("Cluster no. ",index)
|
|
|
1183 htmlheader = gsub("PAGE_TITLE", title, HTMLHEADER)
|
|
|
1184 cat(htmlheader, file = clinfo$html_report_main)
|
|
|
1185 file.copy(paste0(WD,"/style1.css"), clinfo$dir)
|
|
|
1186 HTML.title(title, file = clinfo$html_report_main)
|
|
|
1187 ## print basic info to HTML header:
|
|
|
1188 ## make link back to cluster table
|
|
|
1189 hwrite("Go back to cluster table <br><br>\n", link=HTML_LINKS$CLUSTER_TO_CLUSTER_TABLE) %>%
|
|
|
1190 cat(file = clinfo$html_report_main, append =TRUE)
|
|
|
1191 ## create ling to corresponding supercluster
|
|
|
1192 cat("Cluster is part of ",
|
|
|
1193 hwrite (paste("supercluster: ", clinfo$supercluster),
|
|
|
1194 link=sprintf(HTML_LINKS$CLUSTER_TO_SUPERCLUSTER, clinfo$supercluster)),
|
|
|
1195 "\n<br><br>\n",
|
|
|
1196 file = clinfo$html_report_main, append =TRUE)
|
|
|
1197 HTML.title("Cluster characteristics:", HR = 3, file = clinfo$html_report_main)
|
|
|
1198 cat(format_clinfo(clinfo), file = clinfo$html_report_main, append =TRUE)
|
|
|
1199 #################################################################################
|
|
|
1200 ## add link to image with graph layout with domains
|
|
|
1201 fs = list(
|
|
|
1202 graph_domains = paste0(clinfo$html_report_files,"/graph_domains.png"),
|
|
|
1203 graph_mates = paste0(clinfo$html_report_files,"/graph_mates%04d_%04d.png"),
|
|
|
1204 graph_base = paste0(clinfo$html_report_files,"/graph_base.png"),
|
|
|
1205 graph_comparative = paste0(clinfo$html_report_files,"/graph_comparative.png")
|
|
|
1206 )
|
|
|
1207 fs_relative = list(
|
|
|
1208 graph_domains = paste0(basename(clinfo$html_report_files),"/graph_domains.png"),
|
|
|
1209 graph_comparative = paste0(basename(clinfo$html_report_files),"/graph_comparative.png"),
|
|
|
1210 graph_mates = paste0(basename(clinfo$html_report_files),"/graph_mates%04d_%04d.png")
|
|
|
1211 )
|
|
|
1212
|
|
|
1213 dir.create(clinfo$html_report_files)
|
|
|
1214 load(clinfo$graph_file)
|
|
|
1215 annot = cluster_annotation(index)
|
|
|
1216 annot_summary = summarize_annotation(annot, clinfo$size_real)
|
|
|
1217 vertex_colors = annot2colors(annot,ids = V(GL$G)$name)
|
|
|
1218
|
|
|
1219 color_proportions = sort(table(vertex_colors$color_table)) / clinfo$size_real
|
|
|
1220 ## hits with smaller proportion are not shown!
|
|
|
1221 exclude_colors = names(color_proportions)[color_proportions < 0.001]
|
|
|
1222 vertex_colors$color_table[vertex_colors$color_table %in% exclude_colors] = "#00000050"
|
|
|
1223 vertex_colors$legend = vertex_colors$legend[!vertex_colors$legend %in% exclude_colors]
|
|
|
1224
|
|
|
1225 if (is_comparative()){
|
|
|
1226 comp_codes = get_comparative_codes()$prefix
|
|
|
1227 colors = unique_colors[seq_along(comp_codes)]
|
|
|
1228 names(colors) = comp_codes
|
|
|
1229 species = substring(V(GL$G)$name,1, nchar(comp_codes[1]))
|
|
|
1230 ## create image with highlighted domains
|
|
|
1231 png(fs$graph_comparative, width = PNGWIDTH + LEGENDWIDTH, height = PNGHEIGHT, pointsize = PS)
|
|
|
1232 layout(matrix(1:2, ncol = 2), width = c(10, 3))
|
|
|
1233 plotg(GL$G,GL$L, col = colors[species])
|
|
|
1234 par(mar = c(0, 0, 0, 0))
|
|
|
1235 plot.new()
|
|
|
1236 legend("topleft", col=colors,
|
|
|
1237 legend = names(colors),
|
|
|
1238 pch = 15, cex = 0.7)
|
|
|
1239 dev.off()
|
|
|
1240 HTML.title("comparative analysis:", HR=4, file = clinfo$html_report_main)
|
|
|
1241 html_insert_image(
|
|
|
1242 img_file = fs_relative$graph_comparative,
|
|
|
1243 htmlfile = clinfo$html_report_main)
|
|
|
1244 ## comparative analysis - calculate statistics:
|
|
|
1245 V1V2 = substring(get.edgelist(GL$G),1, nchar(comp_codes[1]))
|
|
|
1246 C_observed = table(data.frame(
|
|
|
1247 V1=factor(V1V2[,1],levels = comp_codes),
|
|
|
1248 V2=factor(V1V2[,2],levels = comp_codes)
|
|
|
1249 ))
|
|
|
1250 C_observed = as.data.frame.matrix(C_observed + t(C_observed))
|
|
|
1251 P_observed = as.data.frame.matrix(C_observed/sum(C_observed))
|
|
|
1252 Edge_propotions = table(factor(V1V2, levels = comp_codes))/(length(V1V2))
|
|
|
1253 P_expected = matrix(Edge_propotions,ncol=1) %*% matrix(Edge_propotions, nrow=1)
|
|
|
1254
|
|
|
1255 spec_counts = df2html(
|
|
|
1256 data.frame(table(factor(species, levels=comp_codes))),
|
|
|
1257 header = c("Species", "Read count")
|
|
|
1258 )
|
|
|
1259
|
|
|
1260 edge_count = df2html(
|
|
|
1261 data.frame(species = comp_codes, (C_observed/2)[,]),
|
|
|
1262 header = c("", comp_codes)
|
|
|
1263 )
|
|
|
1264
|
|
|
1265 P_OE = df2html(
|
|
|
1266 data.frame(species=comp_codes,(P_observed/P_expected[,])),
|
|
|
1267 header = c("", comp_codes)
|
|
|
1268 )
|
|
|
1269
|
|
|
1270 HTML.title("Comparative analysis - species read counts:", HR =5, file = clinfo$html_report_main)
|
|
|
1271 cat(spec_counts,
|
|
|
1272 file = clinfo$html_report_main, append =TRUE)
|
|
|
1273
|
|
|
1274 # show connection only if more than one species in cluster
|
|
|
1275 if (length(unique(species))>1){
|
|
|
1276 HTML.title("comparative analysis - number of edges between species:",
|
|
|
1277 HR =5, file = clinfo$html_report_main)
|
|
|
1278 cat(edge_count,
|
|
|
1279 file = clinfo$html_report_main, append = TRUE)
|
|
|
1280
|
|
|
1281 HTML.title("comparative analysis - observed/expected number of edges between species",
|
|
|
1282 HR =5, file = clinfo$html_report_main)
|
|
|
1283 cat(P_OE,
|
|
|
1284 file = clinfo$html_report_main, append = TRUE)
|
|
|
1285 }
|
|
|
1286
|
|
|
1287 }
|
|
|
1288
|
|
|
1289
|
|
|
1290 ## create image with highlighted domains
|
|
|
1291 png(fs$graph_domains, width = PNGWIDTH + LEGENDWIDTH, height = PNGHEIGHT, pointsize = PS)
|
|
|
1292 layout(matrix(1:2, ncol = 2), width = c(10, 3))
|
|
|
1293 plotg(GL$G,GL$L, col = vertex_colors$color_table)
|
|
|
1294 domains_detected = length(vertex_colors$legend) > 0
|
|
|
1295 par(mar = c(0, 0, 0, 0))
|
|
|
1296 plot.new()
|
|
|
1297 if (domains_detected){
|
|
|
1298 # domains found
|
|
|
1299 legend("topleft", col=vertex_colors$legend,
|
|
|
1300 legend = names(vertex_colors$legend),
|
|
|
1301 pch = 15, cex = 0.7)
|
|
|
1302 }
|
|
|
1303 dev.off()
|
|
|
1304
|
|
|
1305 HTML.title("protein domains:", HR=4, file = clinfo$html_report_main)
|
|
|
1306 if (!domains_detected){
|
|
|
1307 HTML("No protein domains detected", file = clinfo$html_report_main)
|
|
|
1308 }
|
|
|
1309 HTML("protein domains:", HR=4, file = clinfo$html_report_main)
|
|
|
1310 html_insert_image(
|
|
|
1311 img_file = fs_relative$graph_domains,
|
|
|
1312 htmlfile = clinfo$html_report_main)
|
|
|
1313
|
|
|
1314 #############################################################################
|
|
|
1315 if (nrow(annot_summary) == 0){
|
|
|
1316 HTML.title("Reads annotation summary", HR = 3, file = clinfo$html_report_main)
|
|
|
1317 HTML("No similarity hits to repeat databases found", file = clinfo$html_report_main)
|
|
|
1318 }else{
|
|
|
1319 HTML(annot_summary, file = clinfo$html_report_main, align = "left")
|
|
|
1320 }
|
|
|
1321
|
|
|
1322
|
|
|
1323 ## similarity and mate cluster
|
|
|
1324 mate_clusters = get_cluster_connection_info(index, search_type="pair")
|
|
|
1325 similar_clusters = get_cluster_connection_info(index, search_type="similarity")
|
|
|
1326 ## report mate and similarity clusters
|
|
|
1327 if (!is.null(similar_clusters)){
|
|
|
1328 HTML.title("clusters with similarity:", file =clinfo$html_report_main, HR = 3)
|
|
|
1329 cat(df2html(
|
|
|
1330 similar_clusters,
|
|
|
1331 header=c("Cluster","Number of similarity hits"),
|
|
|
1332 sort_col = "N", scroling = TRUE
|
|
|
1333 ),
|
|
|
1334 file =clinfo$html_report_main, append=TRUE)
|
|
|
1335 }
|
|
|
1336 if (!is.null(mate_clusters)){
|
|
|
1337 HTML.title("clusters connected through mates:", file =clinfo$html_report_main, HR = 3)
|
|
|
1338 cat(df2html(
|
|
|
1339 mate_clusters[,c('cl','N','k')],
|
|
|
1340 header = c('Cluster','Number of shared<br> read pairs','k'),
|
|
|
1341 sort_col = "N", scroling = TRUE
|
|
|
1342 ),
|
|
|
1343 file = clinfo$html_report_main,append = TRUE
|
|
|
1344 )
|
|
|
1345
|
|
|
1346 ## create base graph images - it will serve as background for
|
|
|
1347 ## mate clusters plots
|
|
|
1348 png(fs$graph_base,
|
|
|
1349 width = PNGWIDTH, height = PNGHEIGHT, pointsize = PS)
|
|
|
1350 par(mar=c(0,0,0,0),xaxs="i",yaxs="i")
|
|
|
1351 plotg(GL$G,GL$L, col = "#00000050")
|
|
|
1352 dev.off()
|
|
|
1353 ## load base as raster image
|
|
|
1354 base_image = readPNG(fs$graph_base)
|
|
|
1355
|
|
|
1356 for (i in order(mate_clusters$N, decreasing = TRUE)){
|
|
|
1357 mate_ids = unlist(strsplit(mate_clusters$ids[[i]],split=","))
|
|
|
1358 ## print only graph above MAX_N mates
|
|
|
1359 if (length(mate_ids) < MAX_N){ # TODO - use constant
|
|
|
1360 next
|
|
|
1361 }
|
|
|
1362 png(sprintf(fs$graph_mates,index, mate_clusters$cl[i]),
|
|
|
1363 width = PNGWIDTH, height = PNGHEIGHT, pointsize = PS)
|
|
|
1364 color_mate = gsub(".$","",V(GL$G)$name) %in% mate_ids %>%
|
|
|
1365 ifelse("#FF0000FF", "#000000AA")
|
|
|
1366 par(mar=c(0,0,0,0),xaxs="i",yaxs="i")
|
|
|
1367 plot(range(GL$L[,1]), range(GL$L[,2]), type = "n", xlab = "", ylab = "", axes = FALSE,
|
|
|
1368 main = paste0("CL",index," ----> CL",mate_clusters$cl[i]))
|
|
|
1369 rasterImage(base_image,
|
|
|
1370 range(GL$L[,1])[1], range(GL$L[,2])[1],
|
|
|
1371 range(GL$L[,1])[2], range(GL$L[,2])[2]
|
|
|
1372 )
|
|
|
1373 points(GL$L[,1:2], col = color_mate,pch=18,cex=.8)
|
|
|
1374 dev.off()
|
|
|
1375 title = paste0("CL",index," ----> CL",mate_clusters$cl[i])
|
|
|
1376 footer = paste0("No. of shared pairs: :", mate_clusters$N[i])
|
|
|
1377 html_insert_floating_image(
|
|
|
1378 img_file = sprintf(fs_relative$graph_mates, index, mate_clusters$cl[i]),
|
|
|
1379 htmlfile = clinfo$html_report_main, width = 200,
|
|
|
1380 title = title, footer = footer
|
|
|
1381 )
|
|
|
1382 }
|
|
|
1383 }
|
|
|
1384 }
|
