Mercurial > repos > iuc > scater_filter

comparison scater-plot-dist-scatter.R @ 2:7a365ec81b52 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:24:17 +0000
parents b7ea9f09c02f
children
comparison
equal deleted inserted replaced
1:b7ea9f09c02f 2:7a365ec81b52
11 library(ggpubr) 11 library(ggpubr)
12 library(scales) 12 library(scales)
13 13
14 # parse options 14 # parse options
15 15
16 option_list = list( 16 option_list <- list(
17 make_option( 17 make_option(
18 c("-i", "--input-loom"), 18 c("-i", "--input-loom"),
19 action = "store", 19 action = "store",
20 default = NA, 20 default = NA,
21 type = 'character', 21 type = "character",
22 help = "A SingleCellExperiment object file in Loom format." 22 help = "A SingleCellExperiment object file in Loom format."
23 ), 23 ),
24 make_option( 24 make_option(
25 c("-o", "--output-plot-file"), 25 c("-o", "--output-plot-file"),
26 action = "store", 26 action = "store",
27 default = NA, 27 default = NA,
28 type = 'character', 28 type = "character",
29 help = "Path of the PDF output file to save plot to." 29 help = "Path of the PDF output file to save plot to."
30 ), 30 ),
31 make_option( 31 make_option(
32 c("-l", "--log-scale"), 32 c("-l", "--log-scale"),
33 action="store_true", 33 action = "store_true",
34 default=FALSE, 34 default = FALSE,
35 type = 'logical', 35 type = "logical",
36 help = "Plot on log scale (recommended for large datasets)." 36 help = "Plot on log scale (recommended for large datasets)."
37 ) 37 )
38 ) 38 )
39 39
40 opt <- wsc_parse_args(option_list, mandatory = c('input_loom', 'output_plot_file', 'log_scale')) 40 opt <- wsc_parse_args(option_list, mandatory = c("input_loom", "output_plot_file"))
41 41
42 # Check parameter values 42 # Check parameter values
43 43
44 if ( ! file.exists(opt$input_loom)){ 44 if (! file.exists(opt$input_loom)) {
45 stop((paste('File', opt$input_loom, 'does not exist'))) 45 stop((paste("File", opt$input_loom, "does not exist")))
46 } 46 }
47 47
48 # Input from Loom format 48 # Filter out unexpressed features
49 49
50 scle <- import(opt$input_loom, format='loom', type='SingleCellLoomExperiment') 50 sce <- import(opt$input_loom, format = "loom", type = "SingleCellLoomExperiment")
51 51
52 #do the scatter plot of reads vs genes 52 # Do the scatter plot of reads vs genes
53 total_counts <- scle$total_counts 53 total_counts <- sce$total
54 total_features <- scle$total_features_by_counts 54 total_features <- sce$detected
55 count_feats <- cbind(total_counts, total_features) 55 count_feats <- cbind(total_counts, total_features)
56 cf_dm <- as.data.frame(count_feats) 56 cf_dm <- as.data.frame(count_feats)
57 57
58 # Calculate binwidths for reads and features plots. Use 20 bins 58 # Calculate binwidths for reads and features plots. Use 20 bins
59 read_bins <- max(total_counts / 1e6) / 20 59 read_bins <- max(total_counts / 1e6) / 20
60 feat_bins <- max(total_features) / 20 60 feat_bins <- max(total_features) / 20
61 61
62 # Make the plots 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")
63 plot <- ggplot(cf_dm, aes(x=total_counts / 1e6, y=total_features)) + geom_point(shape=1) + geom_smooth() + xlab("Read count (millions)") + 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")
64 ylab("Feature count") + ggtitle("Scatterplot of reads vs features") 64 plot3 <- ggplot(cf_dm, aes(x = total_counts / 1e6, y = total_features)) + geom_point(shape = 1) + geom_smooth() + xlab("Read count (millions)") +
65 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") 65 ylab("Feature count") + ggtitle("Scatterplot of reads vs features")
66 plot2 <- qplot(total_features, geom="histogram", binwidth = feat_bins, ylab="Number of cells", xlab = "Feature counts", fill=I("darkseagreen3")) + ggtitle("Feature counts per cell") 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")
67 plot3 <- plotColData(scle, y="pct_counts_MT", x="total_features_by_counts") + ggtitle("% MT genes") + geom_point(shape=1) + theme(text = element_text(size=15)) + theme(plot.title = element_text(size=15))
68 67
69 if (! opt$log_scale){ 68 if (! opt$log_scale) {
70 final_plot <- ggarrange(plot1, plot2, plot, plot3, ncol=2, nrow=2) 69 final_plot <- ggarrange(plot1, plot2, plot3, plot4, ncol = 2, nrow = 2)
71 ggsave(opt$output_plot_file, final_plot, device="pdf") 70 ggsave(opt$output_plot_file, final_plot, device = "pdf")
72 } else { 71 } else {
73 plot_log_both <- plot + scale_x_continuous(trans = 'log10') + scale_y_continuous(trans = 'log10') 72 plot1_log <- plot1 + scale_x_continuous(trans = "log10") + scale_y_continuous(trans = "log10")
74 plot1_log <- plot1 + scale_y_continuous(trans = 'log10') 73 plot2_log <- plot2 + scale_y_continuous(trans = "log10")
75 plot2_log <- plot2 + scale_y_continuous(trans = 'log10') 74 plot3_log <- plot3 + scale_y_continuous(trans = "log10")
76 plot3_log <- plot3 + scale_y_log10(labels=number) 75 plot4_log <- plot4 + scale_y_log10(labels = number)
77 final_plot_log <- ggarrange(plot1_log, plot2_log, plot_log_both, plot3_log, ncol=2, nrow=2) 76 final_plot_log <- ggarrange(plot1_log, plot2_log, plot3_log, plot4_log, ncol = 2, nrow = 2)
78 ggsave(opt$output_plot_file, final_plot_log, device="pdf") 77 ggsave(opt$output_plot_file, final_plot_log, device = "pdf")
79 } 78 }