utils.R 13.6 KB
Newer Older
1
2
#!/usr/bin/Rscript

Valentina Galata's avatar
Valentina Galata committed
3
## IMPORT
4
5
suppressMessages(library(ggsci)) # colors

Valentina Galata's avatar
Valentina Galata committed
6
7
8
##############################
# INPUT
read_nanostats <- function(fname){
9
    print(sprintf("Reading: %s", fname))
Valentina Galata's avatar
Valentina Galata committed
10
11
12
13
14
15
16
17
18
19
20
21
    df <- read.csv(
        file=fname,
        sep="\t",
        header=TRUE,
        check.names=FALSE,
        stringsAsFactors=FALSE
    )
    df_cols <- c("stat"="Statistic", "value"="Value")
    colnames(df) <- df_cols[colnames(df)]
    return(df)
}

22
read_quast <- function(fname){
23
    print(sprintf("Reading: %s", fname))
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
    df <- read.csv(
        file=fname,
        sep="\t",
        header=TRUE,
        row.names=1,
        check.names=FALSE,
        stringsAsFactors=FALSE
    )
    testit::assert(all(colnames(df) %in% names(ASM_TOOL_NAMES)))
    colnames(df) <- ASM_TOOL_NAMES[colnames(df)]
    df <- df[QUAST_VARS, ASM_TOOL_NAMES]
    return(df)
}

read_plasflow <- function(fname){
39
    print(sprintf("Reading: %s", fname))
40
41
42
43
44
45
46
47
48
49
50
51
52
53
    df <- read.csv(
        file=fname,
        sep="\t",
        header=TRUE,
        stringsAsFactors=FALSE,
        check.names=FALSE
    )
    testit::assert(all(df$tool %in% names(ASM_TOOL_NAMES)))
    df$tool <- ASM_TOOL_NAMES[df$tool]
    df <- df[df$label %in% names(PLASFLOW_NAMES$labels),]
    df$label <- PLASFLOW_NAMES$labels[df$label]
    return(df)
}

54
read_rgi <- function(fname){
55
    print(sprintf("Reading: %s", fname))
56
57
58
59
60
61
62
63
64
65
66
67
68
    df <- read.csv(
        file=fname,
        sep="\t",
        header=TRUE,
        stringsAsFactors=FALSE,
        check.names=FALSE
    )
    colnames(df) <- sapply(colnames(df), function(x){ ifelse(x %in% names(ASM_TOOL_NAMES), ASM_TOOL_NAMES[x], x) })
    testit::assert(all(df$col %in% names(RGI_NAMES$col)))
    df$col <- RGI_NAMES$col[df$col]
    return(df)
}

69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
read_barrnap <- function(fname){
    print(sprintf("Reading: %s", fname))
    df <- read.csv(
        file=fname,
        sep="\t",
        header=TRUE,
        stringsAsFactors=FALSE,
        check.names=FALSE
    )
    testit::assert(all(df$tool %in% names(ASM_TOOL_NAMES)))
    df$tool <- ASM_TOOL_NAMES[df$tool]
    testit::assert(all(df$kingdom %in% names(BARRNAP_KINGDOM_NAMES)))
    df$kingdom <- BARRNAP_KINGDOM_NAMES[df$kingdom]
    return(df)
}

85
read_crispr <- function(fname){
86
    print(sprintf("Reading: %s", fname))
87
88
89
90
91
92
93
94
95
96
97
98
99
100
    df <- read.csv(
        file=fname,
        sep="\t",
        header=TRUE,
        stringsAsFactors=FALSE,
        check.names=FALSE
    )
    testit::assert(all(df$crispr_tool %in% names(CRISPR_TOOL_NAMES)))
    testit::assert(all(df$asm_tool    %in% names(ASM_TOOL_NAMES)))
    df$crispr_tool <- CRISPR_TOOL_NAMES[df$crispr_tool]
    df$asm_tool    <- ASM_TOOL_NAMES[df$asm_tool]
    return(df)
}

