view riboseqr/metagene.py @ 0:e01de823e919 draft default tip

Uploaded
author biotechcoder
date Fri, 01 May 2015 05:41:51 -0400
parents
children
line wrap: on
line source

#!/usr/bin/env python
import os
import sys
import argparse
import logging
import rpy2.robjects as robjects

import utils

rscript = ''
R = robjects.r


def run_rscript(command=None):
    """Run R command, log it, append to rscript"""
    global rscript
    if not command:
        return
    logging.debug(command)
    rscript += '{}\n'.format(command)
    R(command)


def do_analysis(
        rdata_load='Periodicity.rda', selected_lengths='27',
        selected_frames='', hit_mean='10', unique_hit_mean='1',
        ratio_check='TRUE', min5p='-20', max5p='200', min3p='-200', max3p='20',
        cap='', plot_title='', plot_lengths='27', rdata_save='Metagene.rda',
        html_file='Metagene-report.html', output_path=os.getcwd()):
    """Metagene analysis from saved periodicity R data file. """
    run_rscript('library(riboSeqR)')
    run_rscript('load("{}")'.format(rdata_load))

    logging.debug('fS\n{}\nfCs\n{}\n'.format(R['fS'], R['fCs']))
    options = {}
    for key, value, rtype, rmode in (
            ('lengths', selected_lengths, 'int', 'charvector'),
            ('frames', selected_frames, 'int', 'listvector'),
            ('hit_mean', hit_mean, 'int', None),
            ('unique_hit_mean', unique_hit_mean, 'int', None),
            ('ratio_check', ratio_check, 'bool', None),
            ('min5p', min5p, 'int', None), ('max5p', max5p, 'int', None),
            ('min3p', min3p, 'int', None), ('max3p', max3p, 'int', None),
            ('cap', cap, 'int', None),
            ('plot_title', plot_title, 'str', 'charvector'),
            ('plot_lengths', plot_lengths, 'int', 'list')):
            options[key] = utils.process_args(
                value, ret_type=rtype, ret_mode=rmode)

    cmd_args = """fCs, lengths={lengths},
    frames={frames}, hitMean={hit_mean},
    unqhitMean={unique_hit_mean}, fS=fS""".format(**options)

    if ratio_check == 'TRUE':
        cmd_args += ', ratioCheck = TRUE'

    run_rscript('ffCs <- filterHits({})'.format(cmd_args))
    logging.debug("ffCs\n{}\n".format(R['ffCs']))

    cds_args = ('coordinates=ffCs@CDS, riboDat=riboDat, min5p={min5p}, '
                'max5p={max5p}, min3p={min3p}, max3p={max3p}'.format(**options))

    if options['cap']:
        cds_args += ', cap={cap}'.format(**options)

    if options['plot_title']:
        cds_args += ', main={plot_title}'.format(**options)

    html = '<h2>Metagene analysis - results</h2><hr>'
    html += ('<p>Lengths of footprints used in analysis - <strong>'
             '<code>{0}</code></strong><br>Lengths of footprints selected for '
             'the plot - <strong><code>{1}</code></strong>'
             '</p>'.format(selected_lengths, plot_lengths))
    print 'plot lengths', options['plot_lengths']
    for count, length in enumerate(options['plot_lengths']):
        count += 1
        plot_file = 'Metagene-analysis-plot{0}'.format(count)
        for fmat in ('pdf', 'png'):
            run_rscript('{0}(file="{1}")'.format(
                fmat, os.path.join(
                    output_path, '{0}.{1}'.format(plot_file, fmat))))
            run_rscript('plotCDS({0},{1})'.format(
                cds_args, 'lengths={}'.format(length)))
            run_rscript('dev.off()')

        html += ('<h4>Length: {0}</h4><p><img border="1" src="{1}.png" '
                 'alt="Metagene analysis plot" /><br><a href="{1}.pdf">'
                 'PDF version</a></p>'.format(length, plot_file))
    run_rscript(
        'save("ffCs", "riboDat", "fastaCDS", file="{}", compress=FALSE)'.format(rdata_save))

    logging.info('\n{:#^80}\n{}\n{:#^80}\n'.format(
        ' R script for this session ', rscript, ' End R script '))

    with open(os.path.join(output_path, 'metagene.R'), 'w') as r:
        r.write(rscript)

    html += ('<h4>R script for this session</h4>'
             '<p><a href="metagene.R">metagene.R</a></p>'
             '<p>Next step: <em>Plot Ribosome profile</em></p>')

    with open(html_file, 'w') as f:
        f.write(html)

