Mercurial > repos > earlhaminst > gstf_preparation
diff gstf_preparation.py @ 8:92f3966d5bc3 draft
planemo upload for repository https://github.com/TGAC/earlham-galaxytools/tree/master/tools/gstf_preparation commit 88ba62ae8c3d9587a0015c72209242ad0c1df0c2
author | earlhaminst |
---|---|
date | Wed, 16 May 2018 20:03:57 -0400 |
parents | 56bbdbfe3eaa |
children | f4acbfe8d6fe |
line wrap: on
line diff
--- a/gstf_preparation.py Wed Apr 25 11:06:03 2018 -0400 +++ b/gstf_preparation.py Wed May 16 20:03:57 2018 -0400 @@ -1,6 +1,5 @@ from __future__ import print_function -import collections import json import optparse import sqlite3 @@ -9,7 +8,23 @@ version = "0.4.0" gene_count = 0 -Sequence = collections.namedtuple('Sequence', ['header', 'sequence']) + +class Sequence(object): + def __init__(self, header, sequence_parts): + self.header = header + self.sequence_parts = sequence_parts + self._sequence = None + + @property + def sequence(self): + if self._sequence is None: + self._sequence = ''.join(self.sequence_parts) + return self._sequence + + def print(self, fh=sys.stdout): + print(self.header, file=fh) + for line in self.sequence_parts: + print(line, file=fh) def FASTAReader_gen(fasta_filename): @@ -25,8 +40,7 @@ while line and line[0] != '>': sequence_parts.append(line.rstrip()) line = fasta_file.readline() - sequence = "\n".join(sequence_parts) - yield Sequence(header, sequence) + yield Sequence(header, sequence_parts) def create_tables(conn): @@ -360,6 +374,10 @@ # Extract the transcript id by removing everything after the first space and then removing the version if it is an Ensembl id transcript_id = remove_id_version(entry.header[1:].lstrip().split(' ')[0]) + if len(entry.sequence) % 3 != 0: + print("Transcript '%s' in file '%s' has a coding sequence length which is not multiple of 3" % (transcript_id, fasta_arg), file=sys.stderr) + continue + gene_id = fetch_gene_id_for_transcript(conn, transcript_id) if not gene_id: print("Transcript '%s' in file '%s' not found in the gene feature information" % (transcript_id, fasta_arg), file=sys.stderr) @@ -383,6 +401,10 @@ if options.longestCDS and transcript_id not in selected_transcript_ids: continue + if len(entry.sequence) % 3 != 0: + print("Transcript '%s' in file '%s' has a coding sequence length which is not multiple of 3" % (transcript_id, fasta_arg), file=sys.stderr) + continue + species_for_transcript, seq_region_for_transcript = fetch_species_and_seq_region_for_transcript(conn, transcript_id) if not species_for_transcript: print("Transcript '%s' in file '%s' not found in the gene feature information" % (transcript_id, fasta_arg), file=sys.stderr) @@ -391,14 +413,12 @@ if options.headers: # Change the FASTA header to '>TranscriptId_species', as required by TreeBest # Remove any underscore in the species - header = ">%s_%s" % (transcript_id, species_for_transcript.replace('_', '')) - else: - header = entry.header + entry.header = ">%s_%s" % (transcript_id, species_for_transcript.replace('_', '')) if seq_region_for_transcript.lower() in regions: - filtered_fasta_file.write("%s\n%s\n" % (header, entry.sequence)) + entry.print(filtered_fasta_file) else: - output_fasta_file.write("%s\n%s\n" % (header, entry.sequence)) + entry.print(output_fasta_file) conn.close()