101
read_mappability <- function(fname){
102
    print(sprintf("Reading: %s", fname))
103
104
105
106
107
108
109
    df <- read.csv(
        file=fname,
        sep="\t",
        header=TRUE,
        check.names=FALSE,
        stringsAsFactors=FALSE
    )
110
111
    testit::assert(all(df$tool %in% names(ASM_TOOL_NAMES)))
    df$tool <- ASM_TOOL_NAMES[df$tool]
112
113
114
    return(df)
}

115
read_diamond <- function(fname){
116
    print(sprintf("Reading: %s", fname))
117
118
119
120
121
122
123
124
125
126
127
128
    df <- read.csv(
        file=fname,
        sep="\t",
        header=TRUE,
        check.names=FALSE,
        stringsAsFactors=FALSE
    )
    testit::assert(all(df$tool %in% names(ASM_TOOL_NAMES)))
    df$tool <- ASM_TOOL_NAMES[df$tool]
    return(df)
}

129
read_prodigal <- function(fname){
130
    print(sprintf("Reading: %s", fname))
131
132
133
134
135
136
137
138
139
    df <- read.csv(
        file=fname,
        sep="\t",
        header=TRUE,
        check.names=FALSE,
        stringsAsFactors=FALSE
    )
    testit::assert(all(df$tool %in% names(ASM_TOOL_NAMES)))
    df$tool <- ASM_TOOL_NAMES[df$tool]
140
141
142
143
144
    return(df)
}

read_prodigal_gcounts <- function(fname){
    df <- read_prodigal(fname)
145
    df$partial_pct <- 100 * df$partial / df$total
146
147
148
    return(df)
}

149
150
151
152
153
154
read_prodigal_glength <- function(fname){
    df <- read_prodigal(fname)
    return(df)
}


155
read_cdhit <- function(fname){
156
    print(sprintf("Reading: %s", fname))
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
    df <- read.csv(
        file=fname,
        sep="\t",
        header=TRUE,
        check.names=FALSE,
        stringsAsFactors=FALSE
    )
    df <- df[,c("tool1", "tool2", "unique")]
    testit::assert(all(df$tool1 %in% names(ASM_TOOL_NAMES)))
    df$tool1 <- ASM_TOOL_NAMES[df$tool1]
    testit::assert(all(df$tool2 %in% names(ASM_TOOL_NAMES)))
    df$tool2 <- ASM_TOOL_NAMES[df$tool2]
    return(df)
}

172
read_mash_dist_reads <- function(fname){
173
    proc_name <- function(x){
174
175
176
177
178
        rtype <- basename(dirname(x))
        mtype <- basename(dirname(dirname(x)))
        testit::assert(rtype %in% names(READ_TYPES))
        testit::assert(mtype %in% names(META_TYPES))
        return(sprintf("%s %s", META_TYPES[mtype], READ_TYPES[rtype]))
179
180
181
182
183
184
185
    }
    dm <- read.csv(file=fname, sep='\t', header=TRUE, row.names=1, check.names=FALSE)
    colnames(dm) <- sapply(colnames(dm), proc_name)
    rownames(dm) <- sapply(rownames(dm), proc_name)
    return(dm)
}

186
187
188
189
190
191
192
193
194
195
196
197
198
199
read_mash_dist_asm <- function(fname){
    proc_name <- function(x){
        x <- basename(dirname(x))
        testit::assert(x %in% names(ASM_TOOL_NAMES))
        return(ASM_TOOL_NAMES[x])
    }
    dm <- read.csv(file=fname, sep='\t', header=TRUE, row.names=1, check.names=FALSE)
    colnames(dm) <- sapply(colnames(dm), proc_name)
    rownames(dm) <- sapply(rownames(dm), proc_name)
    # testit::assert(all(colnames(dm) %in% names(ASM_TOOL_NAMES)))
    # colnames(dm) <- ASM_TOOL_NAMES[colnames(dm)]
    # rownames(dm) <- ASM_TOOL_NAMES[rownames(dm)]
    return(dm)
}
200

