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 @@
-