From 8082c6f36fd50e51ce642041b2733ea86c592b3e Mon Sep 17 00:00:00 2001 From: Daniel Blankenberg Date: Tue, 23 Feb 2010 16:48:07 -0500 Subject: [PATCH] Add a new FASTQ tool suite. Four FASTQ variants are supported: sanger, illumina, solexa and solid. Tools include: FASTQ Groomer convert between various FASTQ quality formats Combine FASTA and QUAL into FASTQ FASTQ joiner on paired end reads FASTQ splitter on joined paired end reads FASTQ to FASTA converter FASTQ Summary Statistics by column Filter FASTQ reads by quality score and length FASTQ Trimmer by column Manipulate FASTQ reads on various attributes Boxplot of quality statistics (Generic, with outliers) --- datatypes_conf.xml.sample | 7 +- .../interval_to_bedstrict_converter.py | 2 +- lib/galaxy/datatypes/qualityscore.py | 19 +- lib/galaxy/datatypes/sequence.py | 12 + lib/galaxy_utils/__init__.py | 0 lib/galaxy_utils/sequence/__init__.py | 0 lib/galaxy_utils/sequence/fastq.py | 702 ++++++++++++++++++ lib/galaxy_utils/sequence/sequence.py | 61 ++ lib/galaxy_utils/sequence/transform.py | 74 ++ test-data/empty_file.dat | 0 tool_conf.xml.sample | 13 +- tools/fastq/fastq_combiner.py | 45 ++ tools/fastq/fastq_combiner.xml | 55 ++ tools/fastq/fastq_filter.py | 34 + tools/fastq/fastq_filter.xml | 311 ++++++++ tools/fastq/fastq_groomer.py | 37 + tools/fastq/fastq_groomer.xml | 349 +++++++++ tools/fastq/fastq_manipulation.py | 37 + tools/fastq/fastq_manipulation.xml | 386 ++++++++++ tools/fastq/fastq_paired_end_joiner.py | 38 + tools/fastq/fastq_paired_end_joiner.xml | 55 ++ tools/fastq/fastq_paired_end_splitter.py | 33 + tools/fastq/fastq_paired_end_splitter.xml | 56 ++ tools/fastq/fastq_stats.py | 48 ++ tools/fastq/fastq_stats.xml | 72 ++ tools/fastq/fastq_to_fasta.py | 21 + tools/fastq/fastq_to_fasta.xml | 33 + tools/fastq/fastq_trimmer.py | 41 + tools/fastq/fastq_trimmer.xml | 113 +++ tools/metag_tools/split_paired_reads.xml | 4 +- tools/plotting/boxplot.xml | 72 ++ 31 files changed, 2721 insertions(+), 9 deletions(-) create mode 100644 lib/galaxy_utils/__init__.py create mode 100644 lib/galaxy_utils/sequence/__init__.py create mode 100644 lib/galaxy_utils/sequence/fastq.py create mode 100644 lib/galaxy_utils/sequence/sequence.py create mode 100644 lib/galaxy_utils/sequence/transform.py create mode 100644 test-data/empty_file.dat create mode 100644 tools/fastq/fastq_combiner.py create mode 100644 tools/fastq/fastq_combiner.xml create mode 100644 tools/fastq/fastq_filter.py create mode 100644 tools/fastq/fastq_filter.xml create mode 100644 tools/fastq/fastq_groomer.py create mode 100644 tools/fastq/fastq_groomer.xml create mode 100644 tools/fastq/fastq_manipulation.py create mode 100644 tools/fastq/fastq_manipulation.xml create mode 100644 tools/fastq/fastq_paired_end_joiner.py create mode 100644 tools/fastq/fastq_paired_end_joiner.xml create mode 100644 tools/fastq/fastq_paired_end_splitter.py create mode 100644 tools/fastq/fastq_paired_end_splitter.xml create mode 100644 tools/fastq/fastq_stats.py create mode 100644 tools/fastq/fastq_stats.xml create mode 100644 tools/fastq/fastq_to_fasta.py create mode 100644 tools/fastq/fastq_to_fasta.xml create mode 100644 tools/fastq/fastq_trimmer.py create mode 100644 tools/fastq/fastq_trimmer.xml create mode 100644 tools/plotting/boxplot.xml diff --git a/datatypes_conf.xml.sample b/datatypes_conf.xml.sample index 6de6a2cac94..ef93ffe358c 100644 --- a/datatypes_conf.xml.sample +++ b/datatypes_conf.xml.sample @@ -1,6 +1,6 @@ - + @@ -34,6 +34,9 @@ + + + @@ -60,7 +63,9 @@ + + diff --git a/lib/galaxy/datatypes/converters/interval_to_bedstrict_converter.py b/lib/galaxy/datatypes/converters/interval_to_bedstrict_converter.py index de5c7bfda5d..4445a28b103 100644 --- a/lib/galaxy/datatypes/converters/interval_to_bedstrict_converter.py +++ b/lib/galaxy/datatypes/converters/interval_to_bedstrict_converter.py @@ -47,7 +47,7 @@ def __main__(): #does file already conform to bed strict? #if so, we want to keep extended columns, otherwise we'll create a generic 6 column bed file strict_bed = True - if extension == 'bed' and ( chromCol, startCol, endCol, nameCol, strandCol ) == ( 0, 1, 2, 3, 5 ): + if extension in [ 'bed', 'bedstrict' ] and ( chromCol, startCol, endCol, nameCol, strandCol ) == ( 0, 1, 2, 3, 5 ): for count, line in enumerate( open( input_name ) ): line = line.strip() if line == "" or line.startswith("#"): diff --git a/lib/galaxy/datatypes/qualityscore.py b/lib/galaxy/datatypes/qualityscore.py index 138b863bc4d..e51be3267b8 100644 --- a/lib/galaxy/datatypes/qualityscore.py +++ b/lib/galaxy/datatypes/qualityscore.py @@ -9,7 +9,13 @@ from galaxy import util log = logging.getLogger(__name__) -class QualityScoreSOLiD ( data.Text ): +class QualityScore ( data.Text ): + """ + until we know more about quality score formats + """ + file_ext = "qual" + +class QualityScoreSOLiD ( QualityScore ): """ until we know more about quality score formats """ @@ -58,7 +64,7 @@ class QualityScoreSOLiD ( data.Text ): pass return False -class QualityScore454 ( data.Text ): +class QualityScore454 ( QualityScore ): """ until we know more about quality score formats """ @@ -97,9 +103,14 @@ class QualityScore454 ( data.Text ): pass return False -class QualityScoreSolexa ( data.Text ): +class QualityScoreSolexa ( QualityScore ): """ until we know more about quality score formats """ file_ext = "qualsolexa" - \ No newline at end of file + +class QualityScoreIllumina ( QualityScore ): + """ + until we know more about quality score formats + """ + file_ext = "qualillumina" diff --git a/lib/galaxy/datatypes/sequence.py b/lib/galaxy/datatypes/sequence.py index 747c7734aec..8dc665ebb8d 100644 --- a/lib/galaxy/datatypes/sequence.py +++ b/lib/galaxy/datatypes/sequence.py @@ -216,6 +216,18 @@ class FastqSanger( Fastq ): """Class representing a FASTQ sequence ( the Sanger variant )""" file_ext = "fastqsanger" +class FastqSolexa( Fastq ): + """Class representing a FASTQ sequence ( the Solexa variant )""" + file_ext = "fastqsolexa" + +class FastqIllumina( Fastq ): + """Class representing a FASTQ sequence ( the Illumina 1.3+ variant )""" + file_ext = "fastqillumina" + +class FastqSolid( Fastq ): + """Class representing a FASTQ sequence ( the SOLiD (color space) variant )""" + file_ext = "fastqsolid" + try: from galaxy import eggs import pkg_resources; pkg_resources.require( "bx-python" ) diff --git a/lib/galaxy_utils/__init__.py b/lib/galaxy_utils/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/lib/galaxy_utils/sequence/__init__.py b/lib/galaxy_utils/sequence/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/lib/galaxy_utils/sequence/fastq.py b/lib/galaxy_utils/sequence/fastq.py new file mode 100644 index 00000000000..250d85422ce --- /dev/null +++ b/lib/galaxy_utils/sequence/fastq.py @@ -0,0 +1,702 @@ +#Dan Blankenberg +import math +import string +import transform +from sequence import SequencingRead + +class fastqSequencingRead( SequencingRead ): + format = 'sanger' #sanger is default + ascii_min = 33 + ascii_max = 126 + quality_min = 0 + quality_max = 93 + score_system = 'phred' #phred or solexa + sequence_space = 'base' #base or color + @classmethod + def get_class_by_format( cls, format ): + assert format in FASTQ_FORMATS, 'Unknown format type specified: %s' % format + return FASTQ_FORMATS[ format ] + @classmethod + def convert_score_phred_to_solexa( cls, decimal_score_list ): + def phred_to_solexa( score ): + if score <= 0: #can't take log10( 1 - 1 ); make <= 0 into -5 + return -5 + return int( round( 10.0 * math.log10( math.pow( 10.0, ( float( score ) / 10.0 ) ) - 1.0 ) ) ) + return map( phred_to_solexa, decimal_score_list ) + @classmethod + def convert_score_solexa_to_phred( cls, decimal_score_list ): + def solexa_to_phred( score ): + return int( round( 10.0 * math.log10( math.pow( 10.0, ( float( score ) / 10.0 ) ) + 1.0 ) ) ) + return map( solexa_to_phred, decimal_score_list ) + @classmethod + def restrict_scores_to_valid_range( cls, decimal_score_list ): + def restrict_score( score ): + return max( min( score, cls.quality_max ), cls.quality_min ) + return map( restrict_score, decimal_score_list ) + @classmethod + def convert_base_to_color_space( cls, sequence ): + return cls.color_space_converter.to_color_space( sequence ) + @classmethod + def convert_color_to_base_space( cls, sequence ): + return cls.color_space_converter.to_base_space( sequence ) + def is_ascii_encoded( self ): + return ' ' not in self.quality #as per fastq definition only decimal quality strings can have spaces in them (and must have a trailing space) + def get_ascii_quality_scores( self ): + if self.is_ascii_encoded(): + return list( self.quality ) + else: + quality = self.quality.rstrip() #decimal scores should have a trailing space + if quality: + try: + return [ chr( int( val ) + self.ascii_min - self.quality_min ) for val in quality.split( ' ' ) ] + except ValueError, e: + raise ValueError( 'Error Parsing quality String. ASCII quality strings cannot contain spaces (%s): %s' % ( self.quality, e ) ) + else: + return [] + def get_decimal_quality_scores( self ): + if self.is_ascii_encoded(): + return [ ord( val ) - self.ascii_min + self.quality_min for val in self.quality ] + else: + quality = self.quality.rstrip() #decimal scores should have a trailing space + if quality: + return map( int, quality.split( ' ' ) ) + else: + return [] + def convert_read_to_format( self, format, force_quality_encoding = None ): + assert format in FASTQ_FORMATS, 'Unknown format type specified: %s' % format + assert force_quality_encoding in [ None, 'ascii', 'decimal' ], 'Invalid force_quality_encoding: %s' % force_quality_encoding + new_class = FASTQ_FORMATS[ format ] + new_read = new_class() + new_read.identifier = self.identifier + if self.sequence_space == new_class.sequence_space: + new_read.sequence = self.sequence + else: + if self.sequence_space == 'base': + new_read.sequence = self.convert_base_to_color_space( self.sequence ) + else: + new_read.sequence = self.convert_color_to_base_space( self.sequence ) + new_read.description = self.description + if self.score_system != new_read.score_system: + if self.score_system == 'phred': + score_list = self.convert_score_phred_to_solexa( self.get_decimal_quality_scores() ) + else: + score_list = self.convert_score_solexa_to_phred( self.get_decimal_quality_scores() ) + else: + score_list = self.get_decimal_quality_scores() + new_read.quality = "%s " % " ".join( map( str, new_class.restrict_scores_to_valid_range( score_list ) ) ) #need trailing space to be valid decimal fastq + if force_quality_encoding is None: + if self.is_ascii_encoded(): + new_encoding = 'ascii' + else: + new_encoding = 'decimal' + else: + new_encoding = force_quality_encoding + if new_encoding == 'ascii': + new_read.quality = "".join( new_read.get_ascii_quality_scores() ) + return new_read + def get_sequence( self ): + return self.sequence + def slice( self, left_column_offset, right_column_offset ): + new_read = fastqSequencingRead.get_class_by_format( self.format )() + new_read.identifier = self.identifier + new_read.sequence = self.get_sequence()[left_column_offset:right_column_offset] + new_read.description = self.description + if self.is_ascii_encoded(): + new_read.quality = self.quality[left_column_offset:right_column_offset] + else: + quality = map( str, self.get_decimal_quality_scores()[left_column_offset:right_column_offset] ) + if quality: + new_read.quality = "%s " % " ".join( quality ) + else: + new_read.quality = '' + return new_read + def is_valid_format( self ): + if self.is_ascii_encoded(): + for val in self.get_ascii_quality_scores(): + val = ord( val ) + if val < self.ascii_min or val > self.ascii_max: + return False + else: + for val in self.get_decimal_quality_scores(): + if val < self.quality_min or val > self.quality_max: + return False + if not self.is_valid_sequence(): + return False + return True + def is_valid_sequence( self ): + for base in self.get_sequence(): + if base not in self.valid_sequence_list: + return False + return True + def insufficient_quality_length( self ): + return len( self.get_ascii_quality_scores() ) < len( self.sequence ) + def assert_sequence_quality_lengths( self ): + qual_len = len( self.get_ascii_quality_scores() ) + seq_len = len( self.sequence ) + assert qual_len == seq_len, "Invalid FASTQ file: quality score length (%i) does not match sequence length (%i)" % ( qual_len, seq_len ) + def reverse( self, clone = True ): + #need to override how decimal quality scores are reversed + if clone: + rval = self.clone() + else: + rval = self + rval.sequence = transform.reverse( self.sequence ) + if rval.is_ascii_encoded(): + rval.quality = rval.quality[::-1] + else: + rval.quality = reversed( rval.get_decimal_quality_scores() ) + rval.quality = "%s " % " ".join( map( str, rval.quality ) ) + return rval + +class fastqSangerRead( fastqSequencingRead ): + format = 'sanger' + ascii_min = 33 + ascii_max = 126 + quality_min = 0 + quality_max = 93 + score_system = 'phred' + sequence_space = 'base' + +class fastqIlluminaRead( fastqSequencingRead ): + format = 'illumina' + ascii_min = 64 + ascii_max = 126 + quality_min = 0 + quality_max = 62 + score_system = 'phred' + sequence_space = 'base' + +class fastqSolexaRead( fastqSequencingRead ): + format = 'solexa' + ascii_min = 59 + ascii_max = 126 + quality_min = -5 + quality_max = 62 + score_system = 'solexa' + sequence_space = 'base' + +class fastqSolidRead( fastqSequencingRead ): + format = 'solid' #color space + ascii_min = 33 + ascii_max = 126 + quality_min = 0 + quality_max = 93 + score_system = 'phred' + sequence_space = 'color' + valid_sequence_list = map( str, range( 7 ) ) + [ '.' ] + def __len__( self ): + if self.has_adapter_base(): #Adapter base is not counted in length of read + return len( self.sequence ) - 1 + return fastqSequencingRead.__len__( self ) + def has_adapter_base( self ): + if self.sequence and self.sequence[0] in string.letters: #adapter base must be a letter + return True + return False + def insufficient_quality_length( self ): + if self.has_adapter_base(): + return len( self.get_ascii_quality_scores() ) + 1 < len( self.sequence ) + return fastqSequencingRead.insufficient_quality_length( self ) + def assert_sequence_quality_lengths( self ): + if self.has_adapter_base(): + qual_len = len( self.get_ascii_quality_scores() ) + seq_len = len( self.sequence ) + assert qual_len + 1 == seq_len, "Invalid FASTQ file: quality score length (%i) does not match sequence length (%i with adapter base)" % ( qual_len, seq_len ) + else: + return fastqSequencingRead.assert_sequence_quality_lengths( self ) + def get_sequence( self ): + if self.has_adapter_base(): + return self.sequence[1:] + return self.sequence + def reverse( self, clone = True ): + #need to override how color space is reversed + if clone: + rval = self.clone() + else: + rval = self + if rval.has_adapter_base(): + adapter = rval.sequence[0] + #sequence = rval.sequence[1:] + rval.sequence = self.color_space_converter.to_color_space( transform.reverse( self.color_space_converter.to_base_space( rval.sequence ) ), adapter_base = adapter ) + else: + rval.sequence = transform.reverse( rval.sequence ) + + if rval.is_ascii_encoded(): + rval.quality = rval.quality[::-1] + else: + rval.quality = reversed( rval.get_decimal_quality_scores() ) + rval.quality = "%s " % " ".join( map( str, rval.quality ) ) + return rval + def complement( self, clone = True ): + #need to override how color space is complemented + if clone: + rval = self.clone() + else: + rval = self + if rval.has_adapter_base(): #No adapter, color space stays the same + adapter = rval.sequence[0] + sequence = rval.sequence[1:] + if adapter.lower() != 'u': + adapter = transform.DNA_complement( adapter ) + else: + adapter = transform.RNA_complement( adapter ) + rval.sequence = "%s%s" % ( adapter, sequence ) + return rval + def change_adapter( self, new_adapter, clone = True ): + #if new_adapter is empty, remove adapter, otherwise replace with new_adapter + if clone: + rval = self.clone() + else: + rval = self + if rval.has_adapter_base(): + if new_adapter: + if new_adapter != rval.sequence[0]: + rval.sequence = rval.color_space_converter.to_color_space( rval.color_space_converter.to_base_space( rval.sequence ), adapter_base = new_adapter ) + else: + rval.sequence = rval.sequence[1:] + elif new_adapter: + rval.sequence = "%s%s" % ( new_adapter, rval.sequence ) + return rval + + +FASTQ_FORMATS = {} +for format in [ fastqIlluminaRead, fastqSolexaRead, fastqSangerRead, fastqSolidRead ]: + FASTQ_FORMATS[ format.format ] = format + + +class fastqAggregator(): + VALID_FORMATS = FASTQ_FORMATS.keys() + def __init__( self, ): + self.ascii_values_used = [] #quick lookup of all ascii chars used + self.seq_lens = {} #counts of seqs by read len + self.nuc_index_quality = [] #counts of scores by read column + self.nuc_index_base = [] #counts of bases by read column + def consume_read( self, fastq_read ): + #ascii values used + for val in fastq_read.get_ascii_quality_scores(): + if val not in self.ascii_values_used: + self.ascii_values_used.append( val ) + #lengths + seq_len = len( fastq_read ) + self.seq_lens[ seq_len ] = self.seq_lens.get( seq_len, 0 ) + 1 + #decimal qualities by column + for i, val in enumerate( fastq_read.get_decimal_quality_scores() ): + if i == len( self.nuc_index_quality ): + self.nuc_index_quality.append( {} ) + self.nuc_index_quality[ i ][ val ] = self.nuc_index_quality[ i ].get( val, 0 ) + 1 + #bases by column + for i, nuc in enumerate( fastq_read.get_sequence() ): + if i == len( self.nuc_index_base ): + self.nuc_index_base.append( {} ) + nuc = nuc.upper() + self.nuc_index_base[ i ][ nuc ] = self.nuc_index_base[ i ].get( nuc, 0 ) + 1 + def get_valid_formats( self, check_list = None ): + if not check_list: + check_list = self.VALID_FORMATS + rval = [] + sequence = [] + for nuc_dict in self.nuc_index_base: + for nuc in nuc_dict.keys(): + if nuc not in sequence: + sequence.append( nuc ) + sequence = "".join( sequence ) + quality = "".join( self.ascii_values_used ) + for fastq_format in check_list: + fastq_read = fastqSequencingRead.get_class_by_format( fastq_format )() + fastq_read.quality = quality + fastq_read.sequence = sequence + if fastq_read.is_valid_format(): + rval.append( fastq_format ) + return rval + def get_ascii_range( self ): + return ( min( self.ascii_values_used ), max( self.ascii_values_used ) ) + def get_decimal_range( self ): + decimal_values_used = [] + for scores in self.nuc_index_quality: + decimal_values_used.extend( scores.keys() ) + return ( min( decimal_values_used ), max( decimal_values_used ) ) + def get_length_counts( self ): + return self.seq_lens + def get_max_read_length( self ): + return len( self.nuc_index_quality ) + def get_read_count_for_column( self, column ): + if column >= len( self.nuc_index_quality ): + return 0 + return sum( self.nuc_index_quality[ column ].values() ) + def get_read_count( self ): + return self.get_read_count_for_column( 0 ) + def get_base_counts_for_column( self, column ): + return self.nuc_index_base[ column ] + def get_score_list_for_column( self, column ): + return self.nuc_index_quality[ column ].keys() + def get_score_min_for_column( self, column ): + return min( self.nuc_index_quality[ column ].keys() ) + def get_score_max_for_column( self, column ): + return max( self.nuc_index_quality[ column ].keys() ) + def get_score_sum_for_column( self, column ): + return sum( score * count for score, count in self.nuc_index_quality[ column ].iteritems() ) + def get_score_at_position_for_column( self, column, position ): + score_value_dict = self.nuc_index_quality[ column ] + scores = sorted( score_value_dict.keys() ) + for score in scores: + if score_value_dict[ score ] <= position: + position -= score_value_dict[ score ] + else: + return score + def get_summary_statistics_for_column( self, i ): + def _get_med_pos( size ): + halfed = int( size / 2 ) + if size % 2 == 1: + return [ halfed ] + return[ halfed - 1, halfed ] + read_count = self.get_read_count_for_column( i ) + + min_score = self.get_score_min_for_column( i ) + max_score = self.get_score_max_for_column( i ) + sum_score = self.get_score_sum_for_column( i ) + mean_score = float( sum_score ) / float( read_count ) + #get positions + med_pos = _get_med_pos( read_count ) + if 0 in med_pos: + q1_pos = [ 0 ] + q3_pos = [ read_count - 1 ] + else: + q1_pos = _get_med_pos( min( med_pos ) ) + q3_pos = [] + for pos in q1_pos: + q3_pos.append( max( med_pos ) + 1 + pos ) + #get scores at position + med_score = float( sum( [ self.get_score_at_position_for_column( i, pos ) for pos in med_pos ] ) ) / float( len( med_pos ) ) + q1 = float( sum( [ self.get_score_at_position_for_column( i, pos ) for pos in q1_pos ] ) ) / float( len( q1_pos ) ) + q3 = float( sum( [ self.get_score_at_position_for_column( i, pos ) for pos in q3_pos ] ) ) / float( len( q3_pos ) ) + #determine iqr and step + iqr = q3 - q1 + step = 1.5 * iqr + + #Determine whiskers and outliers + outliers = [] + score_list = sorted( self.get_score_list_for_column( i ) ) + left_whisker = q1 - step + for score in score_list: + if left_whisker <= score: + left_whisker = score + break + else: + outliers.append( score ) + + right_whisker = q3 + step + score_list.reverse() + for score in score_list: + if right_whisker >= score: + right_whisker = score + break + else: + outliers.append( score ) + + column_stats = { 'read_count': read_count, + 'min_score': min_score, + 'max_score': max_score, + 'sum_score': sum_score, + 'mean_score': mean_score, + 'q1': q1, + 'med_score': med_score, + 'q3': q3, + 'iqr': iqr, + 'left_whisker': left_whisker, + 'right_whisker': right_whisker, + 'outliers': outliers } + return column_stats + +class fastqReader( object ): + def __init__( self, fh, format = 'sanger' ): + self.file = fh + self.format = format + def close( self ): + return self.file.close() + def next(self): + while True: + fastq_header = self.file.readline() + if not fastq_header: + raise StopIteration + fastq_header = fastq_header.rstrip( '\n\r' ) + #remove empty lines, apparently extra new lines at end of file is common? + if fastq_header: + break + + assert fastq_header.startswith( '@' ), 'Invalid fastq header: %s' % fastq_header + rval = fastqSequencingRead.get_class_by_format( self.format )() + rval.identifier = fastq_header + while True: + line = self.file.readline() + if not line: + raise Exception( 'Invalid FASTQ file: could not parse second instance of sequence identifier.' ) + line = line.rstrip( '\n\r' ) + if line.startswith( '+' ) and ( len( line ) == 1 or line[1:].startswith( fastq_header[1:] ) ): + rval.description = line + break + rval.append_sequence( line ) + while rval.insufficient_quality_length(): + line = self.file.readline() + if not line: + break + rval.append_quality( line ) + rval.assert_sequence_quality_lengths() + return rval + def __iter__( self ): + while True: + yield self.next() + +class fastqNamedReader( object ): + def __init__( self, fh, format = 'sanger' ): + self.file = fh + self.format = format + self.reader = fastqReader( self.file, self.format ) + #self.last_offset = self.file.tell() + self.offset_dict = {} + self.eof = False + def close( self ): + return self.file.close() + def get( self, sequence_id ): + rval = None + if sequence_id in self.offset_dict: + initial_offset = self.file.tell() + seq_offset = self.offset_dict[ sequence_id ].pop( 0 ) + if not self.offset_dict[ sequence_id ]: + del self.offset_dict[ sequence_id ] + self.file.seek( seq_offset ) + rval = self.reader.next() + #assert rval.identifier == sequence_id, 'seq id mismatch' #should be able to remove this + self.file.seek( initial_offset ) + else: + while True: + offset = self.file.tell() + try: + fastq_read = self.reader.next() + except StopIteration: + self.eof = True + break #eof, id not found, will return None + if fastq_read.identifier == sequence_id: + rval = fastq_read + break + else: + if fastq_read.identifier not in self.offset_dict: + self.offset_dict[ fastq_read.identifier ] = [] + self.offset_dict[ fastq_read.identifier ].append( offset ) + return rval + def has_data( self ): + #returns a string representation of remaining data, or empty string (False) if no data remaining + eof = self.eof + count = 0 + rval = '' + if self.offset_dict: + count = sum( map( len, self.offset_dict.values() ) ) + if not eof: + offset = self.file.tell() + try: + fastq_read = self.reader.next() + except StopIteration: + eof = True + self.file.seek( offset ) + if count: + rval = "There were %i known sequence reads not utilized. " + if not eof: + rval = "%s%s" % ( rval, "An additional unknown number of reads exist in the input that were not utilized." ) + return rval + +class fastqWriter( object ): + def __init__( self, fh, format = None, force_quality_encoding = None ): + self.file = fh + self.format = format + self.force_quality_encoding = force_quality_encoding + def write( self, fastq_read ): + if self.format: + fastq_read = fastq_read.convert_read_to_format( self.format, force_quality_encoding = self.force_quality_encoding ) + self.file.write( str( fastq_read ) ) + def close( self ): + return self.file.close() + +class fastaWriter( object ): + def __init__( self, fh ): + self.file = fh + def write( self, fastq_read ): + #this will include SOLiD adapter base if applicable + self.file.write( ">%s\n%s\n" % ( fastq_read.identifier[1:], fastq_read.sequence ) ) + def close( self ): + return self.file.close() + +class fastqJoiner( object ): + def __init__( self, format, force_quality_encoding = None ): + self.format = format + self.force_quality_encoding = force_quality_encoding + def join( self, read1, read2 ): + if read1.identifier.endswith( '/2' ) and read2.identifier.endswith( '/1' ): + #swap 1 and 2 + tmp = read1 + read1 = read2 + read2 = tmp + del tmp + if read1.identifier.endswith( '/1' ) and read2.identifier.endswith( '/2' ): + identifier = read1.identifier[:-2] + else: + identifier = read1.identifier + + #use force quality encoding, if not present force to encoding of first read + force_quality_encoding = self.force_quality_encoding + if not force_quality_encoding: + if read1.is_ascii_encoded(): + force_quality_encoding = 'ascii' + else: + force_quality_encoding = 'decimal' + + new_read1 = read1.convert_read_to_format( self.format, force_quality_encoding = force_quality_encoding ) + new_read2 = read2.convert_read_to_format( self.format, force_quality_encoding = force_quality_encoding ) + rval = FASTQ_FORMATS[ self.format ]() + rval.identifier = identifier + if len( read1.description ) > 1: + rval.description = "+%s" % ( identifier[1:] ) + else: + rval.description = '+' + if rval.sequence_space == 'color': + #need to handle color space joining differently + #convert to nuc space, join, then convert back + rval.sequence = rval.convert_base_to_color_space( new_read1.convert_color_to_base_space( new_read1.sequence ) + new_read2.convert_color_to_base_space( new_read2.sequence ) ) + else: + rval.sequence = new_read1.sequence + new_read2.sequence + if force_quality_encoding == 'ascii': + rval.quality = new_read1.quality + new_read2.quality + else: + rval.quality = "%s %s" % ( new_read1.quality.strip(), new_read2.quality.strip() ) + return rval + def get_paired_identifier( self, fastq_read ): + identifier = fastq_read.identifier + if identifier[-2] == '/': + if identifier[-1] == "1": + identifier = "%s2" % identifier[:-1] + elif identifier[-1] == "2": + identifier = "%s1" % identifier[:-1] + return identifier + +class fastqSplitter( object ): + def split( self, fastq_read ): + length = len( fastq_read ) + #Only reads of even lengths can be split + if length % 2 != 0: + return None, None + half = int( length / 2 ) + read1 = fastq_read.slice( 0, half ) + read1.identifier += "/1" + if len( read1.description ) > 1: + read1.description += "/1" + read2 = fastq_read.slice( half, None ) + read2.identifier += "/2" + if len( read2.description ) > 1: + read2.description += "/2" + return read1, read2 + +class fastaSequence( ): + def __init__( self ): + self.identifier = None + self.sequence = '' #holds raw sequence string: no whitespace + def __len__( self ): + return len( self.sequence ) + def __str__( self ): + return "%s\n%s\n" % ( self.identifier, self.sequence ) + +class fastaReader( object ): + def __init__( self, fh ): + self.file = fh + def close( self ): + return self.file.close() + def next( self ): + line = self.file.readline() + #remove header comment lines + while line and line.startswith( '#' ): + line = self.file.readline() + if not line: + raise StopIteration + assert line.startswith( '>' ), "FASTA headers must start with >" + rval = fastaSequence() + rval.identifier = line.strip() + offset = self.file.tell() + while True: + line = self.file.readline() + if not line or line.startswith( '>' ): + if line: + self.file.seek( offset ) + return rval + #454 qual test data that was used has decimal scores that don't have trailing spaces + #so we'll need to parse and build these sequences not based upon de facto standards + #i.e. in a less than ideal fashion + line = line.rstrip() + if ' ' in rval.sequence or ' ' in line: + rval.sequence = "%s%s " % ( rval.sequence, line ) + else: + rval.sequence += line + offset = self.file.tell() + def __iter__( self ): + while True: + yield self.next() + +class fastaNamedReader( object ): + def __init__( self, fh ): + self.file = fh + self.reader = fastaReader( self.file ) + self.offset_dict = {} + self.eof = False + def close( self ): + return self.file.close() + def get( self, sequence_id ): + rval = None + if sequence_id in self.offset_dict: + initial_offset = self.file.tell() + seq_offset = self.offset_dict[ sequence_id ].pop( 0 ) + if not self.offset_dict[ sequence_id ]: + del self.offset_dict[ sequence_id ] + self.file.seek( seq_offset ) + rval = self.reader.next() + self.file.seek( initial_offset ) + else: + while True: + offset = self.file.tell() + try: + fasta_seq = self.reader.next() + except StopIteration: + self.eof = True + break #eof, id not found, will return None + if fasta_seq.identifier == sequence_id: + rval = fasta_seq + break + else: + if fasta_seq.identifier not in self.offset_dict: + self.offset_dict[ fasta_seq.identifier ] = [] + self.offset_dict[ fasta_seq.identifier ].append( offset ) + return rval + def has_data( self ): + #returns a string representation of remaining data, or empty string (False) if no data remaining + eof = self.eof + count = 0 + rval = '' + if self.offset_dict: + count = sum( map( len, self.offset_dict.values() ) ) + if not eof: + offset = self.file.tell() + try: + fasta_seq = self.reader.next() + except StopIteration: + eof = True + self.file.seek( offset ) + if count: + rval = "There were %i known sequences not utilized. " % count + if not eof: + rval = "%s%s" % ( rval, "An additional unknown number of sequences exist in the input that were not utilized." ) + return rval + +class fastqCombiner( object ): + def __init__( self, format ): + self.format = format + def combine(self, fasta_seq, quality_seq ): + fastq_read = fastqSequencingRead.get_class_by_format( self.format )() + fastq_read.identifier = "@%s" % fasta_seq.identifier[1:] + fastq_read.description = '+' + fastq_read.sequence = fasta_seq.sequence + fastq_read.quality = quality_seq.sequence + return fastq_read diff --git a/lib/galaxy_utils/sequence/sequence.py b/lib/galaxy_utils/sequence/sequence.py new file mode 100644 index 00000000000..6e775fcf1db --- /dev/null +++ b/lib/galaxy_utils/sequence/sequence.py @@ -0,0 +1,61 @@ +#Dan Blankenberg +import transform +import string +from copy import deepcopy + +class SequencingRead( object ): + color_space_converter = transform.ColorSpaceConverter() + valid_sequence_list = string.letters + def __init__( self ): + self.identifier = None + self.sequence = '' #holds raw sequence string: no whitespace + self.description = None + self.quality = '' #holds raw quality string: no whitespace, unless this contains decimal scores + def __len__( self ): + return len( self.sequence ) + def __str__( self ): + return "%s\n%s\n%s\n%s\n" % ( self.identifier, self.sequence, self.description, self.quality ) + def append_sequence( self, sequence ): + self.sequence += sequence.rstrip( '\n\r' ) + def append_quality( self, quality ): + self.quality += quality.rstrip( '\n\r' ) + def is_DNA( self ): + return 'u' not in self.sequence.lower() + def clone( self ): + return deepcopy( self ) + def reverse( self, clone = True ): + if clone: + rval = self.clone() + else: + rval = self + rval.sequence = transform.reverse( self.sequence ) + rval.quality = rval.quality[::-1] + return rval + def complement( self, clone = True ): + if clone: + rval = self.clone() + else: + rval = self + if rval.is_DNA(): + rval.sequence = transform.DNA_complement( rval.sequence ) + else: + rval.sequence = transform.RNA_complement( rval.sequence ) + return rval + def reverse_complement( self, clone = True ): + #need to reverse first, then complement + rval = self.reverse( clone = clone ) + return rval.complement( clone = False ) #already working with a clone if requested + def sequence_as_DNA( self, clone = True ): + if clone: + rval = self.clone() + else: + rval = self + rval.sequence = transform.to_DNA( rval.sequence ) + return rval + def sequence_as_RNA( self, clone = True ): + if clone: + rval = self.clone() + else: + rval = self + rval.sequence = transform.to_RNA( rval.sequence ) + return rval diff --git a/lib/galaxy_utils/sequence/transform.py b/lib/galaxy_utils/sequence/transform.py new file mode 100644 index 00000000000..f2871b8c431 --- /dev/null +++ b/lib/galaxy_utils/sequence/transform.py @@ -0,0 +1,74 @@ +#Dan Blankenberg +#Contains methods to tranform sequence strings +import string + +#FIXME: This should Handle ambiguity codes... +#Translation table for reverse Complement +DNA_COMPLEMENT = string.maketrans( "ACGTacgt", "TGCAtgca" ) +RNA_COMPLEMENT = string.maketrans( "ACGUacgu", "UGCAugca" ) +#Translation table for DNA <--> RNA +DNA_TO_RNA = string.maketrans( "Tt", "Uu" ) +RNA_TO_DNA = string.maketrans( "Uu", "Tt" ) + +#reverse sequence string +def reverse( sequence ): + return sequence[::-1] +#complement DNA sequence string +def DNA_complement( sequence ): + return sequence.translate( DNA_COMPLEMENT ) +#complement RNA sequence string +def RNA_complement( sequence ): + return sequence.translate( RNA_COMPLEMENT ) +#returns the reverse complement of the sequence +def DNA_reverse_complement( self, sequence ): + sequence = reverse( sequence ) + return DNA_complement( sequence ) +def RNA_reverse_complement( self, sequence ): + sequence = reverse( sequence ) + return RNA_complement( sequence ) +def to_DNA( sequence ): + return sequence.translate( DNA_TO_RNA ) +def to_RNA( sequence ): + return sequence.translate( RNA_TO_DNA ) + +class ColorSpaceConverter( object ): + unknown_base = 'N' + unknown_color = '.' + color_to_base_dict = {} + color_to_base_dict[ 'A' ] = { '0':'A', '1':'C', '2':'G', '3':'T', '4':'N', '5':'N', '6':'N', '.':'N' } + color_to_base_dict[ 'C' ] = { '0':'C', '1':'A', '2':'T', '3':'G', '4':'N', '5':'N', '6':'N', '.':'N' } + color_to_base_dict[ 'G' ] = { '0':'G', '1':'T', '2':'A', '3':'C', '4':'N', '5':'N', '6':'N', '.':'N' } + color_to_base_dict[ 'T' ] = { '0':'T', '1':'G', '2':'C', '3':'A', '4':'N', '5':'N', '6':'N', '.':'N' } + color_to_base_dict[ 'N' ] = { '0':'N', '1':'N', '2':'N', '3':'N', '4':'N', '5':'N', '6':'N', '.':'N' } + base_to_color_dict = {} + for base, color_dict in color_to_base_dict.iteritems(): + base_to_color_dict[ base ] = {} + for key, value in color_dict.iteritems(): + base_to_color_dict[ base ][ value ] = key + base_to_color_dict[ base ][ 'N' ] = '4' #force ACGT followed by N to be '4', because this is now 'processed' data; we could force to '.' (non-processed data) also + base_to_color_dict[ 'N' ].update( { 'A':'5', 'C':'5', 'G':'5', 'T':'5', 'N':'6' } ) + def __init__( self, fake_adapter_base = 'G' ): + assert fake_adapter_base in self.base_to_color_dict, 'A bad fake adapter base was provided: %s.' % fake_adapter_base + self.fake_adapter_base = fake_adapter_base + def to_color_space( self, sequence, adapter_base = None ): + if adapter_base is None: + adapter_base = self.fake_adapter_base + last_base = adapter_base #we add a fake adapter base so that the sequence can be decoded properly again + rval = last_base + for base in sequence: + rval += self.base_to_color_dict.get( last_base, self.base_to_color_dict[ self.unknown_base ] ).get( base, self.unknown_color ) + last_base = base + return rval + def to_base_space( self, sequence ): + if not isinstance( sequence, list ): + sequence = list( sequence ) + if sequence: + last_base = sequence.pop( 0 ) + else: + last_base = None + assert last_base in self.color_to_base_dict, 'A valid adapter base must be included when converting to base space from color space. Found: %s' % last_base + rval = '' + for color_val in sequence: + last_base = self.color_to_base_dict[ last_base ].get( color_val, self.unknown_base ) + rval += last_base + return rval diff --git a/test-data/empty_file.dat b/test-data/empty_file.dat new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tool_conf.xml.sample b/tool_conf.xml.sample index 7176cd4d528..d1ba80744ff 100644 --- a/tool_conf.xml.sample +++ b/tool_conf.xml.sample @@ -127,6 +127,7 @@ + @@ -174,7 +175,17 @@
-