Mercurial > repos > vimalkumarvelayudhan > riboseqr_wrapper
comparison riboseqr/ribosome_profile.py @ 5:423ad61697c4
Bugfix 1: [triplet] Lengths (frameCounting) if given should be a range (not zero).
readingFrame function fails with subscript out of bounds.
Bugfix 2: [triplet] Check if transcript name in SAM matches the name in FASTA. Produce
an error message if it's not. This fixes the problem where an empty plot is
produced (no bars).
[ribosome_profile] - A proper error message is now produced if an invalid transcript name is
provided.
Updated test data
| author | Vimalkumar Velayudhan <vimalkumarvelayudhan@gmail.com> |
|---|---|
| date | Tue, 27 Oct 2015 12:21:50 +0000 |
| parents | c34c364ce75d |
| children |
comparison
equal
deleted
inserted
replaced
| 2:b2eb07000039 | 5:423ad61697c4 |
|---|---|
| 7 import rpy2.robjects as robjects | 7 import rpy2.robjects as robjects |
| 8 import utils | 8 import utils |
| 9 | 9 |
| 10 rscript = '' | 10 rscript = '' |
| 11 R = robjects.r | 11 R = robjects.r |
| 12 | |
| 13 | |
| 14 class RibosomeProfileError(Exception): | |
| 15 pass | |
| 12 | 16 |
| 13 | 17 |
| 14 def run_rscript(command=None): | 18 def run_rscript(command=None): |
| 15 """Run R command, log it, append to rscript""" | 19 """Run R command, log it, append to rscript""" |
| 16 global rscript | 20 global rscript |
| 43 </head> | 47 </head> |
| 44 <body> | 48 <body> |
| 45 """ | 49 """ |
| 46 html += '<h2>Plot ribosome profile - results</h2>\n<hr>\n' | 50 html += '<h2>Plot ribosome profile - results</h2>\n<hr>\n' |
| 47 if len(transcript_name): | 51 if len(transcript_name): |
| 52 # get a list of all transcript names as in SAM | |
| 53 refseqs = R('levels(riboDat@riboGR$`RiboSeq file 1`@seqnames)') | |
| 54 refseqs = tuple(refseqs) | |
| 55 if transcript_name not in refseqs: | |
| 56 raise RibosomeProfileError('Transcript "{}" does not exist in SAM alignment. ' | |
| 57 'Please provide the transcript name exactly as in the SAM ' | |
| 58 'alignment.'.format(transcript_name)) | |
| 48 cmd_args = ( | 59 cmd_args = ( |
| 49 '"{transcript_name}", main="{transcript_name}",' | 60 '"{transcript_name}", main="{transcript_name}",' |
| 50 'coordinates=ffCs@CDS, riboData=riboDat,' | 61 'coordinates=ffCs@CDS, riboData=riboDat,' |
| 51 'length={transcript_length}'.format(**options)) | 62 'length={transcript_length}'.format(**options)) |
| 52 if transcript_cap: | 63 if transcript_cap: |