/tools/maf/maf_stats.py

  1#!/usr/bin/env python
  2#Dan Blankenberg
  3"""
  4Reads a list of intervals and a maf. Outputs a new set of intervals with statistics appended.
  5"""
  6
  7import sys
  8from galaxy import eggs
  9import pkg_resources; pkg_resources.require( "bx-python" )
 10import bx.intervals.io
 11from bx.bitset import BitSet
 12from galaxy.tools.util import maf_utilities
 13
 14assert sys.version_info[:2] >= ( 2, 4 )
 15
 16def __main__():
 17    maf_source_type = sys.argv.pop( 1 )
 18    input_maf_filename = sys.argv[1].strip()
 19    input_interval_filename = sys.argv[2].strip()
 20    output_filename = sys.argv[3].strip()
 21    dbkey = sys.argv[4].strip()
 22    try:
 23        chr_col  = int( sys.argv[5].strip() ) - 1
 24        start_col = int( sys.argv[6].strip() ) - 1
 25        end_col = int( sys.argv[7].strip() ) - 1
 26    except:
 27        print >>sys.stderr, "You appear to be missing metadata. You can specify your metadata by clicking on the pencil icon associated with your interval file."
 28        sys.exit()
 29    summary = sys.argv[8].strip()
 30    if summary.lower() == "true": summary = True
 31    else: summary = False
 32
 33    mafIndexFile = "%s/maf_index.loc" % sys.argv[9]
 34    try:
 35        maf_index_filename = sys.argv[10].strip()
 36    except:
 37        maf_index_filename = None
 38    index = index_filename = None
 39    if maf_source_type == "user":
 40        #index maf for use here
 41        index, index_filename = maf_utilities.open_or_build_maf_index( input_maf_filename, maf_index_filename, species = [dbkey] )
 42        if index is None:
 43            print >>sys.stderr, "Your MAF file appears to be malformed."
 44            sys.exit()
 45    elif maf_source_type == "cached":
 46        #access existing indexes
 47        index = maf_utilities.maf_index_by_uid( input_maf_filename, mafIndexFile )
 48        if index is None:
 49            print >> sys.stderr, "The MAF source specified (%s) appears to be invalid." % ( input_maf_filename )
 50            sys.exit()
 51    else:
 52        print >>sys.stdout, 'Invalid source type specified: %s' % maf_source_type 
 53        sys.exit()
 54        
 55    out = open(output_filename, 'w')
 56    
 57    num_region = None
 58    num_bad_region = 0
 59    species_summary = {}
 60    total_length = 0
 61    #loop through interval file
 62    for num_region, region in enumerate( bx.intervals.io.NiceReaderWrapper( open( input_interval_filename, 'r' ), chrom_col = chr_col, start_col = start_col, end_col = end_col, fix_strand = True, return_header = False, return_comments = False ) ):
 63        src = "%s.%s" % ( dbkey, region.chrom )
 64        region_length = region.end - region.start
 65        if region_length < 1:
 66            num_bad_region += 1
 67            continue
 68        total_length += region_length
 69        coverage = { dbkey: BitSet( region_length ) }
 70        
 71        
 72        for block in index.get_as_iterator( src, region.start, region.end ):
 73            for spec in maf_utilities.get_species_in_block( block ):
 74                if spec not in coverage: coverage[spec] = BitSet( region_length )
 75            for block in maf_utilities.iter_blocks_split_by_species( block ):
 76                if maf_utilities.component_overlaps_region( block.get_component_by_src( src ), region ):
 77                    #need to chop and orient the block
 78                    block = maf_utilities.orient_block_by_region( maf_utilities.chop_block_by_region( block, src, region ), src, region, force_strand = '+' )
 79                    start_offset, alignment = maf_utilities.reduce_block_by_primary_genome( block, dbkey, region.chrom, region.start )
 80                    for i in range( len( alignment[dbkey] ) ):
 81                        for spec, text in alignment.items():
 82                            if text[i] != '-':
 83                                coverage[spec].set( start_offset + i )
 84        if summary:
 85            #record summary
 86            for key in coverage.keys():
 87                if key not in species_summary: species_summary[key] = 0
 88                species_summary[key] = species_summary[key] + coverage[key].count_range()
 89        else:
 90            #print coverage for interval
 91            coverage_sum = coverage[dbkey].count_range()
 92            out.write( "%s\t%s\t%s\t%s\n" % ( "\t".join( region.fields ), dbkey, coverage_sum, region_length - coverage_sum ) )
 93            keys = coverage.keys()
 94            keys.remove( dbkey )
 95            keys.sort()
 96            for key in keys:
 97                coverage_sum = coverage[key].count_range()
 98                out.write( "%s\t%s\t%s\t%s\n" % ( "\t".join( region.fields ), key, coverage_sum, region_length - coverage_sum ) )
 99    if summary:
100        out.write( "#species\tnucleotides\tcoverage\n" )
101        for spec in species_summary:
102            out.write( "%s\t%s\t%.4f\n" % ( spec, species_summary[spec], float( species_summary[spec] ) / total_length ) )
103    out.close()
104    if num_region is not None:
105        print "%i regions were processed with a total length of %i." % ( num_region + 1, total_length )
106    if num_bad_region:
107        print "%i regions were invalid." % ( num_bad_region )
108    maf_utilities.remove_temp_index_file( index_filename )
109
110if __name__ == "__main__": __main__()