gmane.comp.lang.perl.bio.general

Re: Speed issues with making IUPAC consensus from alignment

Fields, Christopher J — 2013-05-22T23:17:50


Probably the former, but...


It shouldn't take that long, 7 Mb isn't that large.  Or is that 7 Mb for one genome?


The code isn't really optimized for this, but again this isn't terribly large.  Is the bottleneck reading the alignment in, or is it the consensus_iupac() step?  Hard to say w/o seeing the alignment data itself.


chris

Speed issues with making IUPAC consensus from alignment

Senanu — 2013-05-22T20:15:24

Hi all,

I am wondering if the consensus_iupac method of Bio::Align is known to be extremely slow, or if I'm doing something wrong. 

I have bacterial whole-genome alignments (~7 Mbases) that I made in progressiveMauve and wish to get an IUPAC consensus. (I know that progressiveMauve uses a non-standard XMFA format, but Bio::AlignIO seems to read them just fine.) The code below takes more than all night to make a consensus. It works fine on tiny test alignments. 

Is this a known problem? Is there another way to generate such a consensus?


my $in = Bio::AlignIO->new(-file => $files[0],
                           -format => 'XMFA');
while  (my $aln = $in->next_aln()) {
    foreach  my $seq ($aln->each_seq) {
        $seq->alphabet('dna');
    }
    my $con = $aln->consensus_iupac();
}


Thanks in advance.
Ngwenyama

Phylip format error

Bernard Cohen — 2013-05-20T18:49:50

Hello!

I happen to have checked to see what the PERL webpage says about Phylip format for DNA alignment files and see that it is erroneous. 

I am not a PERL user and do not want to be bothered to register or otherwise learn how to make an official comment, so forward this for someone to pick up.

Phylip format allows up to 10 spaces for taxon names; the data must start in the 11th space. This can be checked on Jo Felsenstein's site.

The PERL page accessed by searching for "Phylip format PERL" allows only 8 spaces for the name. 

B. L. Cohen

Re: Problem compiling Bio::DB::Sam

Miquel Ràmia — 2013-05-21T15:08:18

Compiled correctly!

thank you

Re: Downloading annotated Gene nucleotide fasta files

Smithies, Russell — 2013-05-19T19:26:35

Another option that I've used before is to download the gene2accession, gene2refseq,  and gene_info files from here ftp://ftp.ncbi.nih.gov/gene/DATA/ then parse out the data you require.
It might work for you?

--Russell

-----Original Message-----
From: bioperl-l-bounces< at >lists.open-bio.org [mailto:bioperl-l-bounces< at >lists.open-bio.org] On Behalf Of shalabh sharma
Sent: Saturday, 18 May 2013 4:26 a.m.
To: Francisco J. Ossandón
Cc: bioperl-l
Subject: Re: [Bioperl-l] Downloading annotated Gene nucleotide fasta files

Hey Francisco,
           Thanks a lot. Basically i just wanted gene nucleotide fasta files with GI numbers.
I think i will have to parse it from gbk files.

-Shalabh


On Fri, May 17, 2013 at 11:59 AM, Francisco J. Ossandón < fossandonc< at >hotmail.com> wrote:



--
Shalabh Sharma
Scientific Computing Professional Associate (Bioinformatics Specialist) Department of Marine Sciences University of Georgia Athens, GA 30602-3636

_______________________________________________
Bioperl-l mailing list
Biop

Re: Downloading annotated Gene nucleotide fasta files

shalabh sharma — 2013-05-19T19:33:16

Thanks Russell,
                            Actually i wanted all the Bacterial gene
nucleotide files, so i parsed it from *gbk. But yes these files might help
me for my other parts of my work.

Thanks
Shalabh


On Sun, May 19, 2013 at 3:26 PM, Smithies, Russell <
Russell.Smithies< at >agresearch.co.nz> wrote:

Re: sets of sequences - how to read?

Fields, Christopher J — 2013-05-17T17:37:53

On May 17, 2013, at 12:12 AM, Carnë Draug <carandraug+dev< at >gmail.com>
 wrote:



