Mercurial > repos > jjohnson > find_in_reference
annotate find_in_reference.py @ 2:c4fd2ea4f988
Add the option to test the reversed sequence and the DNA reverse complement of the sequence (ignored if the sequence cannot be interpreted as DNA)
| author | Jim Johnson <jj@umn.edu> |
|---|---|
| date | Thu, 13 Nov 2014 14:09:50 -0600 |
| parents | e83e0ce8fb68 |
| children | 2429b413d90a |
| rev | line source |
|---|---|
| 0 | 1 #!/usr/bin/env python |
| 2 """ | |
| 3 # | |
| 4 #------------------------------------------------------------------------------ | |
| 5 # University of Minnesota | |
| 6 # Copyright 2013, Regents of the University of Minnesota | |
| 7 #------------------------------------------------------------------------------ | |
| 8 # Author: | |
| 9 # | |
| 10 # James E Johnson | |
| 11 # | |
| 12 #------------------------------------------------------------------------------ | |
| 13 """ | |
| 14 | |
| 15 """ | |
| 16 Takes 2 tabular files as input: | |
| 17 1. The file to be filtered | |
| 18 2. The reference file | |
| 19 | |
| 20 The string value of selected column of the input file is searched for | |
| 21 in the string values of the selected column of the reference file. | |
| 22 | |
| 23 The intended purpose is to filter a peptide fasta file in tabular format | |
| 24 by whether those peptide sequences are found in a reference fasta file. | |
| 25 | |
| 26 """ | |
| 27 import sys,re,os.path | |
| 28 import tempfile | |
| 29 import optparse | |
| 30 from optparse import OptionParser | |
| 31 import logging | |
| 32 | |
| 33 | |
| 34 def __main__(): | |
| 35 #Parse Command Line | |
| 36 parser = optparse.OptionParser() | |
| 37 parser.add_option( '-i', '--input', dest='input', help='The input file to filter. (Otherwise read from stdin)' ) | |
| 38 parser.add_option( '-r', '--reference', dest='reference', help='The reference file to filter against' ) | |
| 39 parser.add_option( '-o', '--output', dest='output', help='The output file for input lines filtered by reference') | |
| 40 parser.add_option( '-f', '--filtered', dest='filtered', help='The output file for input lines not in the output') | |
| 41 parser.add_option('-c','--input_column', dest='input_column', default=None, help='The column for the value in the input file. (first column = 1, default to last column)') | |
| 42 parser.add_option('-C','--reference_column', dest='reference_column', default=None, help='The column for the value in the reference file. (first column = 1, default to last column)') | |
| 43 parser.add_option( '-I', '--case_insensitive', dest='ignore_case', action="store_true", default=False, help='case insensitive' ) | |
|
1
e83e0ce8fb68
Add option to reverse the search, find reference field in input field
Jim Johnson <jj@umn.edu>
parents:
0
diff
changeset
|
44 parser.add_option( '-R', '--reverse_find', dest='reverse_find', action="store_true", default=False, help='find the reference string in the input string' ) |
|
2
c4fd2ea4f988
Add the option to test the reversed sequence and the DNA reverse complement of the sequence (ignored if the sequence cannot be interpreted as DNA)
Jim Johnson <jj@umn.edu>
parents:
1
diff
changeset
|
45 parser.add_option( '-B', '--test_reverse', dest='test_reverse', action="store_true", default=False, help='Also search for reversed input string in reference' ) |
|
c4fd2ea4f988
Add the option to test the reversed sequence and the DNA reverse complement of the sequence (ignored if the sequence cannot be interpreted as DNA)
Jim Johnson <jj@umn.edu>
parents:
1
diff
changeset
|
46 parser.add_option( '-D', '--test_dna_reverse_complement', dest='test_reverse_comp', action="store_true", default=False, help='Also search for the DNA reverse complement of input string' ) |
| 0 | 47 parser.add_option( '-k', '--keep', dest='keep', action="store_true", default=False, help='' ) |
| 48 parser.add_option( '-a', '--annotation_columns', dest='annotation_columns', default=None, help='If string is found, add these columns from reference' ) | |
| 49 parser.add_option( '-s', '--annotation_separator', dest='annotation_separator', default=';', help='separator character between annotations from different lines' ) | |
| 50 parser.add_option( '-S', '--annotation_col_sep', dest='annotation_col_sep', default=',', help='separator character between annotation column from the same line' ) | |
| 51 parser.add_option( '-d', '--debug', dest='debug', action='store_true', default=False, help='Turn on wrapper debugging to stdout' ) | |
| 52 (options, args) = parser.parse_args() | |
|
2
c4fd2ea4f988
Add the option to test the reversed sequence and the DNA reverse complement of the sequence (ignored if the sequence cannot be interpreted as DNA)
Jim Johnson <jj@umn.edu>
parents:
1
diff
changeset
|
53 |
|
c4fd2ea4f988
Add the option to test the reversed sequence and the DNA reverse complement of the sequence (ignored if the sequence cannot be interpreted as DNA)
Jim Johnson <jj@umn.edu>
parents:
1
diff
changeset
|
54 revcompl = lambda x: ''.join([{'A':'T','C':'G','G':'C','T':'A','a':'t','c':'g','g':'c','t':'a','N':'N','n':'n'}[B] for B in x][::-1]) |
|
c4fd2ea4f988
Add the option to test the reversed sequence and the DNA reverse complement of the sequence (ignored if the sequence cannot be interpreted as DNA)
Jim Johnson <jj@umn.edu>
parents:
1
diff
changeset
|
55 def test_rcomplement(seq, target): |
|
c4fd2ea4f988
Add the option to test the reversed sequence and the DNA reverse complement of the sequence (ignored if the sequence cannot be interpreted as DNA)
Jim Johnson <jj@umn.edu>
parents:
1
diff
changeset
|
56 if options.test_reverse_comp: |
|
c4fd2ea4f988
Add the option to test the reversed sequence and the DNA reverse complement of the sequence (ignored if the sequence cannot be interpreted as DNA)
Jim Johnson <jj@umn.edu>
parents:
1
diff
changeset
|
57 try: |
|
c4fd2ea4f988
Add the option to test the reversed sequence and the DNA reverse complement of the sequence (ignored if the sequence cannot be interpreted as DNA)
Jim Johnson <jj@umn.edu>
parents:
1
diff
changeset
|
58 comp = revcompl(seq) |
|
c4fd2ea4f988
Add the option to test the reversed sequence and the DNA reverse complement of the sequence (ignored if the sequence cannot be interpreted as DNA)
Jim Johnson <jj@umn.edu>
parents:
1
diff
changeset
|
59 return comp in target |
|
c4fd2ea4f988
Add the option to test the reversed sequence and the DNA reverse complement of the sequence (ignored if the sequence cannot be interpreted as DNA)
Jim Johnson <jj@umn.edu>
parents:
1
diff
changeset
|
60 except: |
|
c4fd2ea4f988
Add the option to test the reversed sequence and the DNA reverse complement of the sequence (ignored if the sequence cannot be interpreted as DNA)
Jim Johnson <jj@umn.edu>
parents:
1
diff
changeset
|
61 pass |
|
c4fd2ea4f988
Add the option to test the reversed sequence and the DNA reverse complement of the sequence (ignored if the sequence cannot be interpreted as DNA)
Jim Johnson <jj@umn.edu>
parents:
1
diff
changeset
|
62 return False |
|
c4fd2ea4f988
Add the option to test the reversed sequence and the DNA reverse complement of the sequence (ignored if the sequence cannot be interpreted as DNA)
Jim Johnson <jj@umn.edu>
parents:
1
diff
changeset
|
63 |
|
c4fd2ea4f988
Add the option to test the reversed sequence and the DNA reverse complement of the sequence (ignored if the sequence cannot be interpreted as DNA)
Jim Johnson <jj@umn.edu>
parents:
1
diff
changeset
|
64 def test_reverse(seq,target): |
|
c4fd2ea4f988
Add the option to test the reversed sequence and the DNA reverse complement of the sequence (ignored if the sequence cannot be interpreted as DNA)
Jim Johnson <jj@umn.edu>
parents:
1
diff
changeset
|
65 return options.test_reverse and seq and seq[::-1] in target |
|
c4fd2ea4f988
Add the option to test the reversed sequence and the DNA reverse complement of the sequence (ignored if the sequence cannot be interpreted as DNA)
Jim Johnson <jj@umn.edu>
parents:
1
diff
changeset
|
66 |
| 0 | 67 # Input files |
| 68 if options.input != None: | |
| 69 try: | |
| 70 inputPath = os.path.abspath(options.input) | |
| 71 inputFile = open(inputPath, 'r') | |
| 72 except Exception, e: | |
| 73 print >> sys.stderr, "failed: %s" % e | |
| 74 exit(2) | |
| 75 else: | |
| 76 inputFile = sys.stdin | |
| 77 # Reference | |
| 78 if options.reference == None: | |
| 79 print >> sys.stderr, "failed: reference file is required" | |
| 80 exit(2) | |
| 81 # Output files | |
| 82 outFile = None | |
| 83 filteredFile = None | |
| 84 if options.filtered == None and options.output == None: | |
| 85 #write to stdout | |
| 86 outFile = sys.stdout | |
| 87 else: | |
| 88 if options.output != None: | |
| 89 try: | |
| 90 outPath = os.path.abspath(options.output) | |
| 91 outFile = open(outPath, 'w') | |
| 92 except Exception, e: | |
| 93 print >> sys.stderr, "failed: %s" % e | |
| 94 exit(3) | |
| 95 if options.filtered != None: | |
| 96 try: | |
| 97 filteredPath = os.path.abspath(options.filtered) | |
| 98 filteredFile = open(filteredPath, 'w') | |
| 99 except Exception, e: | |
| 100 print >> sys.stderr, "failed: %s" % e | |
| 101 exit(3) | |
| 102 incol = -1 | |
| 103 if options.input_column and options.input_column > 0: | |
| 104 incol = int(options.input_column)-1 | |
| 105 refcol = -1 | |
| 106 if options.reference_column and options.reference_column > 0: | |
| 107 refcol = int(options.reference_column)-1 | |
| 108 if options.annotation_columns: | |
| 109 annotate = True | |
| 110 annotation_columns = [int(x) - 1 for x in options.annotation_columns.split(',')] | |
| 111 else: | |
| 112 annotate = False | |
| 113 refFile = None | |
| 114 num_found = 0 | |
| 115 num_novel = 0 | |
| 116 for ln,line in enumerate(inputFile): | |
| 117 annotations = [] | |
| 118 try: | |
| 119 found = False | |
| 120 search_string = line.split('\t')[incol].rstrip('\r\n') | |
| 121 if options.ignore_case: | |
| 122 search_string = search_string.upper() | |
| 123 if options.debug: | |
| 124 print >> sys.stderr, "search: %s" % (search_string) | |
| 125 refFile = open(options.reference,'r') | |
| 126 for tn,fline in enumerate(refFile): | |
| 127 fields = fline.split('\t') | |
|
1
e83e0ce8fb68
Add option to reverse the search, find reference field in input field
Jim Johnson <jj@umn.edu>
parents:
0
diff
changeset
|
128 target_string = fields[refcol].rstrip('\r\n') |
| 0 | 129 if options.ignore_case: |
| 130 target_string = target_string.upper() | |
|
1
e83e0ce8fb68
Add option to reverse the search, find reference field in input field
Jim Johnson <jj@umn.edu>
parents:
0
diff
changeset
|
131 search = search_string if not options.reverse_find else target_string |
|
e83e0ce8fb68
Add option to reverse the search, find reference field in input field
Jim Johnson <jj@umn.edu>
parents:
0
diff
changeset
|
132 target = target_string if not options.reverse_find else search_string |
| 0 | 133 if options.debug: |
|
1
e83e0ce8fb68
Add option to reverse the search, find reference field in input field
Jim Johnson <jj@umn.edu>
parents:
0
diff
changeset
|
134 print >> sys.stderr, "in: %s %s %s" % (search,search in target,target) |
|
2
c4fd2ea4f988
Add the option to test the reversed sequence and the DNA reverse complement of the sequence (ignored if the sequence cannot be interpreted as DNA)
Jim Johnson <jj@umn.edu>
parents:
1
diff
changeset
|
135 if search in target or test_reverse(search,target) or test_rcomplement(search,target): |
| 0 | 136 found = True |
| 137 if annotate: | |
| 138 annotation = options.annotation_col_sep.join([fields[i] for i in annotation_columns]) | |
| 139 annotations.append(annotation) | |
| 140 else: | |
| 141 break | |
| 142 if found: | |
| 143 num_found += 1 | |
| 144 if annotate: | |
| 145 line = '%s\t%s\n' % (line.rstrip('\r\n'),options.annotation_separator.join(annotations)) | |
| 146 if options.keep == True: | |
| 147 if outFile: | |
| 148 outFile.write(line) | |
| 149 else: | |
| 150 if filteredFile: | |
| 151 filteredFile.write(line) | |
| 152 else: | |
| 153 num_novel += 1 | |
| 154 if options.keep == True: | |
| 155 if filteredFile: | |
| 156 filteredFile.write(line) | |
| 157 else: | |
| 158 if outFile: | |
| 159 outFile.write(line) | |
| 160 except Exception, e: | |
| 161 print >> sys.stderr, "failed: Error reading %s - %s" % (options.reference,e) | |
| 162 finally: | |
| 163 if refFile: | |
| 164 refFile.close() | |
| 165 print >> sys.stdout, "found: %d novel: %d" % (num_found,num_novel) | |
| 166 | |
| 167 if __name__ == "__main__" : __main__() | |
| 168 |