201
##############################
202
203
204
205
206
207
208
209
210
211
212
# PLOTS

plot_rgi_overlap <- function(df, ctype, col){
    df <- df[df$type == ctype & df$col == col,]
    df_list <- lapply(ASM_TOOL_NAMES, function(x){ df[df[,x] > 0,"label"] })
    names(df_list) <- ASM_TOOL_NAMES[names(df_list)]
    UpSetR::upset(
        data=UpSetR::fromList(df_list),
        # overlap order
        order.by="degree",
        decreasing=FALSE,
213
214
        # number of sets to plot
        nsets=length(ASM_TOOL_NAMES),
215
216
217
        # y-label title
        mainbar.y.label=sprintf("Intersection size (%s hits, %s)", ctype, col),
        # text size
218
        text.scale = c(1.2, 1.2, 1.2, 1.2, 1.2, 1.2)#,
219
        # colors
220
221
222
223
224
225
226
227
        # set.metadata=list(
        #     data=data.frame(
        #         sets=names(df_list),
        #         Tool=names(df_list),
        #         stringsAsFactors=FALSE
        #     ), 
        #     plots=list(list(type="matrix_rows", column="Tool", colors=ASM_TOOL_COLORS, alpha=0.7))
        # )
228
229
230
    )
}

231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
plot_crispr_overlap <- function(df, asm_tool){
    asm_sets <- lapply(
        CRISPR_TOOL_NAMES,
        function(x){ unique(unlist(df[df$asm_tool==asm_tool & df$crispr_tool == x,"seq_id"])) }
    )
    names(asm_sets) <- CRISPR_TOOL_NAMES[names(asm_sets)]
    print(asm_sets)
    UpSetR::upset(
        data=UpSetR::fromList(asm_sets),
        order.by="degree",
        decreasing=FALSE,
        mainbar.y.label=sprintf("Contig intersection size (%s)", asm_tool),
        sets.x.label="Contigs w/ CRISPR array(s)"
    )
}

Valentina Galata's avatar
Valentina Galata committed
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
# THEMES
mappability_theme <-
    theme_bw() +
    theme(
        # grid
        panel.grid=element_blank(),
        # strip
        strip.background=element_rect(fill="white"),
        strip.text=element_text(size=12, color="black"),
        # axes
        axis.title=element_text(size=12, color="black"),
        axis.text.y=element_text(size=9, color="black"),
        axis.text.x=element_text(size=9, color="black", angle=90, vjust=0.5, hjust=1)
    )

crispr_theme <-
    theme_bw() +
    theme(
        # legend
        legend.title=element_blank(),
        # grid
        panel.grid=element_blank(),
        # strip
        strip.background=element_rect(fill="white"),
        strip.text=element_text(size=12, color="black"),
        # axes
        axis.title=element_text(size=12, color="black"),
        axis.text.y=element_text(size=9, color="black"),
        axis.text.x=element_text(size=9, color="black", angle=90, vjust=0.5, hjust=1)
    )

plasflow_theme <-
    theme_bw() +
    theme(
        # legend
        legend.title=element_blank(),
        # grid
        panel.grid=element_blank(),
        # axes
        axis.title=element_text(size=12, color="black"),
        axis.text.y=element_text(size=9, color="black"),
        axis.text.x=element_text(size=9, color="black", angle=90, vjust=0.5, hjust=1)
    )

prodigal_theme <-
    theme_bw() +
    theme(
        # legend
        legend.title=element_blank(),
        # grid
        panel.grid=element_blank(),
        # axes
        axis.title=element_text(size=12, color="black"),
        axis.text.y=element_text(size=9, color="black"),
        axis.text.x=element_text(size=9, color="black", angle=90, vjust=0.5, hjust=1)
    )

