[Bioperl-l] Bio::Search results
Sean O'Keeffe
limericksean at gmail.com
Fri Sep 2 05:31:42 EDT 2005
I have those four lines (all contained on separate lines) but it didn't work .
I ran "gawk '{if(/^\/\/$/) {print "//\n\nHMMER 2.3.2 (Oct 2003)"} else
{print}}' my_file > whatever" on the pfam file to include a HMMER
header after the // separator for each result.
It now works fine. I still don't understand why I need to include this
header to parse the file. Although from memory don't Blast reports all
include Blast header lines to separate next results?
By the way I think that its the mailer that seems to be altering the
text in these scripts (don't know where the "=3D" came from - it's not
in my original). I'm running a standard linux (suse 9.0) using nedit
as an editor.
Sean.
On 9/2/05, Barry Moore <bmoore at genetics.utah.edu> wrote:
> Sean I see that my mailer messed up the line formatting also. What I
> meant was you should have the top four lines like this:
>
> #Line 1: hmmpfam - search one or more sequences against HMM database
> #Line 2: HMMER 2.3.2(Oct 2003)
> #Line 3: Copyright (C) 1992-2003 HHMI/Washington University School of
> #Line 4: Medicine Freely distributed under the GNU General Public
> License (GPL)
>
> If you have that, and fix the hsp to hit error mentioned it works for
> me.
>
> Barry
>
> -----Original Message-----
> From: bioperl-l-bounces at portal.open-bio.org
> [mailto:bioperl-l-bounces at portal.open-bio.org] On Behalf Of Barry Moore
> Sent: Thursday, September 01, 2005 11:34 AM
> To: Bioperl List
> Subject: RE: [Bioperl-l] Bio::Search results
>
> Sean-
>
>
>
> The top of you input hmmpfam file has these four lines:
>
>
>
> hmmpfam - search one or more sequences against HMM database HMMER 2.3.2
> (Oct 2003) Copyright (C) 1992-2003 HHMI/Washington University School of
> Medicine Freely distributed under the GNU General Public License (GPL)
>
>
>
> When I copied that out of you e-mail and into emacs on my computer I got
> it all as one line like this:
>
>
>
> hmmpfam - search one or more sequences against HMM database HMMER 2.3.2
> (Oct 2003) Copyright (C) 1992-2003 HHMI/Washington University School of
> Medicine Freely distributed under the GNU General Public License (GPL)
>
>
>
> The parser needs to see a line beginning with HMMER to trigger setting
> an internal variable called _reporttype (among other things). Since
> this never happened with a single line header I got the same results you
> did. Your header lines look OK in the e-mail you sent me, but assuming
> you have the same problem with you header lines that I did, fixing that
> should lead you to your next problem...
>
>
>
> This line of your code has an error:
>
>
>
> next unless ($hsp->name =~ /^ig|^lrr|^fn3|^egf|^tsp|^psi/i);
>
>
>
> HSP objects don't have names but their associated hit objects do, so
> this makes your script work as intended:
>
>
>
> next unless ($hit->name =~ /^ig|^lrr|^fn3|^egf|^tsp|^psi/i);
>
>
>
> I wonder if the header line problems could be the result of cross
> platform cutting and pasting (or just cross platform file
> incompatibilities). I noticed also that in you script that you sent the
> '=' turned into '=3D'. What computer platforms do you use in your work?
>
>
>
> HTH
>
>
>
> Barry
>
>
>
> -----Original Message-----
>
> From: Sean O'Keeffe [mailto:limericksean at gmail.com]
>
> Sent: Tuesday, August 30, 2005 9:28 AM
>
> To: Barry Moore
>
> Cc: bioperl-l at portal.open-bio.org
>
> Subject: Re: [Bioperl-l] Bio::Search results
>
>
>
> Hi Barry, thanks for the reply. Below is a snippet of the file (I
> generated it with hmmpfam using the alignment flag set to -A 0, to
> remove alignments - this shouldn't affect the parsing of the file) :
>
>
>
> hmmpfam - search one or more sequences against HMM database HMMER 2.3.2
> (Oct 2003) Copyright (C) 1992-2003 HHMI/Washington University School of
> Medicine Freely distributed under the GNU General Public License (GPL)
>
> - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
>
> HMM file: /usr/local/lib/pfam-tm
>
> Sequence file: Mus_musculus.NCBIM34.jul.pep.fa-short
>
> - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
>
>
>
> Query sequence: ENSMUSP00000089702
>
> Accession: [none]
>
> Description: [none]
>
>
>
> Scores for sequence family classification (score includes all domains):
>
> Model Description Score E-value
> N
>
> -------- ----------- ----- -------
> ---
>
> [no hits above thresholds]
>
>
>
> Parsed for domains:
>
> Model Domain seq-f seq-t hmm-f hmm-t score E-value
>
> -------- ------- ----- ----- ----- ----- ----- -------
>
> [no hits above thresholds]
>
> //
>
>
>
> Query sequence: ENSMUSP00000089701
>
> Accession: [none]
>
> Description: [none]
>
>
>
> Scores for sequence family classification (score includes all domains):
>
> Model Description Score E-value
> N
>
> -------- ----------- ----- -------
> ---
>
> [no hits above thresholds]
>
>
>
> Parsed for domains:
>
> Model Domain seq-f seq-t hmm-f hmm-t score E-value
>
> -------- ------- ----- ----- ----- ----- ----- -------
>
> [no hits above thresholds]
>
> //
>
>
>
> Query sequence: ENSMUSP00000020094
>
> Accession: [none]
>
> Description: [none]
>
>
>
> Scores for sequence family classification (score includes all domains):
>
> Model Description Score E-value
> N
>
> -------- ----------- ----- -------
> ---
>
> LRR_1 Leucine Rich Repeat 55.0 3.7e-15
> 6
>
> LRRNT Leucine rich repeat N-terminal domain 30.1 1.1e-07
> 1
>
>
>
> Parsed for domains:
>
> Model Domain seq-f seq-t hmm-f hmm-t score E-value
>
> -------- ------- ----- ----- ----- ----- ----- -------
>
> LRRNT 1/1 117 142 .. 1 34 [] 30.1 1.1e-07
>
> LRR_1 1/6 168 191 .. 1 25 [] 14.7 0.0049
>
> LRR_1 2/6 192 210 .. 1 25 [] 9.1 0.2
>
> LRR_1 3/6 212 237 .. 1 25 [] 8.4 0.26
>
> LRR_1 4/6 238 257 .. 1 25 [] 10.3 0.1
>
> LRR_1 5/6 259 282 .. 1 25 [] 10.1 0.12
>
> LRR_1 6/6 290 314 .. 1 25 [] 2.4 2
>
> //
>
>
>
> Query sequence: ENSMUSP00000074175
>
> Accession: [none]
>
> Description: [none]
>
>
>
> Scores for sequence family classification (score includes all domains):
>
> Model Description Score E-value
> N
>
> -------- ----------- ----- -------
> ---
>
> CUB CUB domain 232.5 1.3e-68
> 2
>
> Trypsin Trypsin 206.0 1.2e-60
> 1
>
> Sushi Sushi domain (SCR repeat) 88.7 2.6e-25
> 2
>
> EGF_CA Calcium binding EGF domain 29.4 1.8e-07
> 1
>
>
>
> Parsed for domains:
>
> Model Domain seq-f seq-t hmm-f hmm-t score E-value
>
> -------- ------- ----- ----- ----- ----- ----- -------
>
> CUB 1/2 16 137 .. 1 116 [] 70.0 1.1e-19
>
> EGF_CA 1/1 141 188 .. 1 55 [] 29.4 1.8e-07
>
> CUB 2/2 192 301 .. 1 116 [] 162.5 1.5e-47
>
> Sushi 1/2 308 370 .. 1 62 [] 47.5 6.5e-13
>
> Sushi 2/2 375 446 .. 1 62 [] 41.2 5.1e-11
>
> Trypsin 1/1 463 698 .. 1 259 [] 206.0 1.2e-60
>
> //
>
>
>
>
>
> Cheers,
>
> Sean.
>
>
>
>
>
>
>
> On 8/30/05, Barry Moore <bmoore at genetics.utah.edu> wrote:
>
> > Sean,
>
> >
>
> > Don't see anything obviously wrong. If you want to send your input
>
> > file, I'll try to recreate the problem.
>
> >
>
> > Barry
>
> >
>
> > -----Original Message-----
>
> > From: bioperl-l-bounces at portal.open-bio.org
>
> > [mailto:bioperl-l-bounces at portal.open-bio.org] On Behalf Of Sean
>
> > O'Keeffe
>
> > Sent: Tuesday, August 30, 2005 2:55 AM
>
> > To: bioperl-l at portal.open-bio.org
>
> > Subject: [Bioperl-l] Bio::Search results
>
> >
>
> > Hi,
>
> > The following code snippet is something I use to extract information
>
> > from hmmer result files:
>
> >
>
> > use Bio::SearchIO;
>
> >
>
> > my $in =3D new Bio::SearchIO( -format =3D> 'hmmer', -file =3D>
>
> > $ARGV[0] ); while(my $result =3D $in->next_result) {
>
> > print $result->query_name(), "\n",$result->query_description(),"\n";
>
> > while (my $hit =3D $result->next_hit) {
>
> > while(my $hsp =3D $hit->next_domain) {
>
> > next unless ($hsp->name =3D~ /^ig|^lrr|^fn3|^egf|^tsp|^psi/i);
>
> > print $hsp->start(),"\t",$hsp->end(),"\t",$hsp->evalue(),"\n";
>
> > }
>
> > }
>
> > }
>
> >
>
> > The input file is generated by hmmpfam and is given at the command
>
> > line. I use it to scan for specific domain names e.g ig, fn3 lrr etc.
>
> > This code works for the first loop and then ends so I get the name and
>
>
> > description (no hsp values as their are none for this result):
>
> >
>
> > ENSMUSP00000065602=20
>
> > pep:novel supercontig::NT_085813:405:1510:-1 gene:ENSMUSG00000054059
>
> > transcript:ENSMUST00000066517
>
> >
>
> > My question is why does the loop end after one instance. Incidentally
>
> > the outputted name and description above are the last ones in the
>
> > hmmer file (maybe the file is read from the back??? - don't know if
>
> > this means anything). Any thoughts would be appreciated. Thanks,
>
> > Sean.
>
> >
>
> > _______________________________________________
>
> > Bioperl-l mailing list
>
> > Bioperl-l at portal.open-bio.org
>
> > http://portal.open-bio.org/mailman/listinfo/bioperl-l
>
> >
>
>
>
>
>
>
> _______________________________________________
> Bioperl-l mailing list
> Bioperl-l at portal.open-bio.org
> http://portal.open-bio.org/mailman/listinfo/bioperl-l
>
> _______________________________________________
> Bioperl-l mailing list
> Bioperl-l at portal.open-bio.org
> http://portal.open-bio.org/mailman/listinfo/bioperl-l
>
More information about the Bioperl-l
mailing list