1
$\begingroup$ 
0.980634         -                    572         sp|P06672|PDC_ZYMMO                         -                     568                               1.2e-263             874.5          0.0            1        1      6.6e-267  1.4e-263   874.2      0.0       5         562    15     570    12     572    0.98   Uncharacterized                         protein                   OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_005678  PE=3             SV=1

I've asked a similar question before. I have a better understanding of python now and would like to revisit this. Can someone critique my thinking?

I have an HMM output file that contains about 20 fields. There are two sets of IDs, e.g., P006672 and jhhlp_005678

In a separate file are FASTA-formatted sequences linked to these IDs.

Proposed solution:

1

   a) read in fields containing IDs, lets called them fields 8 and 6; 


   b) store field values in an array or dictionary (whatever the proper data structure is for 
  this). Let's called this Fields_8_and_6

  1. a) read proteome file containing FASTA-formatted sequences

    b) Grep and store as ID:sequence key:value pairs; let's call this ID_to_sequence

  2. a) Compare Fields_8_and_6 to ID_to_sequence

    b) Write ID1:sequence1 ID2:sequence2 as pairs into one file

  1. send to ClustalO

Improve this question
$\endgroup$
3
  • $\begingroup$ I'm not sure if this will be helpful to you, but my initial idea would be to use the very helpful pandas module in python to parse the table, then use some sort of regex to recover the required IDs (for example, I might set a regex query to recover 'text between vertical lines' to get P06672, and 'text after =' to get jhhlp_005678). There are tutorials for pandas and regex on the Software Carpentry website. $\endgroup$
    – Laura
    Commented Aug 18, 2020 at 9:49
  • $\begingroup$ Got it, I know I've used both biopython and regex to parse the IDs before, and I've used awk to calculate query coverage (e.g., awk '{$1 = ($17-$16)/$6}1' ). Thanks very much for this $\endgroup$
    – Mike
    Commented Aug 18, 2020 at 10:03
  • $\begingroup$ pandas will do it, but as pointed out BioPython has an HMM package. $\endgroup$
    – M__ ♦
    Commented Apr 5 at 4:03

1 Answer 1

Reset to default
3
$\begingroup$

Biopython can index a fastafile and you can access the sequences by their ID in a dict-like manner.

All that is left for you to do is extract the IDs from your file. (more Biopython)

Improve this answer
$\endgroup$
4
  • $\begingroup$ I didn't know they supported parsing hmmer outputs, thank you very much $\endgroup$
    – Mike
    Commented Aug 18, 2020 at 10:06
  • $\begingroup$ Hi Pallie, I merged my drugbank file with my proteome file; I'm getting duplicated because of the drugbank file. Is there a way around this with additional keys or something of that sort? $\endgroup$
    – Mike
    Commented Aug 24, 2020 at 11:34
  • $\begingroup$ Hi Mike, open a new question and show your code, input, and expected output so people can help :) $\endgroup$
    – Pallie
    Commented Aug 24, 2020 at 12:22
  • $\begingroup$ sounds great thank you! This is a really powerful protocol, automating it and getting it on the web would go a long way $\endgroup$
    – Mike
    Commented Aug 24, 2020 at 14:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.