I've detected homology between targets of ligands in drugbank and proteins in the proteome of a pathogen. I've parsed the output very rudimentary and calculated my query coverage. This exists in an excel file (see attached)

I'd like to further filter these hits. For query coverages > 70% alone I have 1500+ hits. Some redundancy exists and some hits are for drugbank targets for which no data has been deposited in the PDB.

So there are ways of further filtering these hits, but I'd like to do pairwise alignments first before I do anything else, as I think this is the most rational filter I can think of.

I do not want to filter based on lack of PDB data necessarily, as when I check my high query coverage hits, I'm getting predictions that are actually backed up by in vitro data. This includes high coverage hits lacking x-ray data

My question:

I have CLUSTALO installed on my desktop: How do programmatically send it jobs if I have the Drugbank sequences and the pathogen sequences in the same file? I would also like to parse out the pathogen IDs from the excel file. There’s 1500 and they are a bit scattered throughout the columns.

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

0.950704         -                    580         sp|P06672|PDC_ZYMMO                         -                     568                               1.2e-258             858.0          0.0            1        1      6.8e-262  1.5e-258   857.7      0.0       6         546    10     560    7      578    0.98   Uncharacterized                         protein                   OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_005014  PE=3             SV=1

0.950704         -                    691         sp|P06672|PDC_ZYMMO                         -                     568                               8.2e-162             538.3          0.1            1        1      2.4e-164  5.1e-161   535.6      0.1       1         541    89     650    89     658    0.89   Acetolactate                            synthase                  OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_006830  PE=3             SV=1

0.492958         -                    585         sp|P06672|PDC_ZYMMO                         -                     568                               8e-87                290.7          0.0            1        2      3.3e-83   7e-80      267.8      0.0       3         283    2      287    1      292    0.94   Uncharacterized                         protein                   OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_001897  PE=3             SV=1

0.107394         -                    585         sp|P06672|PDC_ZYMMO                         -                     568                               8e-87                290.7          0.0            2        2      1.7e-08   3.6e-05    21.1       0.0       413       474    436    505    376    583    0.72   Uncharacterized                         protein                   OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_001897  PE=3             SV=1
  • $\begingroup$ you will need to structure your questions so it is one question per post by breaking-down your goals into components. At present your question is not answerable I'm afraid, the site rules are designed to enable search engines to easily identify the post. $\endgroup$ – Michael Mar 25 at 10:04
  • $\begingroup$ Is this an automated comment? I shall parse this out into multiple questions (pun intended) $\endgroup$ – Mike Mar 25 at 10:24
  • $\begingroup$ Thanks @Mike (thats not automated "thanks") $\endgroup$ – Michael Mar 25 at 11:31
  • $\begingroup$ Also, make sure to tell us what operating system you are using, and why excel (!?) is involved. Finally, explain if you need every possible combination of pairwise alignments from your sequence file or only some. $\endgroup$ – terdon Mar 25 at 13:04
  • $\begingroup$ Thank you, Michael! Terdon, I just wanted to view the output in a consolidated table so I could send to others for review, though of course the awk one-liners and GREP commands were done on the command line (gnome emulator); I'm using the latest version of ubuntu $\endgroup$ – Mike Mar 25 at 13:17

The question is still in numerous parts,

  • parsing a blast output
  • obtaining sequences
  • feeding this as input into clustalo

If I was automating this I would parse the Blast genbank codes [regex either Perl one-liner (scripts aren't trendy ;-( ), or Pythons "re"] and feed them into NCBI's efetch to hook out the sequence data. Then I would create a subprocess that feed the files into clustalo. In Perl you simply issue the linux command wrapped in "`" either side, in Python use

import os
os.system('clustalo -in file.in -out file.out')

Ensuing your 'bin' is correctly configured to execute clustalo


import subprocess # use popen

To my knowledge there is no wrapper for clustalo unlike several other phylogenetics programs in BioPython, which enormously simplifies pipelines. Others might recommend using shell rather than e.g. Python.

BTW I personally use Muscle having been a long term user of clustalo, thats a long discussion.

To be honest you should look at the power of NBCI's online blast utility. The download facilities are certainly powerful particularly if you are pulling down e.g. 1500 sequences from blast hits, because it automatically can splice the hits of interest from their genes/genomes. However, its context specific.

|improve this answer|||||
  • $\begingroup$ To be clear, I already have the sequence data for these uniprot IDs on my desktop in another file and they are linked to their sequences. I wish to do pairwise alignments on sequences detected by jackhmmer $\endgroup$ – Mike Mar 26 at 8:30
  • $\begingroup$ Hmmmmm ... anyway subprocessing is the specific answer you requested on how to automate clustalo ... unless you are using the bash shell $\endgroup$ – Michael Mar 26 at 9:38
  • $\begingroup$ Currently doing some bipython tutorials at the moment, but yeah I do a lot through the command line $\endgroup$ – Mike Mar 26 at 10:34
  • $\begingroup$ Biopython does include a wrapper for ClustalO (and MAFFT, MUSCLE etc.) $\endgroup$ – Chris_Rands Mar 26 at 13:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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