Mercurial > repos > iuc > scater_create_qcmetric_ready_sce

annotate scater-plot-dist-scatter.R @ 2:b834074a9aff draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
author iuc
date Thu, 09 Sep 2021 12:23:11 +0000
parents fd808de478b1
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
1 #!/usr/bin/env Rscript
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
2
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
3 # Plot the distribution of read counts and feature counts, side by side, then a scatter plot of read counts vs feature counts below
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
4
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
5 # Load optparse we need to check inputs
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
6
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
7 library(optparse)
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
8 library(workflowscriptscommon)
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
9 library(LoomExperiment)
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
10 library(scater)
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
11 library(ggpubr)
1
fd808de478b1 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 61f3899168453092fd25691cf31871a3a350fd3b"
iuc
parents: 0
diff changeset
12 library(scales)
0
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
13
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
14 # parse options
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
15
2
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
16 option_list <- list(
0
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
17 make_option(
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
18 c("-i", "--input-loom"),
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
19 action = "store",
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
20 default = NA,
2
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
21 type = "character",
0
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
22 help = "A SingleCellExperiment object file in Loom format."
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
23 ),
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
24 make_option(
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
25 c("-o", "--output-plot-file"),
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
26 action = "store",
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
27 default = NA,
2
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
28 type = "character",
0
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
29 help = "Path of the PDF output file to save plot to."
1
fd808de478b1 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 61f3899168453092fd25691cf31871a3a350fd3b"
iuc
parents: 0
diff changeset
30 ),
fd808de478b1 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 61f3899168453092fd25691cf31871a3a350fd3b"
iuc
parents: 0
diff changeset
31 make_option(
fd808de478b1 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 61f3899168453092fd25691cf31871a3a350fd3b"
iuc
parents: 0
diff changeset
32 c("-l", "--log-scale"),
2
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
33 action = "store_true",
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
34 default = FALSE,
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
35 type = "logical",
1
fd808de478b1 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 61f3899168453092fd25691cf31871a3a350fd3b"
iuc
parents: 0
diff changeset
36 help = "Plot on log scale (recommended for large datasets)."
0
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
37 )
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
38 )
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
39
2
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
40 opt <- wsc_parse_args(option_list, mandatory = c("input_loom", "output_plot_file"))
0
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
41
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
42 # Check parameter values
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
43
2
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
44 if (! file.exists(opt$input_loom)) {
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
45 stop((paste("File", opt$input_loom, "does not exist")))
0
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
46 }
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
47
2
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
48 # Filter out unexpressed features
0
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
49
2
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
50 sce <- import(opt$input_loom, format = "loom", type = "SingleCellLoomExperiment")
0
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
51
2
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
52 # Do the scatter plot of reads vs genes
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
53 total_counts <- sce$total
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
54 total_features <- sce$detected
0
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
55 count_feats <- cbind(total_counts, total_features)
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
56 cf_dm <- as.data.frame(count_feats)
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
57
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
58 # Calculate binwidths for reads and features plots. Use 20 bins
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
59 read_bins <- max(total_counts / 1e6) / 20
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
60 feat_bins <- max(total_features) / 20
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
61
2
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
62 plot1 <- qplot(total_counts / 1e6, geom = "histogram", binwidth = read_bins, ylab = "Number of cells", xlab = "Read counts (millions)", fill = I("darkseagreen3")) + ggtitle("Read counts per cell")
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
63 plot2 <- qplot(total_features, geom = "histogram", binwidth = feat_bins, ylab = "Number of cells", xlab = "Feature counts", fill = I("darkseagreen3")) + ggtitle("Feature counts per cell")
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
64 plot3 <- ggplot(cf_dm, aes(x = total_counts / 1e6, y = total_features)) + geom_point(shape = 1) + geom_smooth() + xlab("Read count (millions)") +
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
65 ylab("Feature count") + ggtitle("Scatterplot of reads vs features")
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
66 plot4 <- plotColData(sce, y = "subsets_Mito_percent", x = "detected") + ggtitle("% MT genes") + geom_point(shape = 1) + theme(text = element_text(size = 15)) + theme(plot.title = element_text(size = 15)) + xlab("Total features") + ylab("% MT")
0
2d455a7e8a3d planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 5fdcafccb6c645d301db040dfeed693d7b6b4278
iuc
parents:
diff changeset
67
2
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
68 if (! opt$log_scale) {
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
69 final_plot <- ggarrange(plot1, plot2, plot3, plot4, ncol = 2, nrow = 2)
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
70 ggsave(opt$output_plot_file, final_plot, device = "pdf")
1
fd808de478b1 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 61f3899168453092fd25691cf31871a3a350fd3b"
iuc
parents: 0
diff changeset
71 } else {
2
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
72 plot1_log <- plot1 + scale_x_continuous(trans = "log10") + scale_y_continuous(trans = "log10")
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
73 plot2_log <- plot2 + scale_y_continuous(trans = "log10")
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
74 plot3_log <- plot3 + scale_y_continuous(trans = "log10")
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
75 plot4_log <- plot4 + scale_y_log10(labels = number)
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
76 final_plot_log <- ggarrange(plot1_log, plot2_log, plot3_log, plot4_log, ncol = 2, nrow = 2)
b834074a9aff "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 154318f74839a4481c7c68993c4fb745842c4cce"
iuc
parents: 1
diff changeset
77 ggsave(opt$output_plot_file, final_plot_log, device = "pdf")
1
fd808de478b1 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/scater commit 61f3899168453092fd25691cf31871a3a350fd3b"
iuc
parents: 0
diff changeset
78 }