2017-07-20 46 views
0

我有一個嵌入perl的R代碼,其思想是讀取一個目錄,將這些文件放到一個數組中,然後通過數組遍歷每個文件,運行r rscript,基本上解釋了一些統計信息,通過write.table打印出來,並通過ggplot製作一些圖,然後通過ggsave保存,一切正常,至少是perl部分和r腳本的第一部分,它會計算出統計數據並將它們寫在正確的目錄中,但是,當涉及到圖表時,不知何故,ggsave無法正常工作,我也沒有pdf格式的文件。 原因是什麼?在嵌入perl之前,我在r studio上運行r腳本,ggsave工作正常。 任何幫助將不勝感激,我附上下面的代碼:R在perl中嵌入的代碼,ggsave無法保存pdf文件

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

抱歉,下面的代碼的適當的結局,與閉合托架和semicolumn「;然後$ R-> send(「$ Rcmd_pure_OTU」); } –

+0

你爲什麼不在R中完成這一切?你有沒有看過製作的R腳本?如果沒有[可重現的例子],幫助你並不容易(https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) – MrFlick

+0

你有什麼軟件包在運行? ggsave是'ggplot'的一部分,想知道還有什麼可能是衝突的。 –

回答

0

這是非常困難的調試存儲與插值變量的字符串的程序。將R代碼移入它自己的R文件create_plots.R然後直接傳入變量會更好。這裏有兩個方法,然後將工作:

1 - 重寫所有的文件/目錄操作中的R

內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 - 調用RSCRIPT和從你的Perl腳本中傳遞參數

inside create_plot S.R

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

內部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' : $!";