view riboseqr/ribosome_profile.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
# TODO: check first if argument parameters are of the right type if given.
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 plot_transcript(rdata_load='Metagene.rda', transcript_name='',
                    transcript_length='27', transcript_cap='',
                    html_file='Plot-ribosome-profile.html',
                    output_path=os.getcwd()):
    """Plot ribosome profile for a given transcript. """
    options = {}
    for key, value, rtype, rmode in (
            ('transcript_name', transcript_name, 'str', None),
            ('transcript_length', transcript_length, 'int', 'charvector'),
            ('transcript_cap', transcript_cap, 'int', None)):
        options[key] = utils.process_args(value, ret_type=rtype, ret_mode=rmode)

    run_rscript('library(riboSeqR)')
    run_rscript('load("{}")'.format(rdata_load))

    html = '<h2>Plot ribosome profile - results</h2><hr>'
    if len(transcript_name):
        cmd_args = (
            '"{transcript_name}", main="{transcript_name}",'
            'coordinates=ffCs@CDS, riboData=riboDat,'
            'length={transcript_length}'.format(**options))
        if transcript_cap:
            cmd_args += ', cap={transcript_cap}'.format(**options)
        cmd = 'plotTranscript({})'.format(cmd_args)

        for fmat in ('pdf', 'png'):
            run_rscript('{0}(file="{1}")'.format(
                fmat, os.path.join(
                    output_path, 'Ribosome-profile-plot.{}'.format(fmat))))
            run_rscript(cmd)
            run_rscript('dev.off()')

        html += ('<p>Selected ribosome footprint length: <strong>{0}</strong>'
                 '</p><p><img border="1" src="Ribosome-profile-plot.png" '
                 'alt="Ribosome profile plot" /><br>'
                 '<a href="Ribosome-profile-plot.pdf">PDF version</a>'
                 '</p>'.format(transcript_length))
    else:
        msg = 'No transcript name was provided. Did not generate plot.'
        html += '<p>{}</p>'.format(msg)
        logging.debug(msg)

    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, 'ribosome-profile.R'), 'w') as r:
        r.write(rscript)

    html += ('<h4>R script for this session</h4>'
             '<p><a href="ribosome-profile.R">ribosome-profile.R</a></p>')

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


if __name__ == '__main__':
    parser = argparse.ArgumentParser(description='Plot Ribosome profile')

    # required arguments
    flags = parser.add_argument_group('required arguments')
    flags.add_argument('--rdata_load', required=True,
                       help='Saved riboSeqR data from Step 2')
    flags.add_argument('--transcript_name', required=True,
                       help='Name of the transcript to be plotted')
    flags.add_argument(
        '--transcript_length', required=True,
        help='Size class of ribosome footprint data to be plotted',
        default='27')
    flags.add_argument(
        '--transcript_cap', required=True,
        help=('Cap on the largest value that will be plotted as an abundance '
              'of the ribosome footprint data'))
    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='%(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)

    plot_transcript(rdata_load=args.rdata_load,
                    transcript_name=args.transcript_name,
                    transcript_length=args.transcript_length,
                    transcript_cap=args.transcript_cap,
                    html_file=args.html_file, output_path=args.output_path)
    logging.debug('Done!')