[Bioperl-l] Recommended way to download qual files from Genbank?
Phillip San Miguel
pmiguel at purdue.edu
Fri Jan 11 20:46:59 UTC 2008
Hi Malcolm,
Yes that works great!
Well, one caveat:
If you download both the fasta and the qual files:
ncbi_eutil -search db=nucleotide term=AC207960 -fetch rettype=fasta >
AC207960.fasta
ncbi_eutil -search db=nucleotide term=AC207960 -fetch rettype=qual >
AC207960.fasta.qual
The "primary IDs" don't match. The fasta comes out:
>gi|154937460|gb|AC207960.1|
and the qual comes out:
>AC207960.1
which seems to choke most programs that use seq and qual (eg
cross_match) because they want the primary IDs of the seq and qual files
to match.
Otherwise fine, though.
Thanks,
Phillip
Cook, Malcolm wrote:
> Phillip:
>
> Of course - mea culpa - here's the full monty....
>
> Indeed NCBI's eutils can do this:
>
>
>> ncbi_eutil -search db=nucleotide term=AC207960 -fetch rettype=qual >
>>
> AC207960.qual
>
> which uses my script (attached) to wrap NCBI's eutils.
>
> It depends upon NCBI_PowerScripting.pm disributed in PowerFiles_0707.zip
> by NCBI in their "Jul 24-27, 2007" course found at
> http://www.ncbi.nlm.nih.gov/Class/PowerTools/eutils/scripts.html
>
> I made a single edit to NCBI_PowerScripting.pm to 'select STDERR' at the
> very beginning so that trace messages are not printed on STDOUT, such as
> this echoed header:
> Retrieving 1 records from nucleotide...
> ... and footer:
> Received records 1 - 1.
> Wrote data to -.
>
> (otherwise they are interspersed with downloaded qual files)
>
> It also depends on recent version of GetOpt::Long.
>
> Hope it helps.
>
> Malcolm Cook
> Database Applications Manager - Bioinformatics
> Stowers Institute for Medical Research - Kansas City, Missouri
>
>
>
>> -----Original Message-----
>> From: Phillip San Miguel [mailto:pmiguel at purdue.edu]
>> Sent: Friday, January 11, 2008 1:33 PM
>> To: Cook, Malcolm
>> Cc: Chris Fields; bioperl-l
>> Subject: Re: [Bioperl-l] Recommended way to download qual
>> files from Genbank?
>>
>> Hi Malcolm,
>> Looks like your email was (inadvertantly?) redacted in
>> some way. (No attachment and last sentence truncated.) Would
>> it be possible to get a complete version so I can be sure I'm
>> following you?
>> Thanks,
>> Phillip
>>
>> Cook, Malcolm wrote:
>>
>>> Indeed eutil is capable of this
>>>
>>> The following use of my ncbi_eutil (attached) script yeilds what you
>>> want:
>>>
>>> ncbi_eutil -search db=nucleotide term=AC207960 -fetch
>>>
>> rettype=qual >
>>
>>> AC207960.qual
>>>
>>> It depends on the version of NCBI_PowerScripting.pm , such as is
>>> included in
>>>
>>> Malcolm Cook
>>> Database Applications Manager - Bioinformatics Stowers
>>>
>> Institute for
>>
>>> Medical Research - Kansas City, Missouri
>>>
>>>
>>>
>>>
>>>> -----Original Message-----
>>>> From: bioperl-l-bounces at lists.open-bio.org
>>>> [mailto:bioperl-l-bounces at lists.open-bio.org] On Behalf Of Chris
>>>> Fields
>>>> Sent: Friday, January 11, 2008 11:10 AM
>>>> To: Phillip San Miguel
>>>> Cc: bioperl-l
>>>> Subject: Re: [Bioperl-l] Recommended way to download qual
>>>>
>> files from
>>
>>>> Genbank?
>>>>
>>>> I don't think this is possible with the current setup for
>>>> Bio::DB::GenBank (which the script uses). We'll have to
>>>>
>> investigate
>>
>>>> whether it is possible to retrieve this data via NCBI's
>>>>
>> eutils; if so
>>
>>>> we can try adding it in. If you want you can submit this as an
>>>> enhancement request via bugzilla for tracking:
>>>>
>>>> http://bugzilla.open-bio.org/
>>>>
>>>> chris
>>>>
>>>> On Jan 11, 2008, at 10:22 AM, Phillip San Miguel wrote:
>>>>
>>>>
>>>>
>>>>> No problem getting sequence from genbank via a myriad of
>>>>>
>> methods.
>>
>>>>> But as the volume of non-finished sequence in genbank
>>>>>
>> increases the
>>
>>>>> importance of also obtaining quality values for a given sequence
>>>>> increases. Some records include quality values.
>>>>>
>>>>> I typically use bp_fetch.pl to grab a sequence from genbank:
>>>>>
>>>>> bp_fetch.pl -fmt fasta net::genbank:AC207960
>>>>>
>>>>> sends the fasta sequence to STDOUT. But that bp_fetch.pl wasn't
>>>>> designed to pull down quals evidently:
>>>>>
>>>>> bp_fetch.pl -fmt qual net::genbank:AC207960
>>>>>
>>>>> gives:
>>>>>
>>>>> ------------- EXCEPTION: Bio::Root::Exception -------------
>>>>> MSG: You must pass a Bio::Seq::Quality or a Bio::Seq::PrimaryQual
>>>>> object to write_seq() as a parameter named "source"
>>>>> STACK: Error::throw
>>>>> STACK: Bio::Root::Root::throw /usr/local/perl_5.8/lib/site_perl/
>>>>> 5.8.8/Bio/Root/Root.pm:359
>>>>> STACK: Bio::SeqIO::qual::write_seq
>>>>>
>>>>>
>>>> /usr/local/perl_5.8/lib/site_perl/
>>>>
>>>>
>>>>> 5.8.8/Bio/SeqIO/qual.pm:205
>>>>> STACK: /usr/local/perl/bin/bp_fetch.pl:313
>>>>> -----------------------------------------------------------
>>>>>
>>>>> (running under bioperl 1.5.2)
>>>>>
>>>>> The quality values for this accession are in genbank as these URLs
>>>>> demonstrate:
>>>>>
>>>>>
>>>>>
>>>>>
>> http://www.ncbi.nlm.nih.gov/entrez/viewer.fcgi?db=nuccore&id=15493746
>>
>>>> 0
>>>>
>>>>
>>>>>
>>>>>
>> http://www.ncbi.nlm.nih.gov/entrez/viewer.fcgi?db=nuccore&list_uids=1
>>
>>>> 5
>>>>
>>>>
>>>>> 4937460&dopt=fasta
>>>>>
>>>>>
>>>>>
>>>>>
>> http://www.ncbi.nlm.nih.gov/entrez/viewer.fcgi?db=nuccore&list_uids=1
>>
>>>> 5
>>>>
>>>>
>>>>> 4937460&dopt=qual
>>>>>
>>>>> What is the best way to pull down these qual values? They aren't
>>>>> present in "GenBank(Full)" format. They are present in an ASN.1
>>>>> format.
>>>>>
>>>>> Advice would be appreciated.
>>>>>
>>>>> --
>>>>> Phillip
>>>>> Purdue Genomics Core Facility
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> _______________________________________________
>>>>> Bioperl-l mailing list
>>>>> Bioperl-l at lists.open-bio.org
>>>>> http://lists.open-bio.org/mailman/listinfo/bioperl-l
>>>>>
>>>>>
>>>> Christopher Fields
>>>> Postdoctoral Researcher
>>>> Lab of Dr. Robert Switzer
>>>> Dept of Biochemistry
>>>> University of Illinois Urbana-Champaign
>>>>
>>>>
>>>>
>>>> _______________________________________________
>>>> Bioperl-l mailing list
>>>> Bioperl-l at lists.open-bio.org
>>>> http://lists.open-bio.org/mailman/listinfo/bioperl-l
>>>>
>>>>
>>>>
>>>
>>>
>>
>>
More information about the Bioperl-l
mailing list