<div dir="ltr"><div><div>A year or two ago I have already written a numpy-based implementation of distance matrix precisely because biopython's one was both slow and inconvenient. It's available at <a href="https://github.com/SynedraAcus/sampler/blob/master/reductor.py">https://github.com/SynedraAcus/sampler/blob/master/reductor.py</a> (see `DistanceMatrix` class)<br><br></div><div>What's fixed:<br></div><div>- Sequence IDs can be given in any order. Those `dm[j,i] if j>=i else dm[j,i]` calls are annoying, and that gets even uglier (and damn slow, see next point) with str sequence IDs.<br></div><div>- Dict, not list, for the sequence ID/index mapping. The current implementation stores names in a list and iterates over it 4 (four!) times every single time __getitem__ is called with the string IDs.<br></div><div>- Numpy storage for a reduced memory consumption<br></div><div>- A method for creating submatrices given a list of sequence IDs (probably pretty situational)<br></div><div>- Matrix factory, a separate class, takes FASTA file, not an alignment (very situational, needs to be expanded)<br></div><div><br></div><div>It's not compatible with other Biopython modules that use distances (ie Phylo, distance calculations), so I didn't send a pull request. If you ever decide to overhaul this part of project, please feel free to use my code. It's been licensed under CC-BY 4.0, but I can relicense it if you want.<br><br>The overhaul actually is necessary: it's all pretty slow with more than a few tens of sequences and the API sometimes is counterintuitive. Maybe a MatrixIO class and matrix subclasses would be necessary. But I personally haven't even touched other modules and can't promise I'll get to it any time soon.<br></div></div></div><div class="gmail_extra"><br><div class="gmail_quote">2017-08-16 4:31 GMT+08:00 Eric Talevich <span dir="ltr"><<a href="mailto:eric.talevich@gmail.com" target="_blank">eric.talevich@gmail.com</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><div><div><div><div>Hi Andrew,<br><br></div>If you have a DistanceMatrix object handy, you can export it to Phylip's matrix format with the method "format_phylip":<br><br></div><div>dm = DistanceMatrix(...)<br></div>with open("my_file.phylip", "w") as f:<br></div> dm.format_phyilp(f)<br><br></div>Then you can feed this file to another, faster program like Fasttree or Phylip's "neighbor" program.<br><br></div><div>This module doesn't use Numpy internally because it was written at a time when Biopython avoided having Numpy as a dependency. It would probably be a useful improvement to port it to Numpy or Scipy now.<span class="HOEnZb"><font color="#888888"><br><br></font></span></div><span class="HOEnZb"><font color="#888888"><div>-Eric<br></div></font></span></div><div class="HOEnZb"><div class="h5"><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Aug 15, 2017 at 9:26 AM, Andrew Sanchez <span dir="ltr"><<a href="mailto:aas229@nau.edu" target="_blank">aas229@nau.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Here are the results. The run time scales quickly, indeed! It would take days to generate a tree for a 6000 x 6000 matrix. I should have remembered this from my introduction to bioinformatics course.<br>
<br>
In [35]: for num in range(100, 1001, 100):<br>
...: print('Matrix size: %s' % num)<br>
...: %time gbf.estimate_tree_constructor_<wbr>runtime(num)<br>
...:<br>
Matrix size: 100<br>
Wall time: 2.36 s<br>
Matrix size: 200<br>
Wall time: 18.2 s<br>
Matrix size: 300<br>
Wall time: 1min 1s<br>
Matrix size: 400<br>
Wall time: 2min 25s<br>
Matrix size: 500<br>
Wall time: 4min 43s<br>
Matrix size: 600<br>
Wall time: 8min 8s<br>
Matrix size: 700<br>
Wall time: 13min 1s<br>
Matrix size: 800<br>
Wall time: 19min 15s<br>
Matrix size: 900<br>
Wall time: 27min 32s<br>
Matrix size: 1000<br>
Wall time: 37min 50s<br>
<span><br>
> It may be you would be better off with a compiled command line<br>
> phylogenetic tree for work at this scale<br>
<br>
</span>Thanks for the tip! Would you kindly point me in the right direction to begin understand how I can accomplish this?<br>
<br>
Thank you,<br>
Andrew<br>
<div class="m_-3771993471756367316HOEnZb"><div class="m_-3771993471756367316h5"><br>
> On Aug 14, 2017, at 5:31 PM, Peter Cock <<a href="mailto:p.j.a.cock@googlemail.com" target="_blank">p.j.a.cock@googlemail.com</a>> wrote:<br>
><br>
> Hi Andrew,<br>
><br>
> My guess is you are simply seeing the quadratic scaling being a<br>
> problem. Can you try timing a series of subsets, say 10 entries, 50,<br>
> 100, 200, 250, 500, 1000 - that approach ought to be enough to<br>
> estimate how long the full 6000 or so would take.<br>
><br>
> It may be you would be better off with a compiled command line<br>
> phylogenetic tree for work at this scale?<br>
><br>
> Peter<br>
><br>
> On Mon, Aug 14, 2017 at 4:21 PM, Andrew Sanchez <<a href="mailto:aas229@nau.edu" target="_blank">aas229@nau.edu</a>> wrote:<br>
>> I am trying to construct a tree from a DistanceMatrix object with len of 6303 with the following command: `tree = constructor.nj(bio_dmx)`.<br>
>><br>
>> The matrix and constructor were derived like so:<br>
>><br>
>> bio_dmx = _DistanceMatrix(names, nested_dmx)<br>
>> constructor = DistanceTreeConstructor()<br>
>><br>
>> I've tested my workflow on a much smaller distance matrix, just following the examples at <a href="http://biopython.org/wiki/Phylo" rel="noreferrer" target="_blank">http://biopython.org/wiki/Phyl<wbr>o</a> and it worked just fine. When I try to do it with this larger dataset, the process just hangs. I don't know where to begin debugging. First of all, how long should I expect this process to take? From wikipedia: “...typical run times proportional to approximately the square of the number of taxa."<br>
>><br>
>> Maybe it is normal for a tree of this size to take so long to construct? If so, is there a way to run tree = constructor.nj(bio_dmx) so that it produces some output that will allow me to at least see that something is happening?<br>
>><br>
>> I was trying to do this in an IPython session, and eventually I just cancelled the process which had been going for about 48 hours. The result of the keyboard interrupt was:<br>
>><br>
>> /home/aas229/anaconda3/envs/gb<wbr>filter/lib/python3.4/site-pack<wbr>ages/Bio/Phylo/TreeConstructio<wbr>n.py in nj(self, distance_matrix)<br>
>> 697 node_dist[i] = 0<br>
>> 698 for j in range(0, len(dm)):<br>
>> --> 699 node_dist[i] += dm[i, j]<br>
>> 700 node_dist[i] = node_dist[i] / (len(dm) - 2)<br>
>> 701<br>
>><br>
>> /home/aas229/anaconda3/envs/gb<wbr>filter/lib/python3.4/site-pack<wbr>ages/Bio/Phylo/TreeConstructio<wbr>n.py in __getitem__(self, item)<br>
>> 166 raise TypeError("Invalid index type.")<br>
>> 167 # check index<br>
>> --> 168 if row_index > len(self) - 1 or col_index > len(self) - 1:<br>
>> 169 raise IndexError("Index out of range.")<br>
>> 170 if row_index > col_index:<br>
>><br>
>> /home/aas229/anaconda3/envs/gb<wbr>filter/lib/python3.4/site-pack<wbr>ages/Bio/Phylo/TreeConstructio<wbr>n.py in __len__(self)<br>
>> 284 def __len__(self):<br>
>> 285 """Matrix length"""<br>
>> --> 286 return len(self.names)<br>
>> 287<br>
>> 288 def __repr__(self):<br>
>><br>
>> Does this output suggest that the job was in fact running just fine, but just taking a really long time?<br>
>><br>
>> Is there any other info that would be helpful in figuring this out?<br>
>><br>
>> Thank you,<br>
>> Andrew<br>
>> ______________________________<wbr>_________________<br>
>> Biopython mailing list - <a href="mailto:Biopython@mailman.open-bio.org" target="_blank">Biopython@mailman.open-bio.org</a><br>
>> <a href="http://mailman.open-bio.org/mailman/listinfo/biopython" rel="noreferrer" target="_blank">http://mailman.open-bio.org/ma<wbr>ilman/listinfo/biopython</a><br>
<br>
<br>
______________________________<wbr>_________________<br>
Biopython mailing list - <a href="mailto:Biopython@mailman.open-bio.org" target="_blank">Biopython@mailman.open-bio.org</a><br>
<a href="http://mailman.open-bio.org/mailman/listinfo/biopython" rel="noreferrer" target="_blank">http://mailman.open-bio.org/ma<wbr>ilman/listinfo/biopython</a></div></div></blockquote></div><br></div>
</div></div><br>______________________________<wbr>_________________<br>
Biopython mailing list - <a href="mailto:Biopython@mailman.open-bio.org">Biopython@mailman.open-bio.org</a><br>
<a href="http://mailman.open-bio.org/mailman/listinfo/biopython" rel="noreferrer" target="_blank">http://mailman.open-bio.org/<wbr>mailman/listinfo/biopython</a><br></blockquote></div><br><br clear="all"><br>-- <br><div class="gmail_signature" data-smartmail="gmail_signature"><font face="arial,helvetica,sans-serif">Alexey Morozov,<br>LIN SB RAS, bioinformatics group.<br>Irkutsk, Russia.<br></font></div>
</div>