Mercurial > repos > devteam > substitution_rates
annotate substitution_rates.py @ 1:412a03dc24a2 draft default tip
planemo upload commit fa9a2569c9ece4b03f622d0b85b0e5376a25084c-dirty
| author | devteam |
|---|---|
| date | Tue, 29 Sep 2015 11:08:38 -0400 |
| parents | d51dd9e4b517 |
| children |
| rev | line source |
|---|---|
| 0 | 1 #!/usr/bin/env python |
| 2 #guruprasad Ananda | |
| 3 """ | |
| 4 Estimates substitution rates from pairwise alignments using JC69 model. | |
| 5 """ | |
| 6 | |
| 7 from galaxy import eggs | |
|
1
412a03dc24a2
planemo upload commit fa9a2569c9ece4b03f622d0b85b0e5376a25084c-dirty
devteam
parents:
0
diff
changeset
|
8 from galaxy.tools.util.galaxyops import parse_cols_arg |
| 0 | 9 from galaxy.tools.util import maf_utilities |
|
1
412a03dc24a2
planemo upload commit fa9a2569c9ece4b03f622d0b85b0e5376a25084c-dirty
devteam
parents:
0
diff
changeset
|
10 from bx.intervals.io import NiceReaderWrapper |
| 0 | 11 import bx.align.maf |
| 12 import fileinput | |
| 13 import sys | |
| 14 | |
| 15 def stop_err(msg): | |
| 16 sys.stderr.write(msg) | |
| 17 sys.exit() | |
| 18 | |
| 19 | |
| 20 if len(sys.argv) < 3: | |
| 21 stop_err("Incorrect number of arguments.") | |
| 22 | |
| 23 inp_file = sys.argv[1] | |
| 24 out_file = sys.argv[2] | |
| 25 fout = open(out_file, 'w') | |
| 26 int_file = sys.argv[3] | |
| 27 if int_file != "None": #The user has specified an interval file | |
| 28 dbkey_i = sys.argv[4] | |
| 29 chr_col_i, start_col_i, end_col_i, strand_col_i = parse_cols_arg( sys.argv[5] ) | |
| 30 | |
| 31 | |
| 32 def rateEstimator(block): | |
| 33 global alignlen, mismatches | |
| 34 | |
| 35 src1 = block.components[0].src | |
| 36 sequence1 = block.components[0].text | |
| 37 start1 = block.components[0].start | |
| 38 end1 = block.components[0].end | |
| 39 len1 = int(end1)-int(start1) | |
| 40 len1_withgap = len(sequence1) | |
| 41 mismatch = 0.0 | |
| 42 | |
| 43 for seq in range (1, len(block.components)): | |
| 44 src2 = block.components[seq].src | |
| 45 sequence2 = block.components[seq].text | |
| 46 start2 = block.components[seq].start | |
| 47 end2 = block.components[seq].end | |
| 48 len2 = int(end2)-int(start2) | |
| 49 for nt in range(len1_withgap): | |
| 50 if sequence1[nt] not in '-#$^*?' and sequence2[nt] not in '-#$^*?': # Not a gap or masked character | |
| 51 if sequence1[nt].upper() != sequence2[nt].upper(): | |
| 52 mismatch += 1 | |
| 53 | |
| 54 if int_file == "None": | |
| 55 p = mismatch/min(len1, len2) | |
| 56 print >> fout, "%s\t%s\t%s\t%s\t%s\t%s\t%d\t%d\t%.4f" % ( src1, start1, end1, src2, start2, end2, min(len1, len2), mismatch, p ) | |
| 57 else: | |
| 58 mismatches += mismatch | |
| 59 alignlen += min(len1, len2) | |
| 60 | |
| 61 | |
| 62 def main(): | |
| 63 skipped = 0 | |
| 64 not_pairwise = 0 | |
| 65 | |
| 66 if int_file == "None": | |
| 67 try: | |
| 68 maf_reader = bx.align.maf.Reader( open(inp_file, 'r') ) | |
| 69 except: | |
| 70 stop_err("Your MAF file appears to be malformed.") | |
| 71 print >> fout, "#Seq1\tStart1\tEnd1\tSeq2\tStart2\tEnd2\tL\tN\tp" | |
| 72 for block in maf_reader: | |
| 73 if len(block.components) != 2: | |
| 74 not_pairwise += 1 | |
| 75 continue | |
| 76 try: | |
| 77 rateEstimator(block) | |
| 78 except: | |
| 79 skipped += 1 | |
| 80 else: | |
| 81 index, index_filename = maf_utilities.build_maf_index( inp_file, species = [dbkey_i] ) | |
| 82 if index is None: | |
| 83 print >> sys.stderr, "Your MAF file appears to be malformed." | |
| 84 sys.exit() | |
| 85 win = NiceReaderWrapper( fileinput.FileInput( int_file ), | |
| 86 chrom_col=chr_col_i, | |
| 87 start_col=start_col_i, | |
| 88 end_col=end_col_i, | |
| 89 strand_col=strand_col_i, | |
| 90 fix_strand=True) | |
| 91 species = None | |
| 92 mincols = 0 | |
| 93 global alignlen, mismatches | |
| 94 | |
| 95 for interval in win: | |
| 96 alignlen = 0 | |
| 97 mismatches = 0.0 | |
| 98 src = "%s.%s" % ( dbkey_i, interval.chrom ) | |
| 99 for block in maf_utilities.get_chopped_blocks_for_region( index, src, interval, species, mincols ): | |
| 100 if len(block.components) != 2: | |
| 101 not_pairwise += 1 | |
| 102 continue | |
| 103 try: | |
| 104 rateEstimator(block) | |
| 105 except: | |
| 106 skipped += 1 | |
| 107 if alignlen: | |
| 108 p = mismatches/alignlen | |
| 109 else: | |
| 110 p = 'NA' | |
| 111 interval.fields.append(str(alignlen)) | |
| 112 interval.fields.append(str(mismatches)) | |
| 113 interval.fields.append(str(p)) | |
| 114 print >> fout, "\t".join(interval.fields) | |
| 115 #num_blocks += 1 | |
| 116 | |
| 117 if not_pairwise: | |
| 118 print "Skipped %d non-pairwise blocks" % (not_pairwise) | |
| 119 if skipped: | |
| 120 print "Skipped %d blocks as invalid" % (skipped) | |
| 121 | |
| 122 | |
| 123 if __name__ == "__main__": | |
| 124 main() |