Mercurial > repos > jjohnson > find_in_reference

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/find_in_reference.py	Wed Feb 05 08:12:47 2014 -0500
@@ -0,0 +1,149 @@
+#!/usr/bin/env python
+"""
+#
+#------------------------------------------------------------------------------
+#                         University of Minnesota
+#         Copyright 2013, Regents of the University of Minnesota
+#------------------------------------------------------------------------------
+# Author:
+#
+#  James E Johnson
+#
+#------------------------------------------------------------------------------
+"""
+
+"""
+Takes 2 tabular files as input:
+  1. The file to be filtered
+  2. The reference file
+
+The string value of selected column of the input file is searched for
+in the string values of the selected column of the reference file.
+
+The intended purpose is to filter a peptide fasta file in tabular format
+by whether those peptide sequences are found in a reference fasta file.
+
+"""
+import sys,re,os.path
+import tempfile
+import optparse
+from optparse import OptionParser
+import logging
+
+
+def __main__():
+  #Parse Command Line
+  parser = optparse.OptionParser()
+  parser.add_option( '-i', '--input', dest='input', help='The input file to filter. (Otherwise read from stdin)' )
+  parser.add_option( '-r', '--reference', dest='reference', help='The reference file to filter against' )
+  parser.add_option( '-o', '--output', dest='output', help='The output file for input lines filtered by reference')
+  parser.add_option( '-f', '--filtered', dest='filtered', help='The output file for input lines not in the output')
+  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)')
+  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)')
+  parser.add_option( '-I', '--case_insensitive', dest='ignore_case', action="store_true", default=False, help='case insensitive' )
+  parser.add_option( '-k', '--keep', dest='keep', action="store_true", default=False, help='' )
+  parser.add_option( '-a', '--annotation_columns', dest='annotation_columns', default=None, help='If string is found, add these columns from reference' )
+  parser.add_option( '-s', '--annotation_separator', dest='annotation_separator', default=';', help='separator character between annotations from different lines' )
+  parser.add_option( '-S', '--annotation_col_sep', dest='annotation_col_sep', default=',', help='separator character between annotation column from the same line' )
+  parser.add_option( '-d', '--debug', dest='debug', action='store_true', default=False, help='Turn on wrapper debugging to stdout'  )
+  (options, args) = parser.parse_args()
+  # Input files
+  if options.input != None:
+    try:
+      inputPath = os.path.abspath(options.input)
+      inputFile = open(inputPath, 'r')
+    except Exception, e:
+      print >> sys.stderr, "failed: %s" % e
+      exit(2)
+  else:
+    inputFile = sys.stdin
+  # Reference
+  if options.reference == None:
+      print >> sys.stderr, "failed: reference file is required"
+      exit(2)
+  # Output files
+  outFile = None
+  filteredFile = None
+  if options.filtered == None and options.output == None:
+    #write to stdout
+    outFile = sys.stdout
+  else:
+    if options.output != None:
+      try:
+        outPath = os.path.abspath(options.output)
+        outFile = open(outPath, 'w')
+      except Exception, e:
+        print >> sys.stderr, "failed: %s" % e
+        exit(3)
+    if options.filtered != None:
+      try:
+        filteredPath = os.path.abspath(options.filtered)
+        filteredFile = open(filteredPath, 'w')
+      except Exception, e:
+        print >> sys.stderr, "failed: %s" % e
+        exit(3)
+  incol = -1
+  if options.input_column and options.input_column > 0:
+    incol = int(options.input_column)-1
+  refcol = -1
+  if options.reference_column and options.reference_column > 0:
+    refcol = int(options.reference_column)-1
+  if options.annotation_columns:
+    annotate = True
+    annotation_columns = [int(x) - 1 for x in options.annotation_columns.split(',')]
+  else:
+    annotate = False
+  refFile = None
+  num_found = 0
+  num_novel = 0
+  for ln,line in enumerate(inputFile):
+    annotations = []
+    try:
+      found = False
+      search_string = line.split('\t')[incol].rstrip('\r\n')
+      if options.ignore_case:
+        search_string = search_string.upper()
+      if options.debug:
+        print >> sys.stderr, "search: %s" % (search_string)
+      refFile = open(options.reference,'r')
+      for tn,fline in enumerate(refFile):
+        fields = fline.split('\t')
+        target_string =fields[refcol]
+        if options.ignore_case:
+          target_string = target_string.upper()
+        if options.debug:
+          print >> sys.stderr, "in: %s %s %s" % (search_string,search_string in target_string,target_string)
+        if search_string in target_string:
+          found = True
+          if annotate:
+            annotation = options.annotation_col_sep.join([fields[i] for i in annotation_columns])
+            annotations.append(annotation)
+          else:
+            break
+      if found:
+        num_found += 1
+        if annotate:
+          line = '%s\t%s\n' % (line.rstrip('\r\n'),options.annotation_separator.join(annotations))
+        if options.keep == True:
+          if outFile:
+            outFile.write(line)
+        else:
+          if filteredFile:
+            filteredFile.write(line)
+      else:
+        num_novel += 1
+        if options.keep == True:
+          if filteredFile:
+            filteredFile.write(line)
+        else:
+          if outFile:
+            outFile.write(line)
+    except Exception, e:
+      print >> sys.stderr, "failed: Error reading %s - %s" % (options.reference,e)
+    finally:
+      if refFile:
+        refFile.close()
+  print >> sys.stdout, "found: %d novel: %d" % (num_found,num_novel)
+
+if __name__ == "__main__" : __main__()
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/find_in_reference.xml	Wed Feb 05 08:12:47 2014 -0500
@@ -0,0 +1,137 @@
+<?xml version="1.0"?>
+<tool id="find_in_reference" name="find in reference" version="0.0.1">
+  <description>filter peptides that are present in proteins</description>
+  <command interpreter="python">find_in_reference.py  --input "$input"
+  --reference "$reference"
+  #if $column.set == 'yes':
+    --input_column $column.input_column
+    --reference_column $column.reference_column
+  #end if
+  $case_insensitive
+  #if 'novel' in $outputs.__str__ or not 'found' in $outputs.__str__:
+    --output "$novel"
+  #end if
+  #if 'found' in $outputs.__str__:
+    --filtered "$found"
+    #if $annotate.from_ref == 'yes' and str($annotate.annotation_columns) != 'None':
+      --annotation_columns $annotate.annotation_columns
+      #if $annotate.annotation_separator != '':
+        --annotation_separator '$annotate.annotation_separator'
+      #end if
+      #if $annotate.annotation_col_sep != '':
+        --annotation_col_sep '$annotate.annotation_col_sep'
+      #end if
+    #end if
+  #end if
+  </command>
+  <inputs>
+    <param name="input" type="data" format="tabular" label="Input file to be filtered"
+           help="e.g. a peptide fasta converted to tabular"/>
+    <param name="reference" type="data" format="tabular" label="reference file to search"
+           help="e.g. a protein fasta converted to tabular"/>
+    <conditional name="column">
+      <param name="set" type="select" label="select columns to compare">
+        <option value="no" selected="true">Use last column of input and reference</option>
+        <option value="yes">Choose the column of input and reference to compare</option>
+      </param>
+      <when value="no"/>
+      <when value="yes">
+        <param name="input_column" type="data_column" data_ref="input" label="column in input (defaults to last column)"
+           help=""/>
+        <param name="reference_column" type="data_column" data_ref="reference" label="column in reference (defaults to last column)"
+           help=""/>
+      </when>
+    </conditional>
+    <param name="case_insensitive" type="boolean" truevalue="--case_insensitive" falsevalue="" checked="false" label="Ignore case when comparing"/>
+    <param name="outputs" type="select" multiple="true" display="checkboxes" label="Choose outputs">
+      <option value="novel" selected="true">lines with no match in reference</option>
+      <option value="found">lines with match in reference</option>
+    </param>
+    <conditional name="annotate">
+      <param name="from_ref" type="select" label="Annotate found input entries with columns from reference">
+        <option value="no" selected="true">No</option>
+        <option value="yes">Yes</option>
+      </param>
+      <when value="no"/>
+      <when value="yes">
+        <param name="annotation_columns" type="data_column" data_ref="reference" multiple="true" label="columns from reference to append to found input lines"
+           help=""/>
+        <param name="annotation_separator" type="text" value=";" optional="true" label="separator to place between annotations from different reference lines"
+           help="defaults to ;">
+          <validator type="regex" message="Single quote character is not allowed">^[^']*$</validator>
+          <sanitizer>
+            <valid initial="string.printable">
+              <remove value="&apos;"/>
+            </valid>
+            <mapping initial="none">
+              <add source="&apos;" target=""/>
+            </mapping>
+          </sanitizer>
+        </param>
+        <param name="annotation_col_sep" type="text" value="," optional="true" label="separator to place between annotation columns from the same reference line"
+           help="defaults to ,">
+          <validator type="regex" message="Single quote character is not allowed">^[^']*$</validator>
+          <sanitizer>
+            <valid initial="string.printable">
+              <remove value="&apos;"/>
+            </valid>
+            <mapping initial="none">
+              <add source="&apos;" target=""/>
+            </mapping>
+          </sanitizer>
+        </param>
+      </when>
+    </conditional>
+  </inputs>
+  <stdio>
+    <exit_code range="1:" level="fatal" description="Error" />
+  </stdio>
+  <outputs>
+    <data name="found" metadata_source="input" format_source="input" label="${tool.name} on ${on_string}: found">
+      <filter>'found' in str(outputs)</filter>
+    </data>
+    <data name="novel" metadata_source="input" format_source="input" label="${tool.name} on ${on_string}: novel">
+      <filter>'novel' in str(outputs) or not 'found' in str(outputs)</filter>
+    </data>
+  </outputs>
+  <tests>
+    <test>
+      <param name="input" value="human_peptides.tabular" ftype="tabular" dbkey="hg19"/>
+      <param name="reference" value="human_proteins.tabular" ftype="tabular" dbkey="hg19"/>
+      <output name="novel" file="novel_peptides.tabular"/>
+    </test>
+  </tests>
+  <help>
+**Find in Reference**
+
+Filters lines of a tabular input file by checking if the selected input column value
+is a substring of the selected column of any line in the reference file.
+
+This can be used to check if peptides sequences are present in a set of reference proteins,
+as a means of filtering out uninteresting peptide sequences.
+
+For Example with::
+
+  Input
+    >pep1	LIL
+    >pep2	WTF
+    >pep3	ISK
+
+  Reference
+    >prot1	RLET
+    >prot2	LLIL
+    >prot3	LAPSE
+    >prot3	RISKY
+
+  The outputs
+
+  Not found in reference
+    >pep2	WTF
+
+  Found in reference
+    >pep1	LIL
+    >pep3	ISK
+
+
+  </help>
+</tool>
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/found_peptides.tabular	Wed Feb 05 08:12:47 2014 -0500
@@ -0,0 +1,6 @@
+pep_1|sp|Q9BS26|ERP44_HUMAN Endoplasmic reticulum resident protein 44 OS=Homo sapiens GN=ERP44 PE=1 SV=1	MHPAVFLSLPDLRCSLLLL
+pep_2|sp|Q9BS26|ERP44_HUMAN Endoplasmic reticulum resident protein 44 OS=Homo sapiens GN=ERP44 PE=1 SV=1	NENQVVFARVDCDQHSDIAQRYRISKY
+pep_3|sp|Q9BS26|ERP44_HUMAN Endoplasmic reticulum resident protein 44 OS=Homo sapiens GN=ERP44 PE=1 SV=1	FQKLAPSEYRYTLLRDRDEL
+pep_5|sp|P06213|INSR_HUMAN Insulin receptor OS=Homo sapiens GN=INSR PE=1 SV=4	INGQFVERCWT
+pep_7|sp|P08100|OPSD_HUMAN Rhodopsin OS=Homo sapiens GN=RHO PE=1 SV=1	MNGTEGPNFYVPFSNA
+pep_8|sp|P08100|OPSD_HUMAN Rhodopsin OS=Homo sapiens GN=RHO PE=1 SV=1	YNPVIYIMMNKQFRNCMLTTICCGKNPLGDDEASATVSKTETSQVAPA
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/human_peptides.tabular	Wed Feb 05 08:12:47 2014 -0500
@@ -0,0 +1,8 @@
+pep_1|sp|Q9BS26|ERP44_HUMAN Endoplasmic reticulum resident protein 44 OS=Homo sapiens GN=ERP44 PE=1 SV=1	MHPAVFLSLPDLRCSLLLL
+pep_2|sp|Q9BS26|ERP44_HUMAN Endoplasmic reticulum resident protein 44 OS=Homo sapiens GN=ERP44 PE=1 SV=1	NENQVVFARVDCDQHSDIAQRYRISKY
+pep_3|sp|Q9BS26|ERP44_HUMAN Endoplasmic reticulum resident protein 44 OS=Homo sapiens GN=ERP44 PE=1 SV=1	FQKLAPSEYRYTLLRDRDEL
+pep_4|sp|Q9BS26|ERP44_HUMAN Endoplasmic reticulum resident protein 44 OS=Homo sapiens GN=ERP44 PE=1 SV=1	LLILVTWVFTPVTTEITSLDTE
+pep_5|sp|P06213|INSR_HUMAN Insulin receptor OS=Homo sapiens GN=INSR PE=1 SV=4	INGQFVERCWT
+pep_6|sp|P06213|INSR_HUMAN Insulin receptor OS=Homo sapiens GN=INSR PE=1 SV=4	YALVSLSFFRKLRLIRLET
+pep_7|sp|P08100|OPSD_HUMAN Rhodopsin OS=Homo sapiens GN=RHO PE=1 SV=1	MNGTEGPNFYVPFSNA
+pep_8|sp|P08100|OPSD_HUMAN Rhodopsin OS=Homo sapiens GN=RHO PE=1 SV=1	YNPVIYIMMNKQFRNCMLTTICCGKNPLGDDEASATVSKTETSQVAPA
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/human_proteins.tabular	Wed Feb 05 08:12:47 2014 -0500
@@ -0,0 +1,4 @@
+sp|Q9BS26|ERP44_HUMAN Endoplasmic reticulum resident protein 44 OS=Homo sapiens GN=ERP44 PE=1 SV=1	MHPAVFLSLPDLRCSLLLLVTWVFTPVTTEITSLDTENIDEILNNADVALVNFYADWCRFSQMLHPIFEEASDVIKEEFPNENQVVFARVDCDQHSDIAQRYRISKYPTLKLFRNGMMMKREYRGQRSVKALADYIRQQKSDPIQEIRDLAEITTLDRSKRNIIGYFEQKDSDNYRVFERVANILHDDCAFLSAFGDVSKPERYSGDNIIYKPPGHSAPDMVYLGAMTNFDVTYNWIQDKCVPLVREITFENGEELTEEGLPFLILFHMKEDTESLEIFQNEVARQLISEKGTINFLHADCDKFRHPLLHIQKTPADCPVIAIDSFRHMYVFGDFKDVLIPGKLKQFVFDLHSGKLHREFHHGPDPTDTAPGEQAQDVASSPPESSFQKLAPSEYRYTLLRDRDEL
+sp|Q9NSY1|BMP2K_HUMAN BMP-2-inducible protein kinase OS=Homo sapiens GN=BMP2K PE=1 SV=2	MKKFSRMPKSEGGSGGGAAGGGAGGAGAGAGCGSGGSSVGVRVFAVGRHQVTLEESLAEGGFSTVFLVRTHGGIRCALKRMYVNNMPDLNVCKREITIMKELSGHKNIVGYLDCAVNSISDNVWEVLILMEYCRAGQVVNQMNKKLQTGFTEPEVLQIFCDTCEAVARLHQCKTPIIHRDLKVENILLNDGGNYVLCDFGSATNKFLNPQKDGVNVVEEEIKKYTTLSYRAPEMINLYGGKPITTKADIWALGCLLYKLCFFTLPFGESQVAICDGNFTIPDNSRYSRNIHCLIRFMLEPDPEHRPDIFQVSYFAFKFAKKDCPVSNINNSSIPSALPEPMTASEAAARKSQIKARITDTIGPTETSIAPRQRPKANSATTATPSVLTIQSSATPVKVLAPGEFGNHRPKGALRPGNGPEILLGQGPPQQPPQQHRVLQQLQQGDWRLQQLHLQHRHPHQQQQQQQQQQQQQQQQQQQQQQQQQQQHHHHHHHHLLQDAYMQQYQHATQQQQMLQQQFLMHSVYQPQPSASQYPTMMPQYQQAFFQQQMLAQHQPSQQQASPEYLTSPQEFSPALVSYTSSLPAQVGTIMDSSYSANRSVADKEAIANFTNQKNISNPPDMSGWNPFGEDNFSKLTEEELLDREFDLLRSNRLEERASSDKNVDSLSAPHNHPPEDPFGSVPFISHSGSPEKKAEHSSINQENGTANPIKNGKTSPASKDQRTGKKTSVQGQVQKGNDESESDFESDPPSPKSSEEEEQDDEEVLQGEQGDFNDDDTEPENLGHRPLLMDSEDEEEEEKHSSDSDYEQAKAKYSDMSSVYRDRSGSGPTQDLNTILLTSAQLSSDVAVETPKQEFDVFGAVPFFAVRAQQPQQEKNEKNLPQHRFPAAGLEQEEFDVFTKAPFSKKVNVQECHAVGPEAHTIPGYPKSVDVFGSTPFQPFLTSTSKSESNEDLFGLVPFDEITGSQQQKVKQRSLQKLSSRQRRTKQDMSKSNGKRHHGTPTSTKKTLKPTYRTPERARRHKKVGRRDSQSSNEFLTISDSKENISVALTDGKDRGNVLQPEESLLDPFGAKPFHSPDLSWHPPHQGLSDIRADHNTVLPGRPRQNSLHGSFHSADVLKMDDFGAVPFTELVVQSITPHQSQQSQPVELDPFGAAPFPSKQ
+sp|P06213|INSR_HUMAN Insulin receptor OS=Homo sapiens GN=INSR PE=1 SV=4	MATGGRRGAAAAPLLVAVAALLLGAAGHLYPGEVCPGMDIRNNLTRLHELENCSVIEGHLQILLMFKTRPEDFRDLSFPKLIMITDYLLLFRVYGLESLKDLFPNLTVIRGSRLFFNYALVIFEMVHLKELGLYNLMNITRGSVRIEKNNELCYLATIDWSRILDSVEDNYIVLNKDDNEECGDICPGTAKGKTNCPATVINGQFVERCWTHSHCQKVCPTICKSHGCTAEGLCCHSECLGNCSQPDDPTKCVACRNFYLDGRCVETCPPPYYHFQDWRCVNFSFCQDLHHKCKNSRRQGCHQYVIHNNKCIPECPSGYTMNSSNLLCTPCLGPCPKVCHLLEGEKTIDSVTSAQELRGCTVINGSLIINIRGGNNLAAELEANLGLIEEISGYLKIRRSYALVSLSFFRKLRLIRGETLEIGNYSFYALDNQNLRQLWDWSKHNLTITQGKLFFHYNPKLCLSEIHKMEEVSGTKGRQERNDIALKTNGDQASCENELLKFSYIRTSFDKILLRWEPYWPPDFRDLLGFMLFYKEAPYQNVTEFDGQDACGSNSWTVVDIDPPLRSNDPKSQNHPGWLMRGLKPWTQYAIFVKTLVTFSDERRTYGAKSDIIYVQTDATNPSVPLDPISVSNSSSQIILKWKPPSDPNGNITHYLVFWERQAEDSELFELDYCLKGLKLPSRTWSPPFESEDSQKHNQSEYEDSAGECCSCPKTDSQILKELEESSFRKTFEDYLHNVVFVPRKTSSGTGAEDPRPSRKRRSLGDVGNVTVAVPTVAAFPNTSSTSVPTSPEEHRPFEKVVNKESLVISGLRHFTGYRIELQACNQDTPEERCSVAAYVSARTMPEAKADDIVGPVTHEIFENNVVHLMWQEPKEPNGLIVLYEVSYRRYGDEELHLCVSRKHFALERGCRLRGLSPGNYSVRIRATSLAGNGSWTEPTYFYVTDYLDVPSNIAKIIIGPLIFVFLFSVVIGSIYLFLRKRQPDGPLGPLYASSNPEYLSASDVFPCSVYVPDEWEVSREKITLLRELGQGSFGMVYEGNARDIIKGEAETRVAVKTVNESASLRERIEFLNEASVMKGFTCHHVVRLLGVVSKGQPTLVVMELMAHGDLKSYLRSLRPEAENNPGRPPPTLQEMIQMAAEIADGMAYLNAKKFVHRDLAARNCMVAHDFTVKIGDFGMTRDIYETDYYRKGGKGLLPVRWMAPESLKDGVFTTSSDMWSFGVVLWEITSLAEQPYQGLSNEQVLKFVMDGGYLDQPDNCPERVTDLMRMCWQFNPKMRPTFLEIVNLLKDDLHPSFPEVSFFHSEENKAPESEELEMEFEDMENVPLDRSSHCQREEAGGRDGGSSLGFKRSYEEHIPYTHMNGGKKNGRILTLPRSNPS
+sp|P08100|OPSD_HUMAN Rhodopsin OS=Homo sapiens GN=RHO PE=1 SV=1	MNGTEGPNFYVPFSNATGVVRSPFEYPQYYLAEPWQFSMLAAYMFLLIVLGFPINFLTLYVTVQHKKLRTPLNYILLNLAVADLFMVLGGFTSTLYTSLHGYFVFGPTGCNLEGFFATLGGEIALWSLVVLAIERYVVVCKPMSNFRFGENHAIMGVAFTWVMALACAAPPLAGWSRYIPEGLQCSCGIDYYTLKPEVNNESFVIYMFVVHFTIPMIIIFFCYGQLVFTVKEAAAQQQESATTQKAEKEVTRMVIIMVIAFLICWVPYASVAFYIFTHQGSNFGPIFMTIPAFFAKSAAIYNPVIYIMMNKQFRNCMLTTICCGKNPLGDDEASATVSKTETSQVAPA
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/novel_peptides.tabular	Wed Feb 05 08:12:47 2014 -0500
@@ -0,0 +1,2 @@
+pep_4|sp|Q9BS26|ERP44_HUMAN Endoplasmic reticulum resident protein 44 OS=Homo sapiens GN=ERP44 PE=1 SV=1	LLILVTWVFTPVTTEITSLDTE
+pep_6|sp|P06213|INSR_HUMAN Insulin receptor OS=Homo sapiens GN=INSR PE=1 SV=4	YALVSLSFFRKLRLIRLET