[2] | 1 | #!/usr/bin/env python |
---|
| 2 | """ |
---|
| 3 | usage: %prog $input $out_file1 |
---|
| 4 | -1, --cols=N,N,N,N: Columns for start, end, strand in input file |
---|
| 5 | -d, --dbkey=N: Genome build of input file |
---|
| 6 | -o, --output_format=N: the data type of the output file |
---|
| 7 | -g, --GALAXY_DATA_INDEX_DIR=N: the directory containing alignseq.loc |
---|
| 8 | -G, --gff: input and output file, when it is interval, coordinates are treated as GFF format (1-based, half-open) rather than 'traditional' 0-based, closed format. |
---|
| 9 | """ |
---|
| 10 | from galaxy import eggs |
---|
| 11 | import pkg_resources |
---|
| 12 | pkg_resources.require( "bx-python" ) |
---|
| 13 | import sys, string, os, re |
---|
| 14 | from bx.cookbook import doc_optparse |
---|
| 15 | import bx.seq.nib |
---|
| 16 | import bx.seq.twobit |
---|
| 17 | from galaxy.tools.util.galaxyops import * |
---|
| 18 | from galaxy.tools.util.gff_util import * |
---|
| 19 | |
---|
| 20 | assert sys.version_info[:2] >= ( 2, 4 ) |
---|
| 21 | |
---|
| 22 | def stop_err( msg ): |
---|
| 23 | sys.stderr.write( msg ) |
---|
| 24 | sys.exit() |
---|
| 25 | |
---|
| 26 | def reverse_complement( s ): |
---|
| 27 | complement_dna = {"A":"T", "T":"A", "C":"G", "G":"C", "a":"t", "t":"a", "c":"g", "g":"c", "N":"N", "n":"n" } |
---|
| 28 | reversed_s = [] |
---|
| 29 | for i in s: |
---|
| 30 | reversed_s.append( complement_dna[i] ) |
---|
| 31 | reversed_s.reverse() |
---|
| 32 | return "".join( reversed_s ) |
---|
| 33 | |
---|
| 34 | def check_seq_file( dbkey, GALAXY_DATA_INDEX_DIR ): |
---|
| 35 | seq_file = "%s/alignseq.loc" % GALAXY_DATA_INDEX_DIR |
---|
| 36 | seq_path = '' |
---|
| 37 | for line in open( seq_file ): |
---|
| 38 | line = line.rstrip( '\r\n' ) |
---|
| 39 | if line and not line.startswith( "#" ) and line.startswith( 'seq' ): |
---|
| 40 | fields = line.split( '\t' ) |
---|
| 41 | if len( fields ) < 3: |
---|
| 42 | continue |
---|
| 43 | if fields[1] == dbkey: |
---|
| 44 | seq_path = fields[2].strip() |
---|
| 45 | break |
---|
| 46 | return seq_path |
---|
| 47 | |
---|
| 48 | |
---|
| 49 | def __main__(): |
---|
| 50 | options, args = doc_optparse.parse( __doc__ ) |
---|
| 51 | try: |
---|
| 52 | chrom_col, start_col, end_col, strand_col = parse_cols_arg( options.cols ) |
---|
| 53 | dbkey = options.dbkey |
---|
| 54 | output_format = options.output_format |
---|
| 55 | gff_format = options.gff |
---|
| 56 | GALAXY_DATA_INDEX_DIR = options.GALAXY_DATA_INDEX_DIR |
---|
| 57 | input_filename, output_filename = args |
---|
| 58 | except: |
---|
| 59 | doc_optparse.exception() |
---|
| 60 | |
---|
| 61 | includes_strand_col = strand_col >= 0 |
---|
| 62 | strand = None |
---|
| 63 | nibs = {} |
---|
| 64 | twobits = {} |
---|
| 65 | seq_path = check_seq_file( dbkey, GALAXY_DATA_INDEX_DIR ) |
---|
| 66 | if not os.path.exists( seq_path ): |
---|
| 67 | # If this occurs, we need to fix the metadata validator. |
---|
| 68 | stop_err( "No sequences are available for '%s', request them by reporting this error." % dbkey ) |
---|
| 69 | |
---|
| 70 | skipped_lines = 0 |
---|
| 71 | first_invalid_line = 0 |
---|
| 72 | invalid_line = '' |
---|
| 73 | fout = open( output_filename, "w" ) |
---|
| 74 | warnings = [] |
---|
| 75 | warning = '' |
---|
| 76 | twobitfile = None |
---|
| 77 | |
---|
| 78 | for i, line in enumerate( open( input_filename ) ): |
---|
| 79 | line = line.rstrip( '\r\n' ) |
---|
| 80 | if line and not line.startswith( "#" ): |
---|
| 81 | fields = line.split( '\t' ) |
---|
| 82 | try: |
---|
| 83 | chrom = fields[chrom_col] |
---|
| 84 | start = int( fields[start_col] ) |
---|
| 85 | end = int( fields[end_col] ) |
---|
| 86 | if gff_format: |
---|
| 87 | start, end = convert_gff_coords_to_bed( [start, end] ) |
---|
| 88 | if includes_strand_col: |
---|
| 89 | strand = fields[strand_col] |
---|
| 90 | except: |
---|
| 91 | warning = "Invalid chrom, start or end column values. " |
---|
| 92 | warnings.append( warning ) |
---|
| 93 | skipped_lines += 1 |
---|
| 94 | if not invalid_line: |
---|
| 95 | first_invalid_line = i + 1 |
---|
| 96 | invalid_line = line |
---|
| 97 | continue |
---|
| 98 | if start > end: |
---|
| 99 | warning = "Invalid interval, start '%d' > end '%d'. " % ( start, end ) |
---|
| 100 | warnings.append( warning ) |
---|
| 101 | skipped_lines += 1 |
---|
| 102 | if not invalid_line: |
---|
| 103 | first_invalid_line = i + 1 |
---|
| 104 | invalid_line = line |
---|
| 105 | continue |
---|
| 106 | |
---|
| 107 | if strand not in ['+', '-']: |
---|
| 108 | strand = '+' |
---|
| 109 | sequence = '' |
---|
| 110 | |
---|
| 111 | if seq_path and os.path.exists( "%s/%s.nib" % ( seq_path, chrom ) ): |
---|
| 112 | if chrom in nibs: |
---|
| 113 | nib = nibs[chrom] |
---|
| 114 | else: |
---|
| 115 | nibs[chrom] = nib = bx.seq.nib.NibFile( file( "%s/%s.nib" % ( seq_path, chrom ) ) ) |
---|
| 116 | try: |
---|
| 117 | sequence = nib.get( start, end-start ) |
---|
| 118 | except: |
---|
| 119 | warning = "Unable to fetch the sequence from '%d' to '%d' for build '%s'. " %( start, end-start, dbkey ) |
---|
| 120 | warnings.append( warning ) |
---|
| 121 | skipped_lines += 1 |
---|
| 122 | if not invalid_line: |
---|
| 123 | first_invalid_line = i + 1 |
---|
| 124 | invalid_line = line |
---|
| 125 | continue |
---|
| 126 | elif seq_path and os.path.isfile( seq_path ): |
---|
| 127 | if not(twobitfile): |
---|
| 128 | twobitfile = bx.seq.twobit.TwoBitFile( file( seq_path ) ) |
---|
| 129 | try: |
---|
| 130 | sequence = twobitfile[chrom][start:end] |
---|
| 131 | except: |
---|
| 132 | warning = "Unable to fetch the sequence from '%d' to '%d' for build '%s'. " %( start, end-start, dbkey ) |
---|
| 133 | warnings.append( warning ) |
---|
| 134 | skipped_lines += 1 |
---|
| 135 | if not invalid_line: |
---|
| 136 | first_invalid_line = i + 1 |
---|
| 137 | invalid_line = line |
---|
| 138 | continue |
---|
| 139 | else: |
---|
| 140 | warning = "Chromosome by name '%s' was not found for build '%s'. " % ( chrom, dbkey ) |
---|
| 141 | warnings.append( warning ) |
---|
| 142 | skipped_lines += 1 |
---|
| 143 | if not invalid_line: |
---|
| 144 | first_invalid_line = i + 1 |
---|
| 145 | invalid_line = line |
---|
| 146 | continue |
---|
| 147 | if sequence == '': |
---|
| 148 | warning = "Chrom: '%s', start: '%s', end: '%s' is either invalid or not present in build '%s'. " %( chrom, start, end, dbkey ) |
---|
| 149 | warnings.append( warning ) |
---|
| 150 | skipped_lines += 1 |
---|
| 151 | if not invalid_line: |
---|
| 152 | first_invalid_line = i + 1 |
---|
| 153 | invalid_line = line |
---|
| 154 | continue |
---|
| 155 | if includes_strand_col and strand == "-": |
---|
| 156 | sequence = reverse_complement( sequence ) |
---|
| 157 | |
---|
| 158 | if output_format == "fasta" : |
---|
| 159 | l = len( sequence ) |
---|
| 160 | c = 0 |
---|
| 161 | if gff_format: |
---|
| 162 | start, end = convert_bed_coords_to_gff( [ start, end ] ) |
---|
| 163 | fields = [dbkey, str( chrom ), str( start ), str( end ), strand] |
---|
| 164 | meta_data = "_".join( fields ) |
---|
| 165 | fout.write( ">%s\n" % meta_data ) |
---|
| 166 | while c < l: |
---|
| 167 | b = min( c + 50, l ) |
---|
| 168 | fout.write( "%s\n" % str( sequence[c:b] ) ) |
---|
| 169 | c = b |
---|
| 170 | else: # output_format == "interval" |
---|
| 171 | meta_data = "\t".join( fields ) |
---|
| 172 | if gff_format: |
---|
| 173 | format_str = "%s seq \"%s\";\n" |
---|
| 174 | else: |
---|
| 175 | format_str = "%s\t%s\n" |
---|
| 176 | fout.write( format_str % ( meta_data, str( sequence ) ) ) |
---|
| 177 | |
---|
| 178 | fout.close() |
---|
| 179 | |
---|
| 180 | if warnings: |
---|
| 181 | warn_msg = "%d warnings, 1st is: " % len( warnings ) |
---|
| 182 | warn_msg += warnings[0] |
---|
| 183 | print warn_msg |
---|
| 184 | if skipped_lines: |
---|
| 185 | print 'Skipped %d invalid lines, 1st is #%d, "%s"' % ( skipped_lines, first_invalid_line, invalid_line ) |
---|
| 186 | |
---|
| 187 | if __name__ == "__main__": __main__() |
---|