I am running hmmscan on a fasta file containing a proteome using Pfam-A.hmm. I would like to parse the output from this scan to collect a list of domains identified in my sequences, as well as list any novel domain architectures which may have been identified. Currently I am parsing the output file using a regex to identify and collect the domains listed above the inclusion threshold. Unfortunately, I am not sure how to go about identifying novel domain architectures. In the end, I'm hoping to be able to produce a few figures listing the percentage of different Pfam-A domains in the proteome, the percentage number of sequences with a single domain, two domains, three domains etc. and the percentage of domain architectures previously seen before and those which are novel.
Any advice would be greatly appreciated!
Thanks, James
•
written
by
James Ashmore •
2 answers
I would also suggest that you look at the combination of --tblout and --domtblout flags:
From 'hmmcan -h':
--tblout <f> : save parseable table of per-sequence hits to file S
--domtblout <f> : save parseable table of per-domain hits to file S
Even better, most likely, is using pfam_scan.pl, provided by the Pfam group. See: ftp://ftp.sanger.ac.uk/pub/databases/Pfam/Tools/README (about 1/2 way down)
The Bioperl Search::IO module provides parsing for HMMER2 and HMMER3 output. I can provide more details if required.
More details would be really helpful, thank you! A sample of how it works in a Python script would also be great.
Hi neilfws,
I managed to get the Biopython Hmmer parser working. I can print out the query sequence name, and a list of domain hits for that query. However, in one example I have two domains which are both above the inclusion threshold and greatly overlap. How would I go about extracting the domain which is most probable in my sequence? I originally thought about taking the lowest e-value however when there are two genuine domains in my sequence this would only extract one of them. Thanks!
I thought as much, thanks a lot :)
Hello again, my script is working as I wanted thanks to your input. My output is now a list of sequence names with their associated Pfam domains tab delimited and written as they occur across the sequence e.g.
Seq1 PYRIN \t DEAD
I was wondering if there was any high throughput method of analysing whether the domain architectures of these sequences is novel or not? The closest I have found is to use the Pfam online tool to search for sequences with and without specified domains. Thank you.
Yes exactly, it means a combination of domains (in any order or in a specific order) not identified previously. For example on the Pfam database they list the types of architectures where the DDE domain has been found (http://pfam.sanger.ac.uk//family/PF13359.1#tabview=tab1). If you scroll down the architectures on that page you can see that some of those specific architectures have only been found once (There is only 1 sequence with the following architecture).
Log in to answer this question.
-
Can I do a local BLASTp with a PFAM domain as query?
written by Gerardo •
-
Get full output with hmmscan in bio3d
written by lessismore
-
Why do the number of domains in these two orthologs differ? Why is hmmscan annotating one of these …
written by Dunois •
-
How to parse HMMsearch domtblout, resolve overlaps & return PDA?
written by Anand Rao
-
hmmscan using Pfam-A.hmm target database
written by shibl_a •
-
Filtering HmmScan Domain Results
written by hewaholax •
-
Count Hits to Each Unique Domain in Hmmscan Results w/ Python?
written by ethanabaker1 •
-
Network Diagram Of Protein Domain Architectures
written by Pappu
-
Literature Mining For Domains Involved In Regeneration
written by James Ashmore •
-
Search/Align Homologs With The Same Architecture With Hmmer
written by User 9501 •
