2017-07-20 4 views
0

Ich habe einen R-Code in Perl eingebettet, ist die Idee, ein Verzeichnis zu lesen, legen Sie die Dateien in ein Array, dann durchlaufen das Array für jede Datei, führen Sie die r rscript, das im Grunde einige Statistiken berechnet, sie über write.table druckt und auch einige Plots erstellt, über ggplot und dann über ggsave speichert, funktioniert alles gut, zumindest der Perl-Teil und der erste Teil des r-Skripts, es berechnet die stats und schreibt sie korrekt in das richtige Verzeichnis, aber wenn es um die Plots geht, irgendwie funktioniert ggsave nicht und ich habe nirgends eine pdf-Datei. Was kann der Grund sein? Vor dem Einbetten in Perl lief ich das r-Skript auf r Studio und Ggsave funktionierte gut. jede Hilfe zu erkennen ist, befestigen ich den Code unten:R eingebetteten Code in Perl, Ggsave nicht speichern PDF-Datei

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\") 
+0

sorry, folgt das richtige Ende des Codes, mit Schließbügel und semicolumn "; und dann $ R-> send ("$ Rcmd_pure_OTU"); } –

+0

Warum machst du das nicht alles in R? Hast du dir das R-Skript angeschaut? Es ist nicht einfach, Ihnen ohne ein [reproduzierbares Beispiel] zu helfen. (Https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) – MrFlick

+0

Welche Pakete haben Sie ausgeführt? ggsave ist Teil von 'ggplot' und fragt sich, was sonst noch widersprüchlich sein könnte. –

Antwort

0

Es ist sehr schwierig, ein Programm zu debuggen, die als String mit interpolierten Variablen gespeichert wird. Es wäre besser, diesen R-Code in seine eigene R-Datei create_plots.R zu verschieben und dann Variablen direkt weiterzuleiten. Hier sind zwei Ansätze, die dann funktionieren würde:

1 - alle Datei/Verzeichnisoperationen Rewrite in R

innen 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 und übergeben Sie es Parameter aus Ihrem Perl-Skript

innerhalb create_plot S. R

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

innen 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' : $!"; 
Verwandte Themen