source: trunk/brs2010p35/src/jag/jag_main.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: 22.8 KB
Line 
1#!/usr/bin/env python
2
3'''
4Created on Dec 13, 2010
5
6@author: maartenk, estherlips
7
8main class where program is called with from command line
9'''
10
11import getopt
12import sys
13import os
14import logging as log
15import time
16
17sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)) + "/src/")
18#next is only for development reasons (when run as python jag_main.py)
19sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)) + "/../")
20
21import jag.common as common
22
23def _get_starttime():
24    """
25    get system start time
26    """
27    starttime = time.asctime(time.localtime(time.time()))
28    return("Analysis started: " + starttime)
29
30
31def _get_endtime():
32    """
33    get system end time
34    """
35    endtime = time.asctime(time.localtime(time.time()))
36    return("\nFinished analysis: " + endtime + "\n")
37   
38   
39def _get_terminated_time():
40    """
41    get system end time
42    """
43    termtime = time.asctime(time.localtime(time.time()))
44    return("\nAnalysis terminated: " + termtime + "\n")
45   
46   
47def _get_header_of_program():
48    """
49    print the header of the program
50    """
51   
52    return("""
53    +----------------------------------------------------------+
54    |        JAG      |version  SVN290 |  date: 2011-07-27     |
55    |----------------------------------------------------------|
56    |  For documentation, citation &   bug-report instructions:|
57    |        http://ctglab.nl/software/jag/                    |
58    +----------------------------------------------------------+
59    """)
60
61
62def _usage():
63    """
64    print help for using this program
65
66    """
67    return("""
68    === MAPPING SNPs TO GENES ===
69    --annotate         file with genes to map
70    --up               size upstream region in kb
71    --down             size downstream region in kb
72    --gene_boundaries  file with gene boundaries
73    --snp_file         file with snps
74     
75    === ASSOCIATION ANALYSIS ===
76    performs association analysis + permutations
77    --bfile    -b      path of fam, bim and bed file without file extension
78    --group    -g      file with gene groups & snps
79    --perm     -m      number of permutations
80    --no_emp           do not calculate empirical P-value
81   
82    == OPTION FOR ASSOCIATION TEST ==
83    --verbose  -v      print output of every plink file
84    --pheno    -p      add a file with phenotypes
85    --out      -o      set a prefix for output file
86    --logistic         add logistic switch to plink
87    --linear
88    --sex
89    --covar
90    --gene_based       groups are based on single genes
91    --adjust           adjust pvalues of real data with lambda
92    --no_graph         do not create graphs within R
93
94
95    === MERGE PERMUTATIONS ===
96    merge multiple output files when step 2 is split in multiple processes
97    --merge            merge multiple permutation files
98    --out              set a prefix to find the permuted results (optional)
99   
100    === GET RANDOM DRAWS ===
101    draw multiple random SNP or gene sets out of list with all genes and SNPs
102    obligatory
103    --ndraw            number of sets of random genes/SNP sets to be drawn (= integer )
104    --group            group file
105    --snp2gene         list with all genes and SNPs
106   
107    ==switch==
108    --ngenes           select functional set by set name
109    --neff_all         {set name}        select all SNPs default
110    --neff_intergenic  {set name}   only select SNPs in genes
111    --neff_genic       {set name}          only select SNPs not in genes
112
113    obligatory with --neff option:
114    --neff_calc
115    --bfile
116   
117    mandatory:
118    --exclude_group (default enable)
119    --include_group   
120
121    === CALCULATE EMPP OF EMPP ===
122    --orig_empp        file containing the empP on the real data
123    --control_empp     file containing the empP on the random draws
124    --gene_group       select functional set by set name
125
126   ===examples===:
127   association test:         ./jag -bfile hapmap/hapmap_CEU_r23_60 -group hapmap/hapmap_group.txt --perm=100
128                             ./jag -bfile hapmap/hapmap_CEU_r23_60 -group hapmap/hapmap_group.txt -pheno hapmap/hapmap_pheno.txt -m 100
129   self contained test       ./jag --merge
130   random draws              ./jag --ndraw 100 --group hapmap/hapmap_group.txt --neff_genic FG2 --snp2gene=all_coding_genes_mapped_up10.0_down10.0.txt --bfile=hapmap/hapmap_CEU_r23_60 --neff_calc empp.P1.out
131                             ./jag --ndraw 100 --group hapmap/hapmap_group.txt --ngenes FG2 --snp2gene=all_coding_genes_mapped_up10.0_down10.0.txt 
132   competitive test          ./jag --orig_empp test2.empp.merged.P1.out --control_empp draw_neff_genic.empp.P1.out --gene_group FG2
133    """)
134
135
136
137class Main:
138    """
139    set options and runs program
140    """
141   
142    def _step0_annotate_genes(self, out, annotate_file, up_size, down_size, gene_boundaries, snp_file):
143        """
144        --annotate initial groups
145        --up    size upstream region in kb
146        --down     size downstream region in kb
147        --gene-boundaries    file with gene boundaries
148        --snp_file    file with snps
149        """
150   
151        out_all = os.path.abspath(os.path.curdir) + "/all_coding_genes_mapped_up" + str(up_size) + "_down" + str(down_size) + ".txt"
152        import jag.snp2geneset as snp2geneset
153        all_gene2snp, group2snp = snp2geneset.snp2geneset(annotate_file, up_size, down_size, gene_boundaries, snp_file)
154        mapping_as_string = snp2geneset.format_mapping(group2snp)
155       
156        all_genes_as_string = snp2geneset.format_mapping(all_gene2snp)
157        file = out.save_text_to_filename("mapped_up" + str(up_size) + "_down" + str(down_size) + ".txt", mapping_as_string)
158   
159        common.save_text_file(out_all, all_genes_as_string)
160        log.info("Saved snps2all_genes file as " + out_all)
161        log.info("Saved snps2geneset file as " + file)
162
163
164    def _step1_run_association_analysis(self, group, plink, inoutput, gene_based, adjusted):
165        """
166        run association analysis
167        """         
168        log.info("Running association analysis using PLINK with command:")
169        plink.print_command = True
170        from jag.association_analysis import AssociationAnalysis
171        association_analysis = AssociationAnalysis(inoutput)
172        files = association_analysis.run_step1(plink, group, gene_based, adjusted)
173 
174     
175    def _step2_run_permutations(self, group, perm, seed, no_emp, plink, inoutput, gene_based): 
176        """
177        run permutations
178        """   
179       
180        log.info("\nRunning permutations...")
181        plink.print_command = False         
182        from jag.permutation import Permutation
183        permutation = Permutation(plink, group, inoutput)
184        permutation.no_emp = no_emp
185        permutation.localpermutation(perm, seed, gene_based) #run permutations on local machine
186       
187        if perm < 2:
188            log.info("\nNo distribution plot is generated since you should run at least 2 permutations to do so.")
189            inoutput.run_rproject = False
190           
191        if inoutput.run_rproject:
192            import jag.plot_with_r as plot_with_r
193            plotter = plot_with_r.call_r(inoutput)
194            plotter.draw_dist_plot(permutation.files)
195       
196       
197    def _step3_merge_results(self, inoutput):
198        """
199        merge multiple files of permutation results
200       
201        """
202        log.info ("Merging permutations results...")
203        import jag.mergepermutations as mergepermutations
204        merger = mergepermutations.MergePermutation(inoutput)
205        merger.mergeresults()
206       
207        if inoutput.run_rproject:
208            import jag.plot_with_r as plot_with_r
209            plotter = plot_with_r.call_r(inoutput)
210            plotter.draw_dist_plot(merger.files)
211
212        log.info ('Merging is done.')
213        log.info(_get_endtime())
214        sys.exit(0)
215
216
217    def _step4a_draw_genes(self, draw , ngenes, ndraws):
218        """
219        draw random sets of genes, based on number of genes in group of interest
220        """
221   
222        n = str(ndraws)
223        log.info("Drawing " + n + " random gene sets...")
224        import jag.drawrandom as drawrandom
225        genes = drawrandom.Genes(draw)
226        genes.genes(ngenes, ndraws)
227        log.info(_get_endtime())
228        sys.exit(0)
229
230
231    def _step4b_draw_snps(self, draw, neff, ndraws, draw_selection, plink, emp, out):
232        """
233        draw random sets of snps (based on calculated effective snps)
234        """
235        common.getfile_handle(emp, "Empirical Pvalue file not set \
236                                    correctly (this variable is set with --neff_calc )", False)
237
238        import jag.drawrandom as drawrandom
239        log.info ("Drawing random SNP sets (based on nEff)...")
240       
241        snps = drawrandom.Snp(draw)
242        snps.inregion = draw_selection
243
244        snps.snp(neff, ndraws, emp, plink, out)
245        log.info(_get_endtime())
246        sys.exit(0)
247
248   
249    def _calc_emp_p_of_emp_p(self, orig, sims, set_name):
250        """
251        calculation of empirical pvalue for competitive test
252        """
253        #check if all 3 variables are set
254        if sum([bool(item)for item in [orig, sims, set_name]]) == 3:
255            import jag.permutationresults as permutation
256            orig_header, orig_data = permutation.read_permutated_results(orig)
257            try:
258                row_setname = orig_data[0].index(set_name.strip())
259
260            except ValueError:
261                log.info("Could not find " + set_name + " in  the file " + orig)
262                log.info(_get_terminated_time())
263                sys.exit(1)
264
265            pvalue_orig = orig_data[orig_header.index("emp_p")][row_setname]
266
267            sims_header, sims_data = permutation.read_permutated_results(sims)
268            emp_p_sims = False
269            try:
270                emp_p_sims = sims_data[sims_header.index("emp_p")]
271            except ValueError:
272                log.info("Could not find emp_p in the file" + sims)
273                log.info(_get_terminated_time())
274                sys.exit(1)
275
276
277            log.info("Empirical P-value for competitive test = " + str(permutation.calc_emp_emp(pvalue_orig, emp_p_sims)))
278            log.info(_get_endtime())
279
280        else:
281            log.info("Original empp or set is not set")
282            log.info(_get_terminated_time())
283            sys.exit(1)
284           
285   
286    def extract_variables_from_command_line(self, args):
287        """
288        get variables from command line
289        """
290       
291        # if no variables > print usage
292        if len(args) == 0:
293            print _get_header_of_program()
294            print _get_starttime()
295            print "\nNo arguments are given to run JAG."
296            print "Type ./jag --help for usage."
297            print _get_terminated_time()
298            sys.exit(2)
299               
300        else:
301            try:
302                opts, args = getopt.getopt(args, "p:chb:m:vg:o:s:", ["ndraw=", "ngenes=", "neff_calc=",
303                    "out=", "snp2gene=", "neff_all=", "neff_intergenic=",
304                    "neff_genic=", "exclude_group",
305                    "annotate=", "snp_file=", "gene_boundaries=", "down=", "up=",
306                    "include_group", "no_emp", "gene_based", "no_graph", "help",
307                    "merge=",
308                    "verbose", "runplink=",
309                    "output=", "bfile=", "pheno=", "group=", "perm=", "cluster=", "seed=",
310                    "gene_group=",
311                    "orig_empp=", "control_empp=", "gene_group=",
312                    "covar=", "sex", "logistic", "linear", "adjust"])
313           
314            except getopt.GetoptError as err:
315            # print help information and exit:
316            # will print something like "option -assigned_value not recognized"
317            # here is print used and no log.info because this is not initialized
318                print _get_header_of_program()
319                print _get_starttime()
320                print ""
321                print str(err)
322                print "Type ./jag --help for usage."
323                print _get_terminated_time()
324                sys.exit(2)
325               
326        return opts
327
328
329    def enable_logging_with_prefix(self, inoutput):
330        """
331        #log with prefix
332        """
333        log.basicConfig(level=log.INFO, stream=sys.stdout, #format='%(levelname)s\t%(asctime)s   %(filename)s:%(lineno)d\t%(message)s',
334            format='%(message)s', datefmt='%H:%M:%S')
335        file_handler = log.FileHandler(inoutput.out + 'jag.log', mode='w')
336        log.root.addHandler(file_handler)
337       
338
339    def __init__(self, args):
340               
341        opts = self.extract_variables_from_command_line(args)
342       
343        from jag.plink import Plink
344        plink = Plink()
345                       
346        group = None
347        perm = None
348        seed = None
349
350        annotate_file = False
351        up = float(0)
352        down = float(0)
353        gene_boundaries_file = None
354        snp_file = None
355
356        no_emp = False
357        create_graphs = True
358        ngenes = False
359        empirical_p_filename = False
360
361        ndraws = 0
362        draw_selection = None
363        exclude_group = True
364        complete_gene_snp_mapping = None
365        verbose = False
366        gene_based = False
367       
368        merge = False
369        adjusted = False
370        prefix = None
371       
372        orig = False
373        sims = False
374        set_name = False
375       
376        #catch out prefix
377        out_prefix = None
378       
379        # open log file
380        for identifier, assigned_value in opts:
381            if identifier in ("-o", "--out"):
382                out_prefix = assigned_value
383         
384        #set input and output prefix first: this is also used for the log   
385        from jag.file_fetch_and_write import InAndOut
386        inoutput = InAndOut()
387        inoutput.set_outfile(out_prefix)
388       
389        #no logging before these statements
390        self.enable_logging_with_prefix(inoutput)
391       
392        #logging is now enabled
393        log.info(_get_header_of_program())
394        log.info("Writing this text to log file [" + inoutput.out + "jag.log]\n")
395        log.info(_get_starttime())
396       
397        # print command line options
398        log.info("\nJAG options in effect:")
399        for o, a in opts:
400            if o not in ("-h", "--help"):
401                log.info("\t" + o + " " + a)
402            else:
403                log.info("\t" + o + " " + a)
404                log.info("\nPrinting help documentation...")
405                log.info(_usage())
406                print(_get_terminated_time())
407                sys.exit(2)
408     
409        if len(opts) == 0:
410            _usage()
411           
412        else:
413           
414            for identifier, assigned_value in opts:
415
416                if identifier in ("-o", "--out"):
417                    #this stuff is printed to catch the out for the location of the log file
418                    out_prefix = assigned_value
419                                       
420                elif identifier in ("-g", "--group"):
421                    assert common.getfile_handle(assigned_value, \
422                    "Group file (set with -g or --group) is not set correct", verbose)
423                    group = assigned_value
424                   
425                elif identifier in ("-m", "--perm"):
426                    perm = int(assigned_value)
427                   
428                elif identifier in ("-v", "--verbose"):
429                    plink.verbose = True
430                    verbose = True
431                   
432                elif identifier in ("--snp_file"):
433                    snp_file = assigned_value
434                   
435                    assert common.getfile_handle(assigned_value, \
436                                                 "SNP file (--snp_file) is not set correctly.", verbose)
437
438                elif identifier in ("--control_empp"):
439                    assert common.getfile_handle(assigned_value, \
440                                                  "sims file not set correct(is set by --control_empp)", verbose)
441                    sims = assigned_value
442                   
443                elif identifier in ("--orig_empp"):
444                    assert common.getfile_handle(assigned_value, \
445                                                  "orig file not set correct(is set by --orig_empp)", verbose)
446                    orig = assigned_value
447                   
448                elif identifier in ("--gene_group"):
449                    set_name = assigned_value
450                   
451                elif identifier in ("--no_emp"):
452                    no_emp = True
453                    create_graphs = False
454                   
455                elif identifier in ("--no_graph"):
456                    create_graphs = False
457                   
458                elif identifier in ("--ndraw"):
459                    ndraws = int(assigned_value)
460                   
461                elif identifier in ("--ngenes"):
462                    ngenes = assigned_value
463                   
464                elif identifier in ("--annotate"):
465                    assert common.getfile_handle(assigned_value, \
466                                                 "Annotation file (--annotate) is not set correctly.", verbose)
467                    annotate_file = assigned_value
468                   
469                elif identifier in ("--up"):
470                    try:
471                        up = float(assigned_value)
472                       
473                    except ValueError:
474                        log.info("Value after --up parameter should be a number")
475                        sys.exit(1)
476                                       
477                elif identifier in ("--down"):
478                    try:
479                        down = float(assigned_value)
480                       
481                    except ValueError:
482                        log.info("Value after --down parameter should be a number")
483                        sys.exit(1)
484               
485                elif identifier in ("--gene_boundaries"):
486                    assert common.getfile_handle(assigned_value, "Gene boundary file is not set correct(is set by --gene_boundaries)", verbose)
487                    gene_boundaries_file = assigned_value
488                 
489                elif identifier in ("--snp2gene"):
490                    complete_gene_snp_mapping = assigned_value
491
492                elif identifier in ("--neff_genic"):
493                    draw_selection = "genic"
494                    set_name = assigned_value
495                   
496                elif identifier in ("--neff_intergenic"):
497                    draw_selection = "intergenic"
498                    set_name = assigned_value
499                   
500                elif identifier in ("--neff_all"):
501                    draw_selection = "all"
502                    set_name = assigned_value
503                   
504                elif identifier in ("--covar"):
505                    assert common.getfile_handle(assigned_value, \
506                         "covar file is not specified correctly (use --covar;", verbose)
507                    plink.covar_file = assigned_value
508                   
509                elif identifier in ("--sex"):
510                    plink.switches += "--sex "
511                   
512                elif identifier in ("--linear"):
513                    plink.switches += "--linear "
514                   
515                elif identifier in ("--logistic"):
516                    plink.switches += "--logistic "
517                   
518                elif identifier in ("--adjust"):
519                    plink.switches += "--adjust "
520                    adjusted = True
521
522                elif identifier in ("--exclude_group"):
523                    exclude_group = True
524                   
525                elif identifier in ("--include_group"):
526                    exclude_group = False
527
528                elif identifier in ("--neff_calc"):
529                    assert common.getfile_handle(assigned_value, "empirical_p_filename file is not \
530                    specified correctly (use --neff_calc;\
531                     default file iempirical_p_filenamemp.P1.out", verbose)
532                    empirical_p_filename = assigned_value
533                   
534                elif identifier in ("-b", "--bfile"):
535                    plink.bfile = common.check_bim_bed_fam(assigned_value)
536                   
537                elif identifier in ("-p", "--pheno"):
538                    assert common.getfile_handle(assigned_value, \
539                     "phenotype file not set correct(is set by -p or --pheno)", verbose)
540                    plink.pheno_file = assigned_value
541                   
542                elif identifier in ("-s", "--seed"):
543                    seed = assigned_value
544                   
545                elif identifier in ("--gene_based"):
546                    gene_based = True
547
548                elif identifier in ("--merge"):
549                    prefix = assigned_value
550                    inoutput = InAndOut()
551                    inoutput.set_outfile(prefix)
552                    merge = True
553                                   
554                else:
555                    assert False, "unhandled option"
556                   
557            log.info("")       
558                 
559            inoutput.run_rproject = create_graphs
560                       
561            if merge:
562                self._step3_merge_results(inoutput)
563
564            if orig or sims:
565                self._calc_emp_p_of_emp_p(orig , sims , set_name)
566                sys.exit(0)
567
568            if annotate_file and gene_boundaries_file and snp_file:
569                self._step0_annotate_genes(inoutput, annotate_file, up, down, \
570                                            gene_boundaries_file, snp_file)
571
572            elif ngenes or draw_selection:
573                common.getfile_handle(complete_gene_snp_mapping, "--snp2gene is not set correctly.", verbose)
574                common.getfile_handle(group, "--group is not set.", verbose)
575
576                import jag.drawrandom as drawrandom
577                draw = drawrandom.DrawRandom(complete_gene_snp_mapping, group, inoutput)
578
579                draw.exclude = exclude_group
580
581                if seed:
582                    draw.setseed(seed)
583
584                if ngenes:
585                    self._step4a_draw_genes(draw, ngenes, ndraws)
586               
587                if draw_selection:
588                    self._step4b_draw_snps(draw, set_name, ndraws, draw_selection, \
589                                           plink, empirical_p_filename, out_prefix)
590
591            elif(perm > 0 and plink.bfile and group):
592                #Run Step 1 + Step 2
593                self._step1_run_association_analysis(group, plink, inoutput, gene_based, adjusted)
594                self._step2_run_permutations(group, perm, seed, no_emp, plink, inoutput, gene_based)
595
596            elif(perm == 0):
597                # Run only association analysis
598                self._step1_run_association_analysis(group, plink, inoutput, gene_based, adjusted)
599                log.info("\nNo permutations will be proceeded since number of permutations is zero.")
600               
601            else:
602                log.info("You are using an invalid combination of parameters.\nCheck the help file for the correct combination.")
603
604        log.info(_get_endtime())
605           
606if __name__ == '__main__':
607    Main(sys.argv[1:])
Note: See TracBrowser for help on using the repository browser.