|
0
|
1 #!/usr/bin/env Rscript
|
|
|
2 suppressPackageStartupMessages(library(DBI))
|
|
|
3 suppressPackageStartupMessages(library(RSQLite))
|
|
|
4
|
|
|
5 CONNECTED = FALSE
|
|
|
6 if (FALSE) {
|
|
|
7 ## for testing
|
|
|
8 seqdb = "/mnt/raid/spolecny/petr/RE2/comparative_test/sequences.db"
|
|
|
9 hitsortdb = "/mnt/raid/spolecny/petr/RE2/comparative_test/hitsort.db"
|
|
|
10 class_file = "/mnt/raid/users/petr/workspace/repex_tarean/databases/classification_tree.rds"
|
|
|
11 ## connect to sqlite databases
|
|
|
12 SEQDB = dbConnect(RSQLite::SQLite(), seqdb)
|
|
|
13 HITSORTDB = dbConnect(RSQLite::SQLite(), hitsortdb)
|
|
|
14 CLS_TREE = readRDS(class_file)
|
|
|
15 }
|
|
|
16
|
|
|
17 connect_to_databases = function(seqdb, hitsortdb,classification_hierarchy_file = NULL){
|
|
|
18 if (!CONNECTED){
|
|
|
19 SEQDB <<- dbConnect(RSQLite::SQLite(), seqdb)
|
|
|
20 HITSORTDB <<- dbConnect(RSQLite::SQLite(), hitsortdb)
|
|
|
21 if (!is.null(classification_hierarchy_file)){
|
|
|
22 CLS_TREE <<- readRDS(classification_hierarchy_file)
|
|
|
23 }
|
|
|
24 CONNECTED <<- TRUE
|
|
|
25 }
|
|
|
26 }
|
|
|
27
|
|
|
28 disconnect_database = function(){
|
|
|
29 if (CONNECTED){
|
|
|
30 dbDisconnect(SEQDB)
|
|
|
31 dbDisconnect(HITSORTDB)
|
|
|
32 CONNECTED <<- FALSE
|
|
|
33 }
|
|
|
34 }
|
|
|
35
|
|
|
36 nested2named_list = function(x){
|
|
|
37 y = as.list(unlist(x[[1]]))
|
|
|
38 names(y) = unlist(x[[2]])
|
|
|
39 return(y)
|
|
|
40 }
|
|
|
41
|
|
|
42 is_comparative = function(){
|
|
|
43 prefix_codes = dbGetQuery(SEQDB,"SELECT * FROM prefix_codes")
|
|
|
44 if (nrow(prefix_codes) == 0){
|
|
|
45 return(FALSE)
|
|
|
46 }else{
|
|
|
47 return(TRUE)
|
|
|
48 }
|
|
|
49 }
|
|
|
50
|
|
|
51 get_comparative_codes = function(){
|
|
|
52 prefix_codes = dbGetQuery(SEQDB,"SELECT * FROM prefix_codes")
|
|
|
53 return(prefix_codes)
|
|
|
54 }
|
|
|
55
|
|
|
56 add_preamble = function(html_file, preamble){
|
|
|
57 html_content=readLines(html_file)
|
|
|
58 modified_html_content = gsub("<body>",
|
|
|
59 paste("<body>\n", preamble,"\n"),
|
|
|
60 html_content)
|
|
|
61 cat(modified_html_content, file = html_file, sep="\n")
|
|
|
62 }
|
|
|
63
|
|
|
64
|
|
|
65 df2html = function(df, header = NULL, sort_col = NULL, digits = 3, rounding_function=signif, decreasing = TRUE, scroling = FALSE, width = 300){
|
|
|
66 if (!is.null(sort_col)){
|
|
|
67 df = df[order(df[,sort_col], decreasing = decreasing),]
|
|
|
68 }
|
|
|
69 if (!is.null(digits)){
|
|
|
70 for (i in seq_along(df)){
|
|
|
71 if(is.numeric(df[,i])){
|
|
|
72 df[,i] = rounding_function(df[,i], digits)
|
|
|
73 }
|
|
|
74 }
|
|
|
75 }
|
|
|
76 if (is.null(header)){
|
|
|
77 h = ""
|
|
|
78 }else{
|
|
|
79 h = paste0(" <th>",header,"</th>\n", collapse="") %>%
|
|
|
80 paste0(" <tr>\n", .," </tr>\n")
|
|
|
81 }
|
|
|
82 x = apply(df,1,function(x)paste0(" <td>",x,"</td>\n", collapse="")) %>%
|
|
|
83 paste0(" <tr>\n", .," </tr>\n", collapse = "")
|
|
|
84 if (scroling){
|
|
|
85 cols = paste0('<col width="',rep(round(100/ncol(df)),ncol(df)),'%">\n',collapse ="")
|
|
|
86 height = min(200, 22 * nrow(df))
|
|
|
87 out = paste0(
|
|
|
88 '<table cellspacing="0" cellpadding="0" border="0" width="',width,'">\n',
|
|
|
89 ' <tr>\n',
|
|
|
90 ' <td>\n',
|
|
|
91 ' <table cellspacing="0" cellpadding="1" border="1" width="', width,'" >\n',
|
|
|
92 cols,
|
|
|
93 h,
|
|
|
94 ' </table>\n',
|
|
|
95 ' </td>\n',
|
|
|
96 ' </tr>\n',
|
|
|
97 ' <tr>\n',
|
|
|
98 ' <td>\n',
|
|
|
99 ' <div style="width:',width,'px; height:',height,'px; overflow:auto;">\n',
|
|
|
100 ' <table cellspacing="0" cellpadding="1" border="1" width="',width,'" >\n',
|
|
|
101 cols,
|
|
|
102 x,
|
|
|
103 ' </table>\n',
|
|
|
104 ' </div>\n',
|
|
|
105 ' </td>\n',
|
|
|
106 ' </tr>\n',
|
|
|
107 '</table>\n'
|
|
|
108 )
|
|
|
109
|
|
|
110 }else{
|
|
|
111 out = paste ("<table>\n", h,x, "</table>\n")
|
|
|
112 }
|
|
|
113 return(out)
|
|
|
114 }
|
|
|
115
|
|
|
116 start_html = function(filename, header){
|
|
|
117 cat(header, file = filename)
|
|
|
118 html_writer = function(content, fn=HTML, ...){
|
|
|
119 fn(content, append = TRUE, file = filename, ...)
|
|
|
120 }
|
|
|
121 }
|
|
|
122
|
|
|
123 preformatted = function(x){
|
|
|
124 ## make preformatted html text
|
|
|
125 return(
|
|
|
126 paste(
|
|
|
127 "<pre>\n",
|
|
|
128 x,
|
|
|
129 "</pre>"
|
|
|
130 ,sep="")
|
|
|
131 )
|
|
|
132 }
|
|
|
133
|
|
|
134
|
|
|
135 summary_histogram = function(fn, N_clustering, N_omit=0, size_adjusted=NULL, top_clusters_prop){
|
|
|
136 ## assume connection do databases
|
|
|
137 communities = dbGetQuery(
|
|
|
138 HITSORTDB,
|
|
|
139 "SELECT DISTINCT cluster, size, supercluster, supercluster_size FROM communities ORDER BY supercluster"
|
|
|
140 )
|
|
|
141 if (N_omit != 0){
|
|
|
142 ## adjust communities and cluster sizes:
|
|
|
143 cluster_to_adjust = which(
|
|
|
144 communities$size[order(communities$cluster)][1:length(size_adjusted)] != size_adjusted
|
|
|
145 )
|
|
|
146 ## keep original value:
|
|
|
147 communities$size_original = communities$size
|
|
|
148 superclusters_to_adjust = unique(communities$supercluster[communities$cluster %in% cluster_to_adjust])
|
|
|
149 for (cl in cluster_to_adjust){
|
|
|
150 communities[communities$cluster == cl,'size'] = size_adjusted[cl]
|
|
|
151 }
|
|
|
152 for (cl in superclusters_to_adjust){
|
|
|
153 communities[communities$supercluster == cl,'supercluster_size'] =
|
|
|
154 sum(communities[communities$supercluster == cl,'size'])
|
|
|
155 }
|
|
|
156 }else{
|
|
|
157 cluster_to_adjust=NULL
|
|
|
158 }
|
|
|
159 singlets = N_clustering - sum(communities$size)
|
|
|
160
|
|
|
161 supercluster_size = sort(unique(communities[, c('supercluster', 'supercluster_size')])$supercluster_size, decreasing = TRUE)
|
|
|
162
|
|
|
163
|
|
|
164 clid2size = sort(communities$size, decreasing = TRUE)
|
|
|
165
|
|
|
166 cluster_id = split(communities$cluster, communities$supercluster)
|
|
|
167 cluster_id_sort = lapply(cluster_id, function(x)x[order(clid2size[x], decreasing = FALSE)])
|
|
|
168
|
|
|
169 cluster_size_unsorted = split(communities$size, communities$supercluster)
|
|
|
170 cluster_size_sort = lapply(cluster_size_unsorted, function(x) (sort(x)))
|
|
|
171 ## reorder by size of superclusters
|
|
|
172 cluster_size_sort_sort = cluster_size_sort[order(sapply(cluster_size_sort, sum), decreasing = TRUE)]
|
|
|
173 cluster_id_sort_sort = cluster_id_sort[order(sapply(cluster_size_sort, sum), decreasing = TRUE)]
|
|
|
174
|
|
|
175
|
|
|
176 Nmax = max(sapply(cluster_size_sort_sort, length))
|
|
|
177 M = cbind(
|
|
|
178 sapply(cluster_size_sort_sort, function(x)y = c(x, rep(0, 1 + Nmax - length(x)))),
|
|
|
179 c(1, rep(0, Nmax))
|
|
|
180 )
|
|
|
181
|
|
|
182 Mid = cbind(
|
|
|
183 sapply(cluster_id_sort_sort, function(x)y = c(x, rep(0, 1 + Nmax - length(x)))),
|
|
|
184 c(1, rep(0, Nmax))
|
|
|
185 )
|
|
|
186
|
|
|
187 recolor = matrix(ifelse(Mid %in% cluster_to_adjust,TRUE,FALSE), ncol=ncol(Mid))
|
|
|
188 indices = which(recolor, arr.ind = TRUE)
|
|
|
189
|
|
|
190
|
|
|
191 png(fn, width = 1200, height = 700, pointsize = 20)
|
|
|
192 barplot(M,
|
|
|
193 width = c(supercluster_size, singlets),
|
|
|
194 space = 0, ylim = c(0, max(supercluster_size) * 1.2),
|
|
|
195 col = c("#CCCCCC"), names.arg = rep("", ncol(M)))
|
|
|
196 plot(0,
|
|
|
197 xlim = c(0, sum(c(supercluster_size, singlets))),
|
|
|
198 ylim = c(0, max(supercluster_size) * 1.2),
|
|
|
199 type = "n", yaxs = 'i', axes = FALSE, xlab = "", ylab = "")
|
|
|
200
|
|
|
201 rect(0, 0,
|
|
|
202 sum(supercluster_size),
|
|
|
203 max(supercluster_size) * 1.2,
|
|
|
204 col = "#0000FF10")
|
|
|
205
|
|
|
206 rect(sum(supercluster_size), 0,
|
|
|
207 sum(supercluster_size) + singlets,
|
|
|
208 max(supercluster_size) * 1.2,
|
|
|
209 col = "#FFAAFF10")
|
|
|
210
|
|
|
211 barplot(M,
|
|
|
212 width = c(supercluster_size, singlets),
|
|
|
213 space = 0, ylim = c(0, max(supercluster_size) * 1.2),
|
|
|
214 col = "#AAAAAA", names.arg = rep("", ncol(M)), add = TRUE,
|
|
|
215 xlab = "Proportion of reads [%]", ylab = "Number of reads",
|
|
|
216 main = paste(N_clustering, "reads total"))
|
|
|
217
|
|
|
218
|
|
|
219 for (i in seq_along(indices[,1])){
|
|
|
220 y1 = sum(M[1:indices[i,'row'],indices[i,'col']])
|
|
|
221 x1 = sum(M[,1:indices[i,'col']])
|
|
|
222 if(indices[i,'row'] == 1){
|
|
|
223 y0=0
|
|
|
224 }else{
|
|
|
225 y0 = sum(M[1:(indices[i,'row']-1),indices[i,'col']])
|
|
|
226 }
|
|
|
227 if (indices[i,'col']==1){
|
|
|
228 x0=0
|
|
|
229 }else{
|
|
|
230 x0 = sum(M[,1:(indices[i,'col']-1)])
|
|
|
231 }
|
|
|
232 rect(x0,y0,x1,y1, col="#88FF88")
|
|
|
233 }
|
|
|
234 abline(v=top_clusters_prop * sum(c(supercluster_size, singlets)), col="#00000088", lwd=3, lty=3)
|
|
|
235
|
|
|
236 text(sum(supercluster_size) / 2,
|
|
|
237 max(supercluster_size) * 1.05,
|
|
|
238 labels = paste0(sum(supercluster_size), " reads in\n",
|
|
|
239 length(supercluster_size), " supeclusters (", nrow(communities), " clusters)")
|
|
|
240 )
|
|
|
241
|
|
|
242
|
|
|
243 text(sum(supercluster_size) + singlets / 2,
|
|
|
244 max(supercluster_size) * 1.05,
|
|
|
245 labels = paste(singlets, "singlets"))
|
|
|
246
|
|
|
247 axis(1,at=seq(0,N_clustering,length.out=11),label=seq(0,100,by=10))
|
|
|
248 dev.off()
|
|
|
249 clustering_info = list(
|
|
|
250 Number_of_reads_in_clusters = sum(supercluster_size),
|
|
|
251 Number_of_clusters = nrow(communities),
|
|
|
252 Number_of_superclusters = length(supercluster_size),
|
|
|
253 Number_of_singlets = singlets
|
|
|
254 )
|
|
|
255 return(clustering_info)
|
|
|
256 }
|
|
|
257
|
|
|
258
|
|
|
259 rectMap=function(x,scale.by='row',col=1,xlab="",ylab="",grid=TRUE,axis_pos=c(1,4),cexx=NULL,cexy=NULL){
|
|
|
260 if (scale.by=='row'){
|
|
|
261 #x=(x)/rowSums(x)
|
|
|
262 x=(x)/apply(x,1,max)
|
|
|
263 }
|
|
|
264 if (scale.by=='column'){
|
|
|
265 x=t(t(x)/apply(x,2,max))
|
|
|
266 }
|
|
|
267 nc=ncol(x)
|
|
|
268 nr=nrow(x)
|
|
|
269 coords=expand.grid(1:nr,1:nc)
|
|
|
270 plot(coords[,1],coords[,2],type='n',axes=F,xlim=range(coords[,1])+c(-.5,.5),ylim=range(coords[,2])+c(-.5,.5),xlab=xlab,ylab=ylab)
|
|
|
271 axis(axis_pos[1],at=1:nr,labels=rownames(x),lty=0,tick=FALSE,line=0,cex.axis=0.5/log10(nr))
|
|
|
272 axis(axis_pos[2],at=1:nc,labels=colnames(x),lty=0,tick=FALSE,las=2,line=0 ,hadj=0, cex.axis=0.7)
|
|
|
273 axis(2,at=1:nc,labels=colnames(x),lty=0,tick=FALSE,las=2,line=0 ,hadj=1, cex.axis=0.7)
|
|
|
274
|
|
|
275 mtext(side = 1, "Cluster id", las=1, line = 3, cex = 0.5)
|
|
|
276 line = 1.5 + log10(nr)
|
|
|
277 mtext(side = 2, "Proportions of individual samples", las =0, line = line, cex = 0.5)
|
|
|
278 s=c(x)/2 # to get it proportional
|
|
|
279 w = c(x)/2
|
|
|
280 rect(coords[,1]-0.5,coords[,2]-s,coords[,1]+0.5,coords[,2]+s,col=col,border=NA)
|
|
|
281 if (grid){
|
|
|
282 abline(v=0:(nr)+.5,h=0:(nc)+.5,lty=2,col="#60606030")
|
|
|
283 }
|
|
|
284 box(col="#60606030",lty=2)
|
|
|
285 }
|
|
|
286
|
|
|
287 plot_rect_map = function(read_counts,cluster_annotation, output_file,Xcoef=1,Ycoef=1){
|
|
|
288 counts = read.table(read_counts,header=TRUE,as.is=TRUE)
|
|
|
289 annot = read.table(cluster_annotation, sep="\t",header=FALSE,as.is=TRUE)
|
|
|
290 N = nrow(annot)
|
|
|
291 colnames(annot) = c("cluster", "Automatic.classification")
|
|
|
292 annot$number.of.reads = rowSums(counts[1 : nrow(annot) ,-1])
|
|
|
293 unique_repeats = names(sort(table(c(annot$Automatic.classification,rep('nd',N))),decreasing = TRUE))
|
|
|
294
|
|
|
295 M = as.matrix(counts[1:N,-(1:2)])
|
|
|
296 rownames(M) = paste0("CL",rownames(M))
|
|
|
297 Mn1=(M)/apply(M,1,max)
|
|
|
298 Mn2=M/max(M)
|
|
|
299 Mn2=M/apply(M,1,sum)
|
|
|
300
|
|
|
301 ord1 = hclust(dist(Mn1),method = "ward.D")$order
|
|
|
302 ord2 = hclust(dist(t(Mn2)))$order
|
|
|
303 wdth = (400 + N*10 ) * Xcoef
|
|
|
304 hgt = (600 + ncol(M)*50) * Ycoef
|
|
|
305 ptsize = round((wdth*hgt)^(1/4))
|
|
|
306 png(output_file, width=wdth,height=hgt, pointsize = ptsize) # was 50
|
|
|
307 ploting_area_width = 3 + log10(N)*3
|
|
|
308 ploting_area_sides = 1
|
|
|
309 layout(matrix(c(4,2,3,4,1,3),ncol=3,byrow = TRUE),
|
|
|
310 width=c(ploting_area_sides,ploting_area_width,ploting_area_sides),
|
|
|
311 height=c(3,ncol(M)*0.5))
|
|
|
312 par(xaxs='i', yaxs = 'i')
|
|
|
313 par(las=2,mar=c(4,0,0,0),cex.axis=0.5)
|
|
|
314 rectMap(Mn2[ord1,ord2],scale.by='none',col=1, grid=TRUE)
|
|
|
315 par(las=2,mar=c(1,0,1,0), mgp = c(2,0.5,0))
|
|
|
316 barplot(annot$number.of.reads[ord1], col = 1)
|
|
|
317 mtext(side = 2, "Cluster size", las = 3, line = 2, cex = 0.5)
|
|
|
318 par(mar=c(0,0,10,0))
|
|
|
319 plot.new()
|
|
|
320 st = dev.off()
|
|
|
321 ## calculate coordinated if boxes to create hyperlink
|
|
|
322 X0 = wdth/(ploting_area_sides * 2 + ploting_area_width)* ploting_area_sides
|
|
|
323 X1 = wdth/(ploting_area_sides * 2 + ploting_area_width)*(ploting_area_sides + ploting_area_width)
|
|
|
324 L = round(seq(X0,X1, length.out = N + 1)[1:N])
|
|
|
325 R = round(seq(X0,X1, length.out = N + 1)[2:(N + 1)])
|
|
|
326 cn = rownames(Mn2[ord1,ord2])
|
|
|
327 cluster_links = paste0(
|
|
|
328 "seqclust/clustering/clusters/dir_CL",
|
|
|
329 sprintf("%04d", as.integer(substring(cn,3 ))),
|
|
|
330 "/index.html")
|
|
|
331 coords = paste0(L, ",", 1, ",", R, ",", hgt)
|
|
|
332 clustermap = paste0(
|
|
|
333 '\n<map name="clustermap"> \n',
|
|
|
334 paste0(
|
|
|
335 '<area shape="rect"\n coords="',coords, '"\n',
|
|
|
336 ' href="', cluster_links, '"\n',
|
|
|
337 ' title="', cn, '"/>\n',
|
|
|
338 collapse = ""),
|
|
|
339 "</map>\n")
|
|
|
340 return(clustermap)
|
|
|
341 }
|
