annotate plot_comparative_clustering_summary.R @ 29:53dc6aef5441 draft

planemo upload commit 20bdf879b52796d3fb251a20807191ff02084d3c-dirty
author petr-novak
date Thu, 03 Aug 2023 07:32:40 +0000
parents 5a05925340b0
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
17
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
1 #!/usr/bin/env Rscript
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
2 library(optparse)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
3 ## TODO - add scale to legend!
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
4 twenty_colors = c(
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
5 '#e6194b', '#3cb44b', '#ffe119', '#4363d8', '#f58231',
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
6 '#911eb4', '#46f0f0', '#f032e6', '#bcf60c', '#fabebe',
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
7 '#008080', '#e6beff', '#9a6324', '#fffac8', '#800000',
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
8 '#aaffc3', '#808000', '#ffd8b1', '#000075', "#000000"
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
9 )
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
10
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
11 get_color = function(classification, size){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
12 ## 20 of unique colors, first is black
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
13 unique_colors = twenty_colors[1:opt$number_of_colors]
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
14 Ncol = length(unique_colors)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
15 ## rest wil be grey:
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
16 grey_color = "#a9a9a9"
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
17 ## unique repeats without All
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
18 include = !classification %in% "All"
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
19 unique_repeats = names(c(sort(by(size[include], INDICES = classification[include], FUN = sum), decreasing = TRUE)))
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
20 color_table = unique_colors[1:min(Ncol,length(unique_repeats))]
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
21 names(color_table) = unique_repeats[1:min(Ncol,length(unique_repeats))]
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
22 color = rep(grey_color, length(classification))
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
23 names(color) = classification
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
24 for (ac in names(color_table)){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
25 color[names(color) %in% ac] = color_table[ac]
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
26 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
27 return(color)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
28 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
29
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
30
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
31 make_legend = function(color){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
32 ## simplify description:
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
33 names(color) = gsub(".+/","",names(color))
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
34 description = sapply(split(names(color), color), function(x) paste(unique(x), collapse=";"))
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
35 description = gsub(".+;.+", "Other", description)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
36 description = gsub("All", "Other", description)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
37 if ("Other" %in% description & length(description) > 1){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
38 description = c(description[! description %in% "Other"], description[description %in% "Other"])
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
39 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
40 ord = order(factor(names(description), levels = twenty_colors))
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
41 legend_info = list(name = gsub("All", "NA", description)[ord], color = names(description)[ord])
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
42 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
43
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
44 plot_rect_map = function(read_counts,cluster_annotation, output_file,GS, RL, Xcoef=1,Ycoef=1){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
45 ## read_counts : correspond to COMPARATIVE_ANALYSIS_COUNTS.csv
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
46 ## cluster annotation : CLUSTER_TABLE.csv
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
47 counts = read.table(read_counts,header=TRUE,as.is=TRUE)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
48 input_read_counts = unlist(read.table(read_counts, nrows = 1, comment.char = "",sep="\t")[-(1:2)])
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
49
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
50 counts_file_valid = ncol(counts) == (length(input_read_counts) + 2) & all(colnames(input_read_counts)[1:2]==c("cluster", "supercluster"))
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
51 ## find which line is header
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
52 header_line = grep(".*Cluster.*Supercluster.*Size", readLines(cluster_annotation))
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
53 annot = read.table(cluster_annotation, sep="\t",header=TRUE,as.is=TRUE, skip = header_line - 1)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
54 ## validate
20
5a05925340b0 Uploaded
petr-novak
parents: 17
diff changeset
55 annot_file_valid = all(c("Cluster","Supercluster","Size","Size_adjusted","Automatic_annotation","TAREAN_annotation","Final_annotation") %in% colnames(annot))
17
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
56
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
57
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
58 if (!annot_file_valid | !counts_file_valid){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
59 pdf(output_file)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
60 plot.new()
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
61 text(0.5,0.5,"Input is not valid, check input files!")
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
62 dev.off()
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
63 stop("Input files are not valid!")
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
64 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
65 print(annot_file_valid)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
66 print(counts_file_valid)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
67 ## remove counts which are not in annotation - only clusters in annot will be plotted!
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
68 counts = counts[annot$Cluster,]
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
69 N = nrow(annot)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
70
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
71 counts_automatic = counts
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
72 annot_automatic = annot
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
73 input_read_counts_automatic = input_read_counts
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
74 # remove organelar and contamination if required make count correction
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
75 if (opt$nuclear_only){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
76 exclude=grep("contamination|organelle",annot$Automatic_annotation)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
77 if (length(exclude)>0){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
78 counts_automatic = counts[-exclude, , drop=FALSE]
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
79 annot_automatic = annot[-exclude, ,drop=FALSE]
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
80 input_read_counts_automatic = input_read_counts - colSums(counts[exclude,-c(1:2) , drop=FALSE])
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
81 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
82 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
83 color_auto = get_color(annot_automatic$Automatic_annotation, annot_automatic$Size)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
84
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
85 legend_info = make_legend(color_auto)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
86 params = list(Automatic_annotation = list(
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
87 color = color_auto,
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
88 legend = legend_info,
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
89 counts = counts_automatic,
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
90 annot = annot_automatic,
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
91 input_read_counts = input_read_counts_automatic
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
92 )
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
93 )
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
94
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
95
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
96 if (!is.null(annot$Final_annotation)){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
97
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
98 ## column with manual annotation exist - check if correct
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
99 if (any(annot$Final_annotation %in% "" | any(is.na(annot$Final_annotation)))){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
100 message("Final annotation is not complete, skipping")
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
101 }else{
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
102
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
103 counts_manual = counts
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
104 annot_manual = annot
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
105 input_read_counts_manual = input_read_counts
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
106 ## correction must be done idependetly in case manual and automatic classification differ in annotation
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
107 if (opt$nuclear_only){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
108 exclude=grep("contamination|organelle",annot$Final_annotation)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
109 if (length(exclude)>0){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
110 counts_manual = counts[-exclude, , drop=FALSE]
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
111 annot_manual = annot[-exclude, ,drop=FALSE]
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
112 input_read_counts_manual = input_read_counts - colSums(counts[exclude,-c(1:2) , drop=FALSE])
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
113 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
114 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
115 ## append
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
116 color_manual = get_color(annot_manual$Final_annotation, annot_manual$Size)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
117 legend_info_manual = make_legend(color_manual)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
118
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
119 params$Final_annotation = list(
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
120 color = color_manual,
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
121 legend = legend_info_manual,
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
122 counts = counts_manual,
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
123 annot = annot_manual,
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
124 input_read_counts = input_read_counts_manual
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
125
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
126 )
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
127 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
128 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
129
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
130 ## set size of pdf output
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
131 wdth = (3 + N*0.03 ) * Xcoef
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
132 hgt = (2.2 + ncol(counts)*0.20) * Ycoef
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
133 if (!any(is.na(GS))){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
134 hgt = hgt + ncol(counts)*0.20 * Ycoef
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
135 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
136 ptsize = round((wdth*hgt)^(1/4))*5
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
137
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
138
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
139 pdf(output_file, width=wdth,height=hgt, pointsize = ptsize) # was 50
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
140 ## originaly - printing of both automatic and final annotation
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
141 ## now - print only final_annotation if available
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
142 if (length(params) == 2){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
143 ## remove automatic
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
144 params$Automatic_annotation = NULL
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
145 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
146 ##
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
147
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
148 for (j in seq_along(params)){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
149 Nclust = nrow(params[[j]]$annot)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
150 ##prepare matrix for plotting
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
151 M = as.matrix(params[[j]]$counts[1:Nclust,-(1:2)])
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
152 rownames(M) = paste0("CL",rownames(M))
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
153 Mn1=(M)/apply(M,1,max)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
154 Mn2=M/max(M)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
155 ord1 = hclust(dist(Mn1),method = "ward.D")$order
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
156 ord2 = hclust(dist(t(Mn2)))$order
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
157
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
158 ploting_area_width = 3 + log10(Nclust)*3
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
159 ploting_area_sides = 1.5
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
160 legend_width = 3
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
161 title_height = 0.5
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
162 if (any(is.na(GS))){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
163 layout(matrix(c(0,0,0,3,0,2,0,3,0,1,0,3),ncol=4,byrow = TRUE),
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
164 width=c(ploting_area_sides,ploting_area_width,ploting_area_sides, legend_width),
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
165 height=c(title_height, 3,ncol(M)*0.8))
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
166 }else{
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
167 ## extra row for legends
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
168
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
169
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
170 layout(matrix(c(0,0,0,3,0,2,0,3,0,1,0,3,0,0,0,4),ncol=4,byrow = TRUE),
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
171 width=c(ploting_area_sides,ploting_area_width,ploting_area_sides, legend_width),
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
172 height=c(title_height, 3,ncol(M)*0.8,ncol(M)*0.8 ))
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
173 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
174
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
175
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
176 par(xaxs='i', yaxs = 'i')
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
177 par(las=2,mar=c(4,0,0,0),cex.axis=0.5)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
178
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
179 if (any(is.na(GS))){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
180 rectMap(Mn2[ord1,ord2],scale.by='row',col=params[[j]]$color[ord1], grid=TRUE)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
181 }else{
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
182 # use genomic sizes
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
183 Mn3 = t(t(M) * (GS[colnames(M),] / params[[j]]$input_read_counts))[ord1,ord2]
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
184 ## rescale
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
185 MaxGS = max(Mn3)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
186 Mn3 = Mn3/max(Mn3)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
187 rectMap(Mn3,scale.by='none',col=params[[j]]$color[ord1], grid=TRUE)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
188 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
189 par(las=2,mar=c(1,0,1,0), mgp = c(2,1,0))
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
190 barplot(params[[j]]$annot$Size[ord1], col = 1)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
191 mtext(side = 2, "Cluster size", las = 3, line = 2.5, cex = 0.5)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
192 mtext(side=3, names(params)[j], las=0, line=1)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
193 plot.new()
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
194 legend("topleft", col= params[[j]]$legend$color, legend=params[[j]]$legend$name, pch=15, cex=0.7, bty="n", pt.cex=1)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
195 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
196
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
197 if (!any(is.na(GS))){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
198 ## plot GS scale
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
199 par(xaxs='i', yaxs = 'i')
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
200 print(log(nrow(Mn3)))
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
201 par(las=2,mar=c(4,0,0,log(nrow(Mn3))),cex.axis=0.5) # same par as recplot above to keep the scale
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
202 Mn3scale = Mn3
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
203 Mn3scale[,] = 0
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
204 colnames(Mn3scale)=rep("", ncol(Mn3scale))
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
205 rownames(Mn3scale)=rep("", nrow(Mn3scale))
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
206 Mn3scale[,1] = seq(0,1, length.out = nrow(Mn3))
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
207 rectMap(Mn3scale,scale.by='none',col="grey", grid=FALSE, boxlab="", draw_box=FALSE, center=FALSE)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
208 slabels = pretty(c(0,MaxGS), n = 10)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
209 sat = slabels/MaxGS * nrow(Mn3scale)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
210 axis(side=1, at= sat, labels = slabels, line = 0)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
211 mtext(side = 1, text = "Repeat abundance", las=1, line=2.5,cex=0.4)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
212 mtext(side = 2, text = "Rectangle\n height", las=1, line=2,cex=0.4, at=1)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
213
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
214 axis(2, at=c(0.5, 1, 1.5), labels=c(0,0.5,1),line=0)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
215 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
216 st = dev.off()
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
217 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
218
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
219 rectMap=function(x,scale.by='row',col=1,xlab="",ylab="",grid=TRUE,axis_pos=c(1,4),boxlab = "Cluster Id", cexx=NULL,cexy=NULL, draw_box=TRUE, center=TRUE){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
220 if (scale.by=='row'){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
221 #x=(x)/rowSums(x)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
222 x=(x)/apply(x,1,sum)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
223 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
224 if (scale.by=='column'){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
225 x=t(t(x)/apply(x,2,max))
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
226 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
227 nc=ncol(x)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
228 nr=nrow(x)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
229 coords=expand.grid(1:nr,1:nc)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
230 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)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
231 axis(axis_pos[1],at=1:nr,labels=rownames(x),lty=0,tick=FALSE,line=0,cex.axis=0.5/log10(nr))
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
232 axis(axis_pos[2],at=1:nc,labels=colnames(x),lty=0,tick=FALSE,las=2,line=0 ,hadj=0, cex.axis=0.7)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
233 axis(2,at=1:nc,labels=colnames(x),lty=0,tick=FALSE,las=2,line=0 ,hadj=1, cex.axis=0.7)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
234
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
235 mtext(side = 1, boxlab, las=1, line = 3, cex = 0.5)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
236 line = 1.5 + log10(nr)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
237 #mtext(side = 2, "Proportions of individual samples", las =0, line = line, cex = 0.5)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
238 s=x/2
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
239 w = c(x)/2
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
240 if(center){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
241 rect(coords[,1]-0.5,coords[,2]-s,coords[,1]+0.5,coords[,2]+s,col=col,border=NA)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
242 }else{
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
243 rect(coords[,1]-0.5,coords[,2]-0.5,coords[,1]+0.5,coords[,2]+x-0.5,col=col,border=NA)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
244 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
245 if (grid){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
246 abline(v=0:(nr)+.5,h=0:(nc)+.5,lty=2,col="#60606030",lwd=0.2)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
247 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
248 if(draw_box){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
249 box(col="#60606030",lty=2, lwd=0.2)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
250 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
251 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
252
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
253 option_list <- list(
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
254 make_option(c("-c", "--cluster_table"), default=NA, type = "character",
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
255 help="file from RepeatExplorer2 clustering - CLUSTER_TABLE.csv"),
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
256
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
257 make_option(c("-m", "--comparative_counts"),default = NA,type = "character",
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
258 help="file from RepeatExplorer2 output - COMPARATIVE_ANALYSIS_COUNTS.csv"),
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
259
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
260 make_option(c("-o", "--output"), type="character",
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
261 default="comparative_analysis_summary.pdf",
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
262 help="File name for output figures (pdf document)"),
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
263 make_option(c("-N", "--number_of_colors"), type="integer", default=10,
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
264 help="Number of unique colors used from plotting (2-20, default is 10)"),
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
265
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
266 make_option(c("-g", "--genome_size"),default = NA,type = "character",
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
267 help="file from genome sizes of species provided in tab delimited file in the format:
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
268
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
269 species_code1 GenomeSize1
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
270 species_code2 GenomeSize2
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
271 species_code3 GenomeSize3
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
272 species_code4 GenomeSize4
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
273
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
274 provide the same codes for species as in file COMPARATIVE_ANALYSIS_COUNTS.csv. The use of genome
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
275 sizes file imply the --nuclear_only option. If genome sizes are used, genomic abundance scale is added.
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
276 "),
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
277 make_option(c("-n", "--nuclear_only"), default = FALSE, type="logical",
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
278 action = "store_true",
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
279 help="remove all non-nuclear sequences (organelle and contamination). ")
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
280 )
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
281
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
282
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
283 opt = parse_args(OptionParser(option_list = option_list))
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
284
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
285 if (any(is.na(c(opt$cluster_table, opt$comparative_counts)))){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
286 message("\nBoth files: CLUSTER_TABLE.csv and COMPARATIVE_ANALYSIS_COUNTS.csv must be provided\n")
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
287 q()
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
288 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
289
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
290 if (!opt$number_of_colors %in% 1:20){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
291 message("number of color must be in range 1..20")
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
292 stop()
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
293 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
294
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
295 if (!is.na(opt$genome_size)){
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
296 GS = read.table(opt$genome_size, header=FALSE, as.is=TRUE, row.names = 1)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
297 opt$nuclear_only=TRUE
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
298 }else{
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
299 GS = NA
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
300 RL = NA
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
301 }
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
302
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
303 plot_rect_map(opt$comparative_counts, opt$cluster_table, opt$output, GS, RL)
d14b68e9fd1d Uploaded - new tools added
petr-novak
parents:
diff changeset
304