diamond_theme2 <-
    theme_bw() +
    theme(
        # legend
        legend.title=element_blank(),
        # grid
        panel.grid=element_blank(),
        # strip
        strip.background=element_rect(fill="white"),
        strip.text=element_text(size=9, color="black"),
        # axes
        axis.title=element_text(size=12, color="black"),
        axis.text.y=element_text(size=9, color="black"),
        axis.text.x=element_text(size=9, color="black", angle=90, vjust=0.5, hjust=1)
    )

rgi_theme <- 
    theme_bw() +
    theme(
        # grid
        panel.grid=element_blank(),
        # strip
        strip.background=element_rect(fill="white"),
        strip.text=element_text(size=9, color="black"),
        # axes
        axis.title=element_text(size=12, color="black"),
        axis.text.y=element_text(size=9, color="black"),
        axis.text.x=element_text(size=9, color="black", angle=90, vjust=0.5, hjust=1)
    )

334
335
barrnap_theme <- rgi_theme

Valentina Galata's avatar
Valentina Galata committed
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
quast_theme <-
    theme_bw() +
    theme(
        # grid
        panel.grid=element_blank(),
        # strip
        strip.background=element_rect(fill="white"),
        strip.text=element_text(size=12, color="black"),
        # axes
        axis.title=element_text(size=12, color="black"),
        axis.text.y=element_text(size=9, color="black"),
        axis.text.x=element_text(size=9, color="black", angle=90, vjust=0.5, hjust=1)
    )

cdhit_theme <-
    theme_bw() +
    theme(
        # grid
        panel.grid=element_blank(),
        # strip
        strip.background=element_rect(fill="white"),
        strip.text=element_text(size=9, color="black"),
        # axes
        axis.title=element_text(size=12, color="black"),
        axis.text.y=element_text(size=9, color="black"),
        axis.text.x=element_text(size=9, color="black", angle=90, vjust=0.5, hjust=1)
    )

364
365
366
##############################
# CONST

367
368
369
370
371
372
373
374
375
376
377
###############
# Reads
META_TYPES <- c(
    "metag"="metaG",
    "metat"="metaT"
)
READ_TYPES <- c(
    "sr"="SR",
    "lr"="LR"
)

378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
###############
# Diamond
DIAMOND_VARS <- c(
    "qs_len_ratio"="Query/subject length ratio",
    "qcov"="Query coverage",
    "scov"="Subject coverage"
)
DIAMOND_VAR_LABELLER <- function(x){
    if(x %in% names(DIAMOND_VARS)){
        return(DIAMOND_VARS[x])
    } else {
        return(x)
    }
}

393
###############
394
# Assemblers
395
396
397

