source: trunk/brs2010p35/src/jag/plot_with_r.py @ 2

Last change on this file since 2 was 2, checked in by tim.te.beek@…, 8 years ago

Add revision 338 snapshot of previous SVN instance, omitting PLINK binaries

File size: 4.4 KB
Line 
1'''
2Created on Apr 12, 2011
3
4@author: maartenk, esther lips
5'''
6import os
7import subprocess
8import tempfile
9import logging as log
10
11class call_r(object):
12    '''
13    Class to group methods to call R.
14    '''
15
16    def __init__(self, inout):
17        '''Initialize Object to call R, with argument inout to
18        direct produced files to correct output folder.'''
19        #IO handle
20        self.inout = inout
21
22        #Find R binary using 'which'
23        self.r_binary = None
24        process = subprocess.Popen("which R", shell=True, stdout=subprocess.PIPE)
25        exit_status_which = process.wait()
26        if exit_status_which == 0:
27            self.r_binary = process.communicate()[0].strip()
28        else:
29            log.info ("R not found" + str(exit_status_which))
30
31
32    def execute(self, r_script, var_dict, outfile=None):
33        '''
34        Executing function
35        '''
36        #Skip if we could not find an R binary
37        if not self.r_binary:
38            return
39
40        args_text = ""
41        if not outfile:
42            outfile = tempfile.mkstemp()[1]
43
44        for key, val in var_dict.iteritems():
45            args_text += str(key) + "=\"" + str(val) + "\" "
46        command = self.r_binary + " CMD BATCH --no-save --no-restore '--args "\
47         + args_text + "' " + r_script + " " + outfile
48
49        process = subprocess.Popen(command, shell=True, stdout=subprocess.PIPE)
50        exit_status = process.wait()
51        if exit_status == 0:
52            pass
53        else:
54            stdout, stderr = process.communicate()
55            log.info("ERROR: executing R software did not succeed. Plots cannot be generated\n"\
56                      + str(stdout) + str(stderr) + str(exit_status))
57
58
59    def draw_dist_plot(self, merged_files):
60        '''
61        Draw distribution plots
62        '''
63        for (pheno_key, files) in merged_files.iteritems():
64            log.info("Saved distribution plot as " + self.inout.out + "distribution_sumlogs.P" + pheno_key + ".pdf")
65            r_variables = {"emp_file":files["empp"], "out_file":files["perm"], \
66                            "prefix" :self.inout.out, "pheno":pheno_key}
67            self.execute(_get_package_path() + "distribution_plot.R", r_variables)
68
69   
70    def draw_qq_plots(self, files, adjusted):
71        """Draw Quantile-Quantile plots from assoc files."""
72       
73        assoc_files_per_pheno = [values["assoc_files"] for values in files.itervalues()]
74        assoc_files = [ assoc_file for assoc_file in from_iterable(assoc_files_per_pheno)]
75         
76        pheno_nr =0
77        for assoc_file in assoc_files:
78            pheno_nr +=1
79            image_filename = self.inout.out + "assoc_qq_plot.P" +str(pheno_nr)+ ".pdf"
80            log.info("Saved QQ plot from association analysis as " + image_filename)
81            r_variables = {"filename":image_filename, "assocfile":assoc_file, \
82                           "header":os.path.basename(assoc_file), \
83                            "r_path":_get_package_path(), "adjusted":adjusted}
84            self.execute(_get_package_path() + "generate_single_qq.R", r_variables)
85
86
87    def draw_qq_plots_for_sets(self, files):
88        """Draw quantile-quantile plots for P values of gene sets."""
89        raw_p_values = dict([(pheno_key, values["raw_p_files"]) for pheno_key, values in files.iteritems()])
90        for (pheno_name, filename) in raw_p_values.iteritems():
91            image_filename = self.inout.out + "qq-plot_all_sets.P" + str(pheno_name) + ".pdf"
92            log.info("Saved QQ plot(s) from gene group(s) as " + image_filename)
93            r_variables = {"p_values_file" : filename , \
94                            "image_filename": image_filename, \
95                            "r_path":_get_package_path()}
96            self.execute(_get_package_path() + "multiple_qq_plot_for_sets.R", r_variables)
97
98
99def convert_list_to_arrayray(list_to_convert):
100    """
101    Convert a python list to a R array (this array is an string in python)
102    """
103    r_array_text = str(list_to_convert).replace("[", "c(").replace(']', ')')
104    return r_array_text
105
106
107def from_iterable(iterables):
108    """
109    this function is only provided in python 2.6 and higher
110    in the itertools.chain.from_iterable module.
111    This is a drop in replacement for this function
112    """
113    for iterabl in iterables:
114        for element in iterabl:
115            yield element
116           
117
118
119def _get_package_path():
120    """
121    get the absolute path of this file
122    """
123    return os.path.dirname(os.path.abspath(__file__)) + "/"   
Note: See TracBrowser for help on using the repository browser.