Mercurial > repos > jjohnson > find_in_reference
changeset 0:fe0327a3ba81
Uploaded
author | jjohnson |
---|---|
date | Sat, 04 Jan 2014 09:03:57 -0500 |
parents | |
children | 856033fb26e8 |
files | find_in_reference.py find_in_reference.xml test-data/found_peptides.tabular test-data/human_peptides.tabular test-data/human_proteins.tabular test-data/novel_peptides.tabular |
diffstat | 6 files changed, 211 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/find_in_reference.py Sat Jan 04 09:03:57 2014 -0500 @@ -0,0 +1,123 @@ +#!/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( '-k', '--keep', dest='keep', action="store_true", default=False, help='' ) + 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 + refFile = None + for ln,line in enumerate(inputFile): + try: + found = False + search_string = line.split('\t')[incol].rstrip('\r\n') + if options.debug: + print >> sys.stderr, "search: %s" % (search_string) + refFile = open(options.reference,'r') + for tn,fline in enumerate(refFile): + target_string = fline.split('\t')[refcol] + 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 + break + if found: + if options.keep == True: + if outFile: + outFile.write(line) + else: + if filteredFile: + filteredFile.write(line) + else: + 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() + +if __name__ == "__main__" : __main__() +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/find_in_reference.xml Sat Jan 04 09:03:57 2014 -0500 @@ -0,0 +1,68 @@ +<?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 + #if 'novel' in $outputs.__str__ or not 'found' in $outputs.__str__: + --output "$novel" + #end if + #if 'found' in $outputs.__str__: + --filtered "$found" + #end if + </command> + <inputs> + <param name="input" type="data" format="tabular" label="Input file to be filtered" + help="e.g. tophat junctions.bed run without GTF option or no-novel-junctions"/> + <param name="reference" type="data" format="tabular" label="reference file to search" + help="e.g. tophat junctions.bed run with GTF option and no-novel-junctions"/> + <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="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> + </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 +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. + + </help> +</tool>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/found_peptides.tabular Sat Jan 04 09:03:57 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 Sat Jan 04 09:03:57 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 Sat Jan 04 09:03:57 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 Sat Jan 04 09:03:57 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