Mercurial > repos > iuc > scater_plot_dist_scatter

comparison scater-plot-dist-scatter.R @ 1:2e41b35b5bdd 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 61f3899168453092fd25691cf31871a3a350fd3b"
author iuc
date Tue, 03 Sep 2019 14:25:32 -0400
parents 4887c4c69847
children 81e5bdff4853
comparison
equal deleted inserted replaced
0:4887c4c69847 1:2e41b35b5bdd
7 library(optparse) 7 library(optparse)
8 library(workflowscriptscommon) 8 library(workflowscriptscommon)
9 library(LoomExperiment) 9 library(LoomExperiment)
10 library(scater) 10 library(scater)
11 library(ggpubr) 11 library(ggpubr)
12 library(scales)
12 13
13 # parse options 14 # parse options
14 15
15 option_list = list( 16 option_list = list(
16 make_option( 17 make_option(
24 c("-o", "--output-plot-file"), 25 c("-o", "--output-plot-file"),
25 action = "store", 26 action = "store",
26 default = NA, 27 default = NA,
27 type = 'character', 28 type = 'character',
28 help = "Path of the PDF output file to save plot to." 29 help = "Path of the PDF output file to save plot to."
30 ),
31 make_option(
32 c("-l", "--log-scale"),
33 action="store_true",
34 default=FALSE,
35 type = 'logical',
36 help = "Plot on log scale (recommended for large datasets)."
29 ) 37 )
30 ) 38 )
31 39
32 opt <- wsc_parse_args(option_list, mandatory = c('input_loom', 'output_plot_file')) 40 opt <- wsc_parse_args(option_list, mandatory = c('input_loom', 'output_plot_file', 'log_scale'))
33 41
34 # Check parameter values 42 # Check parameter values
35 43
36 if ( ! file.exists(opt$input_loom)){ 44 if ( ! file.exists(opt$input_loom)){
37 stop((paste('File', opt$input_loom, 'does not exist'))) 45 stop((paste('File', opt$input_loom, 'does not exist')))
49 57
50 # Calculate binwidths for reads and features plots. Use 20 bins 58 # Calculate binwidths for reads and features plots. Use 20 bins
51 read_bins <- max(total_counts / 1e6) / 20 59 read_bins <- max(total_counts / 1e6) / 20
52 feat_bins <- max(total_features) / 20 60 feat_bins <- max(total_features) / 20
53 61
54 #make the plots 62 # Make the plots
55 plot <- ggplot(cf_dm, aes(x=total_counts / 1e6, y=total_features)) + geom_point(shape=1) + geom_smooth() + xlab("Read count (millions)") + 63 plot <- ggplot(cf_dm, aes(x=total_counts / 1e6, y=total_features)) + geom_point(shape=1) + geom_smooth() + xlab("Read count (millions)") +
56 ylab("Feature count") + ggtitle("Scatterplot of reads vs features") 64 ylab("Feature count") + ggtitle("Scatterplot of reads vs features")
57 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 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")
58 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 plot2 <- qplot(total_features, geom="histogram", binwidth = feat_bins, ylab="Number of cells", xlab = "Feature counts", fill=I("darkseagreen3")) + ggtitle("Feature counts per cell")
59 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)) 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))
60 68
61 final_plot <- ggarrange(plot1, plot2, plot, plot3, ncol=2, nrow=2) 69 if (! opt$log_scale){
62 ggsave(opt$output_plot_file, final_plot, device="pdf") 70 final_plot <- ggarrange(plot1, plot2, plot, plot3, ncol=2, nrow=2)
71 ggsave(opt$output_plot_file, final_plot, device="pdf")
72 } else {
73 plot_log_both <- plot + scale_x_continuous(trans = 'log10') + scale_y_continuous(trans = 'log10')
74 plot1_log <- plot1 + scale_y_continuous(trans = 'log10')
75 plot2_log <- plot2 + scale_y_continuous(trans = 'log10')
76 plot3_log <- plot3 + scale_y_log10(labels=number)
77 final_plot_log <- ggarrange(plot1_log, plot2_log, plot_log_both, plot3_log, ncol=2, nrow=2)
78 ggsave(opt$output_plot_file, final_plot_log, device="pdf")
79 }