(note: critical point in this is Bio::ASN1::Entrezgene would allow this, I'm not sure it would.  Otherwise this is all really hand-wavy)

To me a 'set of stuff', particularly when the 'stuff' is stored sequentially in a flat file, is a simple 'database' or 'store' of similar items, where the class allows one the ability to look up particular members in the set, but also could store higher level information about the set as a whole if needed.  If it were me, I would implement a method particular to Bio::SeqIO::entrezgene that specifically creates and returns this ( next_geneset(), for instance ); next_seq() could then be implemented to iterate through the items in that database/store.  

Two useful things come out of this.  First, if the data for the Entrez Gene file/chunk are parsed to store offsets per ID, one would only need to parse out the chunks needed (offset of ID to next offset), then pass that into the parser and create

Re: Downloading annotated Gene nucleotide fasta files

shalabh sharma — 2013-05-17T16:26:26

Hey Francisco,
           Thanks a lot. Basically i just wanted gene nucleotide fasta
files with GI numbers.
I think i will have to parse it from gbk files.

-Shalabh


On Fri, May 17, 2013 at 11:59 AM, Francisco J. Ossandón <
fossandonc< at >hotmail.com> wrote:

Re: Downloading annotated Gene nucleotide fasta files

Francisco J. Ossandón — 2013-05-17T15:59:04

Hi,
You can get the annotations from here:
ftp://ftp.ncbi.nih.gov/genomes/Bacteria/

The ".ffn" are the genes nucleotide fasta files but it does not show the
product name, on the other hand the ".faa" are the genes aminoacid fasta
files and shows the product name, but if you want both product and
nucleotide is much better to use the Genbank ".gbk" files that contains the
complete data and you can parse it easily using BioPerl to obtain all genes,
and then print the /protein_id, /product, and the nucleotide sequences in a
new fasta file. Check these to see how to do it:
http://www.bioperl.org/wiki/HOWTO:SeqIO
http://www.bioperl.org/wiki/HOWTO:Feature-Annotation

Cheers,

Francisco J. Ossandon

-----Mensaje original-----
De: bioperl-l-bounces< at >lists.open-bio.org
[mailto:bioperl-l-bounces< at >lists.open-bio.org] En nombre de shalabh sharma
Enviado el: viernes, 17 de mayo de 2013 10:55
Para: bioperl-l
Asunto: [Bioperl-l] Downloading annotated Gene nucleotide fasta files

HI,
                 First of all i am really

Downloading annotated Gene nucleotide fasta files

shalabh sharma — 2013-05-17T14:54:55

HI,
                 First of all i am really sorry for sending this mail here,
i am not sure if this is the right forum. I know lot of people work on
similar stuff.
I wrote to NCBI but nobody replied.

Actually i am looking for all bacterial/microbial gene annotation
nucleotide fasta files.
Does anyone knows where to download these kind of files.
I tried *ffn files but they are not annotated.
Or is there any module in bioperl that i can use ?
I would really appreciate your help.

Thanks
Shalabh

Re: sets of sequences - how to read?

Carnë Draug — 2013-05-17T05:12:24

:s I'm not sure I understood your suggestion. I think the problem is
just the introduction of a new concept, a "set" of stuff (genes in
this case), and how should SeqIO handle multiple sets.

Carnë

Re: Problem compiling Bio::DB::Sam

Fields, Christopher J — 2013-05-17T04:16:12

For the record, this is now fixed in the latest Bio::Samtools (via Lincoln).

chris

On May 16, 2013, at 8:12 AM, "Fields, Christopher J" <cjfields< at >illinois.edu> wrote:

Re: sets of sequences - how to read?

Fields, Christopher J — 2013-05-17T04:08:04

This doesn't surprise me too much; I know there have been some changes brewing, but didn't know when they would land.  I guess that would be... <looks at watch>... now.

My feeling is this will require writing some code for a higher-level layer of abstraction, say a Bio::DB::* (which would allow some internal indexing of the files maybe using a Bio::Index::*, look ups for specific gene IDs, etc).  How hard that would be to implement is another thing, have no idea w/o seeing what the data look like beyond they are in ASN1.

chris

On May 15, 2013, at 8:53 PM, Carnë Draug <carandraug+dev< at >gmail.com> wrote:

Re: [Biojava-l] Workshop on Sustainable Software for Science: Practice and Experiences

Andreas Prlic — 2013-05-16T04:31:34

Thanks Hilmar, you were faster than me in sending this out..

You are right, it would be very interesting to hear what some of the long
running open-bio projects have to say on the topic of sustainability. Let
me know if anybody is interested in a submission!

Andreas




On Wed, May 15, 2013 at 1:44 PM, Hilmar Lapp <hlapp< at >drycafe.net> wrote:

Re: Workshop on Sustainable Software for Science: Practice and Experiences

Fields, Christopher J — 2013-05-16T13:09:45

Yes, though we need to make sure others (e.g. those not subscribed to open-bio-l) are in the loop.

November is a possibility for me.

chris

On May 16, 2013, at 4:10 AM, Peter Cock <p.j.a.cock< at >googlemail.com>
 wrote:

Re: Problem compiling Bio::DB::Sam

Fields, Christopher J — 2013-05-16T13:12:40

It may be due to the new samtools release (v 0.1.19).  I know Heng Li has been working on the code over the last year for threading support (notice the undefined functions).  Have you tried v 0.1.18?

chris

On May 16, 2013, at 5:42 AM, Miquel Ràmia <miquel.ramia< at >uab.cat> wrote:

Problem compiling Bio::DB::Sam

Miquel Ràmia — 2013-05-16T10:42:29

Hi all,

I get this message when compiling Bio::DB::Sam:

Building Bio-SamTools

gcc -g -Wall -O2 -fPIC -o bam2bedgraph bam2bedgraph.o  -L/var/lib/gbrowse2/databases/samtools/samtools-0.1.19 -lbam -lm -lz

/var/lib/gbrowse2/databases/samtools/samtools-0.1.19/libbam.a(bgzf.o): In function `mt_destroy':

/var/lib/gbrowse2/databases/samtools/samtools-0.1.19/bgzf.c:458: undefined reference to `pthread_join'

/var/lib/gbrowse2/databases/samtools/samtools-0.1.19/libbam.a(bgzf.o): In function `bgzf_mt':

/var/lib/gbrowse2/databases/samtools/samtools-0.1.19/bgzf.c:445: undefined reference to `pthread_create'

collect2: ld returned 1 exit status

make: *** [bam2bedgraph] Error 1


Is this error related to the module or some dependencies? or maybe a 
problem with my system?

Any help appreciated, thanks!

Re: Workshop on Sustainable Software for Science: Practice and Experiences

Peter Cock — 2013-05-16T09:10:25

On Thu, May 16, 2013 at 5:43 AM, Fields, Christopher J
<cjfields< at >illinois.edu> wrote:

This sounds like a good idea, although given the time and place I am
unlikely to be able to attend in person:

First Workshop on Sustainable Software for Science: Practice and
Experiences (WSSSPE)
(to held in conjunction with SC13, Sunday, 17 November 2013, Denver, CO, USA)
http://wssspe.researchcomputing.org.uk/

Rather than trying to discuss this over four mailing lists should we switch
to the cross project list open-bio-l, or continue off-list?
http://lists.open-bio.org/mailman/listinfo/open-bio-l

Thanks,

Peter

Re: Workshop on Sustainable Software for Science: Practiceand Experiences

Fields, Christopher J — 2013-05-16T04:43:22

Jason and I have discussed looking into opportunity's like this, I think it makes sense to try a joint submission.  

chris

On May 15, 2013, at 3:44 PM, Hilmar Lapp <hlapp< at >drycafe.net> wrote:

sets of sequences - how to read?

Carnë Draug — 2013-05-16T01:53:55

Hi

when accessing entrez gene using eutils to get multiple genes, NCBI
now returns an Entrezgene-Set[1] rather than a list of EntrezGene.
This change must have happened sometime on the last 2 months. Compare:

use Bio::DB::EUtilities;

my %sets = (
  eutil   => 'efetch',
  db      => 'gene',
  retmode => 'text',
  rettype => 'asn1',
  email   => 'bioperl-l< at >lists.open-bio.org',
);

## this mimics the previous behaviour of the NCBI server but the
multiple requests will annoy their servers
my < at >ids = (3014, 85235);
my $response;
foreach (< at >ids) {
  my $fetcher = Bio::DB::EUtilities->new(%sets, id => $_);
  $response .= $fetcher->get_Response->content;
}
print $fetcher->get_Response->content;

## this used to be the right way to do it, but now returns an Entrezgene-Set
my $fetcher = Bio::DB::EUtilities->new(%sets, id => \< at >ids);
$response .= $fetcher->get_Response->content;
print $fetcher->get_Response->content;

There is no module to read these Entrezgene-Set in Perl at the moment,
since Bio::ASN1::EntrezGene; is

blastxml

Wenlan Tian — 2013-04-29T18:17:08

Hi,all,
I had done a blastall search and got a xml output file. How could i parse 
it to other format? I've found many scripts online, but no one worked. I am 
new to linux. Can anyone give me a good script to use?

Thanks,
Vivi