if __name__ == '__main__':

    parser = argparse.ArgumentParser(description='Metagene analysis')

    # required arguments
    flags = parser.add_argument_group('required arguments')
    flags.add_argument('--rdata_load', required=True,
                       help='Saved riboSeqR data from Periodicity step.')
    flags.add_argument('--selected_lengths', required=True,
                       help='Select frame lengths to filter. Comma-separated',
                       default='27')
    flags.add_argument(
        '--selected_frames', required=True,
        help='Select frames corresponding to frame lengths. Comma-separated')

    flags.add_argument(
        '--hit_mean', required=True,
        help='Mean number of hits within the replicate group for filtering',
        default='10')

    flags.add_argument(
        '--unique_hit_mean', required=True,
        help='Mean number of unique sequences within the replicate group '
             'for filtering', default='1')

    parser.add_argument(
        '--rdata_save', help='File to write R data to (default: %(default)s)',
        default='Metagene.rda')

    parser.add_argument(
        '--ratio_check',
        help='Check the ratios of the expected phase to maximal phase '
             'within the putative coding sequence (default: %(default)s)',
        choices=['TRUE', 'FALSE'], default='TRUE')

    parser.add_argument(
        '--plot_lengths',
        help='Length of footprints to be plotted. Multiple values should be '
             'comma-separated. In that case, multiple plots will be produced'
             '(default: %(default)s)', default='27')

    parser.add_argument(
        '--min5p',
        help='The distance upstream of the translation start to be plotted '
             '(default: %(default)s)', default='-20')

    parser.add_argument(
        '--max5p',
        help='The distance downstream of the translation start to be plotted '
             '(default: %(default)s)', default='200')

    parser.add_argument(
        '--min3p',
        help='The distance upstream of the translation end to be plotted '
             '(default: %(default)s)', default='-200')

    parser.add_argument(
        '--max3p',
        help='The distance downtream of the translation end to be plotted '
             '(default: %(default)s)', default='20')

    parser.add_argument(
        '--cap', help='If given, caps the height of plotted values '
                      '(default: %(default)s)')

    parser.add_argument('--plot_title', help='Title of the plot', default='')
    parser.add_argument('--html_file', help='HTML file with reports')
    parser.add_argument('--output_path', help='Directory to save output files')
    parser.add_argument(
        '--debug', help='Produce debug output', action='store_true')

    args = parser.parse_args()
    if args.debug:
        level = logging.DEBUG
    else:
        level = logging.INFO

    logging.basicConfig(format='%(module)s: %(levelname)s - %(message)s',
                        level=level, stream=sys.stdout)
    logging.debug('Supplied Arguments\n{}\n'.format(vars(args)))

    if not os.path.exists(args.output_path):
        os.mkdir(args.output_path)

    do_analysis(
        rdata_load=args.rdata_load, selected_lengths=args.selected_lengths,
        selected_frames=args.selected_frames, hit_mean=args.hit_mean,
        unique_hit_mean=args.unique_hit_mean, ratio_check=args.ratio_check,
        min5p=args.min5p, max5p=args.max5p, min3p=args.min3p, max3p=args.max3p,
        cap=args.cap, plot_title=args.plot_title,
        plot_lengths=args.plot_lengths, rdata_save=args.rdata_save,
        html_file=args.html_file, output_path=args.output_path)

    logging.info('Done!')