Perl-R-Code

#!/usr/bin/perl -w

use strict;
use Getopt::Long;
use Statistics::R;

my ($temp_dir,$file_name,$file_type,$pval_col,$chart_title,$output_format,$plot_name);

usage() if (@ARGV < 1 or !GetOptions('tempdir=s' => \$temp_dir, 'filename=s' => \$file_name, 'filetype=s' => \$file_type, 'pvalcol=s' => \$pval_col, 'charttitle:s' => \$chart_title, 'outputformat=s' => \$output_format, 'plotname=s' => \$plot_name));

chdir($temp_dir);

my $R = Statistics::R->new();

$R -> run(qq(name="$chart_title"));

$R -> run(qq(ylabel = "Observed P-Vlaues (-log10)"));
$R -> run(qq(xlabel = "Theoretical P-Values (-log10)"));


if ($file_type eq "text/csv") {
        $R -> run(qq(traitfile <- read.csv("$file_name", header=TRUE, sep=",", as.is="$pval_col", na.strings="nan", stringsAsFactors=FALSE)));
} elsif ($file_type eq "text/plain") {
        $R -> run(qq(traitfile <- read.table("$file_name", header=TRUE, as.is="$pval_col", na.strings="nan", stringsAsFactors=FALSE)));
}

$R -> run(qq(ylimit = -log10(min(traitfile\$$pval_col ,na.rm=TRUE))+1));
$R -> run(qq(tlen <- length(traitfile\$$pval_col)));
$R -> run(qq(degf <- length(traitfile\$$pval_col) - 2));
$R -> run(qq(a <- rt(tlen, degf)));
$R -> run(qq(b <- pt(a, degf)));
$R -> run(qq(c <- -log10(traitfile\$$pval_col)));
$R -> run(qq(d <- -log10(b)));
$R -> run(qq(c[c > ylimit] <- ylimit));
$R -> run(qq(traitfile.fit <- lm(c ~ d)));

if ($output_format eq "png") {
        $R -> run(qq(png(filename="$plot_name", width=800, height=800, units="px", pointsize=12, bg="white", type="cairo")));
} elsif ($output_format eq "jpeg") {
        $R -> run(qq(jpeg(filename="$plot_name", width=800, height=800, units="px", pointsize=12, bg="white", type="cairo")));
} elsif ($output_format eq "pdf") {
        $R -> run(qq(pdf("$plot_name", width = 10, height = 10)));
}

$R -> run(qq(qqplot(d, c, main=name, ylab = "Observed P-Values (-log10)", xlab ="Theoretical P-Values (-log10)", ylim=c(0,ylimit), xlim=c(0,ylimit), pch=21, cex =0.5)));
$R -> run(qq(deno<-0.456));
$R -> run(qq(w1<-1.0-traitfile\$$pval_col));
$R -> run(qq(w2<-qchisq(w1,df=1)));
$R -> run(qq(lam1<-median(w2)/deno));
$R -> run(qq(text(0.1, ylimit, paste("GC Lambda= ", format(round(lam1,2), nsmall = 2)), pos=4, cex=1.5, col="darkgreen")));
$R -> run(qq(abline(a=0, b=1, col="blue")));
$R -> run(qq(dev.off()));

my $ret= $R -> read;
$R -> stopR();

exit;

sub usage {
        print "Unknown option: @_\n" if (@_);
        print "Usage: test_qqPlot.pl [options]
                      [--tempdir ----- full path of temp directory]
                      [--filename ----- input file name. this file contains data for Q-Q plot]
                      [--filetype ----- type of uploaded file]
                      [--pvalcol ------ column name in which p-values are stored]
                      [--charttitle ----- title for chart]
                      [--outputformat ----- format for plot]
                      [--plotname ----- name for plot]\n";
        exit;
}