R embedded code in perl, ggsave failing to save pdf file

Question

I have a R code embedded in perl, the idea is to read a directory, put the files into an array, then loop through the array for each file, run the r rscript, that basically caluclates some stats, print them via write.table, and also make some plots, via ggplot and then save them via ggsave, everything works fine, at least the perl part, and the first part of the r script, it calculates the stats and write them correclty in the correct directory, however, when it comes to the plots, somehow ggsave doesn't work and i have no pdf file nowhere. what can the reason be? Before embedding in perl i was running the r script on r studio and ggsave was working fine. any help will be appreciated, i attach the code below:

my $mapfile = shift;
my $OTUdir = shift;
opendir my $OpOTUdir, $OTUdir;
my @OTUtab = grep {/\.txt/} readdir $OpOTUdir;
my $OTUtab;

foreach $OTUtab (@OTUtab){
my @splitname = split (/\_/, $OTUtab);
my $primer = $splitname[0];
my $dirout = "Subplot_dir_OTU".$primer.".out";
my $dir = "./".$primer."_plots/";
my $SubPlotdir = "bsub -o $dirout  mkdir $dir";
my $run1 = system ('bash','-c',"$bsub && $SubPlotdir ") == 0
 or die "Can't create the directories per primer set";
sleep 1 until -e "$dirout";

my $Rcmd_pure_OTU = "      #here starts the r script
library (\'vegan\')
raw_data <- read.csv(\"".$OTUdir."$OTUtab\", row.names =1, sep =\"\t\",   dec=\".\", header =T, skip =1)
pre_OTU_tab <- raw_data[,-which(names(raw_data) == \"taxonomy\")]
OTU_tab <- t(pre_OTU_tab)
log_OTU <- log10(OTU_tab + 1)

map <- read.csv(\"$mapfile\", sep =\"\", row.names =1, header = T)
ordered_Map <- map[match(row.names(OTU_tab), row.names(map)),]
re_ordered_Map <- ordered_Map[complete.cases(ordered_Map),]
dist_bray <- vegdist(log_OTU, method= \"bray\", binary=FALSE)
dist_bray_binary <- vegdist(log_OTU, method=\"bray\", binary=TRUE)
cap_bray <- capscale(dist_bray ~ 1)
cap_bray_bin <- capscale(dist_bray_binary ~ 1)
cap_bray_year <- capscale(dist_bray ~ as.factor(re_ordered_Map\$Origin))
cap_bray_bin_year <- capscale(dist_bray_binary ~     as.factor(re_ordered_Map\$Origin))
CA1perc_log <- (cap_bray\$CA\$eig[1]/sum(cap_bray\$CA\$eig))*100
CA1perc_log <- round(CA1perc_log, digits = 1)
CA2perc_log <- (cap_bray\$CA\$eig[2]/sum(cap_bray\$CA\$eig))*100
CA2perc_log <- round(CA2perc_log, digits = 1)
CA1perc_log_bin <- (cap_bray_bin\$CA\$eig[1]/sum(cap_bray_bin\$CA\$eig))*100
CA1perc_log_bin <- round(CA1perc_log_bin, digits = 1)
CA2perc_log_bin <- (cap_bray_bin\$CA\$eig[2]/sum(cap_bray_bin\$CA\$eig))*100
CA2perc_log_bin <- round(CA2perc_log_bin, digits = 1)
CAperc_log_year <- round((sum(cap_bray_year\$CCA\$eig)/sum(cap_bray_year\$CA\$eig, cap_bray_year\$CCA\$eig))*100, digits = 1)
CAperc_log_bin_year <- round((sum(cap_bray_bin_year\$CCA\$eig)/sum(cap_bray_bin_year\$CA\$eig, cap_bray_bin_year\$CCA\$eig))*100, digits = 1)
#print the stats 
#Year
SigTest1 <- anova(cap_bray_year, by=\"term\", step=9999, perm.max=9999)
CapResults1 <- c(CAperc_log_year, SigTest1\$'Pr(>F)'[1])
Results1 <- rbind(CapResults1) 
colnames(Results1) <- c(\"Percent_Correlated\", \"pval\")
write.table(Results1,      file=\"".$dir."Year_pvalue_constraints_".$primer.".txt\", sep=\"\t\", quote=F, col.names=NA)
#binary year
SigTest2 <- anova(cap_bray_bin_year, by=\"term\", step=9999, perm.max=9999)
CapResults2 <- c(CAperc_log_bin_year, SigTest2\$'Pr(>F)'[1])
Results2 <- rbind(CapResults2) 
colnames(Results2) <- c(\"Percent_Correlated\", \"pval\")
write.table(Results2,    file=\"".$dir."Binary_year_p_value_constraints_".$primer.".txt\", sep=\"\t\", quote=F, col.names=NA)

data_for_plot <- cbind(cap_bray\$CA\$u,re_ordered_Map)
data_for_plot_bin <- cbind(cap_bray_bin\$CA\$u,re_ordered_Map)

data_for_plot_year <- cbind(cap_bray_year\$CA\$u,re_ordered_Map)
data_for_plot_bin_year <- cbind(cap_bray_bin_year\$CA\$u,re_ordered_Map)
cbbPalette <- c(\"#000000\", \"#E69F00\", \"#56B4E9\", \"#009E73\", \"#F0E442\", \"#0072B2\", \"#D55E00\", \"#CC79A7\")
ColorCount <- length(unique(re_ordered_Map\$Origin))
GetPalette <- colorRampPalette(cbbPalette, bias =3, interpolate = \"spline\", alpha = TRUE)
plot_unco <- ggplot(data_for_plot) + 
geom_point( size=4, aes(x=MDS1, y=MDS2, shape= inf_uni, color = Origin),  position = position_jitter(w = 0.1, h = 0.1))+
geom_vline(xintercept = 0, size = 0.3) +
geom_hline(yintercept = 0, size = 0.3) +
scale_color_manual(values = GetPalette(ColorCount)) +
theme (
panel.background = element_blank(),
legend.key = element_rect (fill = \"white\"),
legend.text = element_text (size = 15),
legend.title = element_text (size = 17, face = \"bold\"),
axis.text = element_text(size =17),
axis.line = element_line(color= \"black\", size =0.6),
axis.title = element_text(size =19, face = \"bold\")    
)
ggsave(\"".$dir."My_".$primer."_Uncostrained.pdf\", plot = plot_unco, width=12, height=12, units = \"in\")
plot_bin_unco <- ggplot(data_for_plot_bin) + 
geom_point( size=4, aes(x=MDS1, y=MDS2, shape= inf_uni, color = Origin),   position = position_jitter(w = 0.1, h = 0.1))+
geom_vline(xintercept = 0, size = 0.3) +
geom_hline(yintercept = 0, size = 0.3) +
scale_color_manual(values = GetPalette(ColorCount)) +
theme (
panel.background = element_blank(),
legend.key = element_rect (fill = \"white\"),
legend.text = element_text (size = 15),
legend.title = element_text (size = 17, face = \"bold\"),
axis.text = element_text(size =17),
axis.line = element_line(color= \"black\", size =0.6),
axis.title = element_text(size =19, face = \"bold\")    
)
ggsave(\"".$dir."My_".$primer."_Uncostrained_bin.pdf\", plot =plot_bin_unco, width=12, height=12, units = \"in\")

sorry, follows the proper ending of the code, with closing brackets and semicolumn " ; and then $R-> send("$Rcmd_pure_OTU"); } — Alfredo Mari, Jul 20 '17 at 15:05
Why don't you just do this all in R? Have you actually looked at the R script that's produced? It's not easy to help you without a [reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) — MrFlick, Jul 20 '17 at 15:07
What packages do you have running? ggsave is part of `ggplot`, wondering what else may be conflicting. — Geochem B, Jul 20 '17 at 15:08
@MrFlick, well, i am not super good in r concerning file handling, i mostly use it for stats, whereas for directories and files i normally use perl, but any suggestion how to convert in r would be really appreciated, i am not sure what you mean by looking at the product of the r script, for the tables i have the right tables in the right directories, but not the pdf. Whereas the r script alone was running ok @ — Alfredo Mari, Jul 20 '17 at 19:50

score 0 · Answer 1 · answered Jul 24 '17 at 08:49

It is very difficult to debug a program that is stored as a string with interpolated variables. It would be better to move that R code into its own R file create_plots.R and then just pass in variables directly. Here are 2 approaches that would then work:

1 - Rewrite all of the file / directory operations in R

inside create_plots.R

# Find files:
files <- list.files(pattern = "\\.txt$")

# Create a directory:
dir.create(file.path(mainDir, subDir))
setwd(file.path(mainDir, subDir))

# Other stuff ...

2 - Invoke Rscript and pass it parameters from your perl script

inside create_plots.R

# Collect command line arguments
args  <- commandArgs( trailingOnly = TRUE )
foo   <- args[1]
bar   <- args[2]
baz   <- args[3]

inside create_plots.pl

# Invoke from perl script:
my $cmd = '/path/to/Rscript create_plots.R param1 param2 param3'
(system( $cmd ) == 0)
    or die "Unable to run '$cmd' : $!";

R embedded code in perl, ggsave failing to save pdf file

1 Answers1