# names
ASM_TOOL_NAMES <- c(
398
    # LR, SR, hybrid
399
400
    "flye"="Flye",
    "megahit"="MEGAHIT",
401
    "metaspades"="metaSPAdes",
402
403
404
405
406
407
408
    "operams"="OPERA-MS",
    "metaspadeshybrid"="metaSPAdes (H)",
    # polishing w/ metaT (w/o LR)
    "megahitmetatracon"="Racon(MEGAHIT + metaT)",
    "metaspadesmetatracon"="Racon(metaSPAdes + metaT)",
    "operamsmetatracon"="Racon(OPERA-MS + metaT)",
    "metaspadeshybridmetatracon"="Racon(metaSPAdes (H) + metaT)"
409
)
Valentina Galata's avatar
Valentina Galata committed
410
411
412
413
414
ASM_TOOL_NAMES <- ASM_TOOL_NAMES[c(
    snakemake@config$assemblers$lr,
    snakemake@config$assemblers$sr,
    snakemake@config$assemblers$hy
)]
415
# colors
Valentina Galata's avatar
Valentina Galata committed
416
417
418
419
ASM_TOOL_PALETTE1 <- ggsci::pal_nejm("default", alpha=1)(5)
ASM_TOOL_PALETTE2 <- ggsci::pal_nejm("default", alpha=0.8)(5)
ASM_TOOL_PALETTE3 <- ggsci::pal_nejm("default", alpha=0.6)(5)
ASM_TOOL_PALETTE4 <- ggsci::pal_nejm("default", alpha=0.4)(5)
420
ASM_TOOL_COLORS <- c(
Valentina Galata's avatar
Valentina Galata committed
421
422
423
424
425
426
427
428
429
430
431
    # LR, SR, hybrid
    "Flye"=ASM_TOOL_PALETTE1[1],
    "MEGAHIT"=ASM_TOOL_PALETTE1[2],
    "metaSPAdes"=ASM_TOOL_PALETTE1[3],
    "OPERA-MS"=ASM_TOOL_PALETTE1[4],
    "metaSPAdes (H)"=ASM_TOOL_PALETTE1[5],
    # polishing w/ metaT (w/o LR)
    "Racon(MEGAHIT + metaT)"=ASM_TOOL_PALETTE2[2],
    "Racon(metaSPAdes + metaT)"=ASM_TOOL_PALETTE2[3],
    "Racon(OPERA-MS + metaT)"=ASM_TOOL_PALETTE2[4],
    "Racon(metaSPAdes (H) + metaT)"=ASM_TOOL_PALETTE2[5]
432
)
Valentina Galata's avatar
Valentina Galata committed
433
ASM_TOOL_COLORS <- ASM_TOOL_COLORS[ASM_TOOL_NAMES]
434

435
###############
436
# Gene tools
437
438
439
440
441
442
# GENE_TOOL_NAMES <- c(
#     "prodigal_partial"="Prodigal (partial)",
#     "prodigal_total"="Prodigal (total)",
#     "cdhit_unique"="CD-HIT (unique)",
#     "cdhit_total"="CD-HIT (total)"
# )
443

444
###############
445
# CRISPR tools
446

447
448
449
CRISPR_TOOL_NAMES <- c(
    "minced"="MinCED",
    "casc"="CasC"
450
451
)

452
###############
453
# PlasFlow
454
455

# names
456
457
458
459
460
461
PLASFLOW_NAMES <- list(
    statstype=c(
        count="Sequence count",
        sum="Cumulative sequence length [bp]",
        count_pct="Sequence count [%]",
        sum_pct="Cumulative sequence length [%]"
462
463
464
465
466
    ),
    labels=c(
        chromosome="Chromosome",
        plasmid="Plasmid",
        unclassified="Unclassified"
467
468
    )
)
469
# colors
470
PLASFLOW_COLORS <- list(
471
    labels=ggsci::pal_nejm("default", alpha=1)(4)[c(2,3,4)]
472
)
473
names(PLASFLOW_COLORS$labels) <- PLASFLOW_NAMES$labels
474

475
###############
476
477
# RGI
RGI_NAMES <- list(
478
    col=c(
479
480
481
482
483
484
        "Best_Hit_ARO"="ARO term",
        "ARO"="ARO",
        "Drug Class"="Drug class",       
        "Resistance Mechanism"="Resistance mechanism",
        "AMR Gene Family"="Gene family"
    )
485
486
)

487
488
489
490
491
492
493
494
495
###############
# barrnap
BARRNAP_KINGDOM_NAMES <- c(
    "arc"="Archea",
    "bac"="Bacteria",
    "euk"="Eukaryota",
    "mito"="Metazoan mitochondria"
)

496
###############
497
498
499
500
501
502
503
504
# QUAST
QUAST_VARS <- c(
    "# contigs",
    "Largest contig",            
    "Total length",
    "N50",
    "L50",                       
    "# N's per 100 kbp"
505
)