Converted phylo from muscle to biotite

Author: msauria

README.md

@@ -24,13 +24,12 @@ VCF-kit is a command-line based collection of utilities for performing analysis
 
 __VCF-Kit has been upgraded to Python 3__
 
-VCF-kit has been tested with Python 3.6. VCF-kit makes use of additional software for a variety of tasks:
+VCF-kit has been tested with Python 3.10. VCF-kit makes use of additional software for a variety of tasks:
 
 * bwa (v 0.7.12)
 * samtools (v 1.3)
 * bcftools (v 1.3)
 * blast (v 2.2.31+)
-* muscle (v 3.8.31)
 * primer3 (v 2.5.0)
 
 You can install these dependencies and VCF-kit using conda, or you can use a Docker image.
@@ -46,7 +45,6 @@ conda create -n vcf-kit \
   "samtools>=1.10" \
   "bcftools>=1.10" \
   "blast>=2.2.31" \
-  "muscle>=3.8.31" \
   "primer3>=2.5.0"
 
 conda activate vcf-kit

docs/phylo.md

@@ -5,7 +5,7 @@ vk phylo fasta <vcf> [<region>]
 vk phylo tree (nj|upgma) [--plot] <vcf> [<region>]
 ```
 
-The `phylo` command can be used to generate dendrograms, tree files, or a fasta file of variants concatenated together (equivelent to a multiple sequence alignment) from a VCF. Tree files are generated in [Newick format](http://evolution.genetics.washington.edu/phylip/newicktree.html)) with [MUSCLE](http://drive5.com/muscle/) using [UPGMA](https://en.wikipedia.org/wiki/UPGMA) or [neighbor-joining](https://en.wikipedia.org/wiki/Neighbor_joining). VCF-kit can use the output tree file to  generate a plot of the tree/phylogeny.
+The `phylo` command can be used to generate dendrograms, tree files, or a fasta file of variants concatenated together (equivelent to a multiple sequence alignment) from a VCF. Tree files are generated in [Newick format](http://evolution.genetics.washington.edu/phylip/newicktree.html)) with [biotite](https://www.biotite-python.org/latest/index.html) using [UPGMA](https://en.wikipedia.org/wiki/UPGMA) or [neighbor-joining](https://en.wikipedia.org/wiki/Neighbor_joining). VCF-kit can use the output tree file to  generate a plot of the tree/phylogeny.
 
 `phylo` can read a VCF directly or from stdin by using `-`.
 
@@ -38,7 +38,7 @@ bcftools filter --set-GTs . --exclude 'FMT/DP < 20'  data/test.vcf.gz | vk phylo
 
 ## Generating a tree/phylogeny
 
-`vk phylo tree` can be used to generate a tree/phylogeny from a vcf file. This command uses a fasta file (identical to what is produced using `vk phylo fasta`), and uses [MUSCLE](https://en.wikipedia.org/wiki/MUSCLE_(alignment_software) to produce a tree file. 
+`vk phylo tree` can be used to generate a tree/phylogeny from a vcf file. This command uses a fasta file (identical to what is produced using `vk phylo fasta`), and uses [biotite](https://www.biotite-python.org/latest/index.html) to produce a tree file. 
 
 ### Generate a UPGMA tree
 

requirements.txt

@@ -8,6 +8,7 @@ cyvcf2
 docopt
 biopython==1.77
 logzero
+biotite
 pomegranate
 clint
 requests

vcfkit/phylo.py

@@ -18,6 +18,11 @@
 
 import numpy as np
 from jinja2 import Template
+import biotite.sequence.align as align, phylo
+import biotite.sequence.phylo as phylo
+from biotite.sequence import NucleotideSequence
+import networkx as nx
+import matplotlib.pyplot as plt
 
 from clint.textui import colored, indent, puts_err
 from docopt import docopt
@@ -52,56 +57,71 @@ def first(s):
         """
             Generate an aligned fasta from a VCF file.
         """
-        gt_set = np.chararray((0,len(v.samples)))
         gt_set = []
         for line in variant_set:
             if line.is_snp:
                 gt_set.append(firstv(line.gt_bases))
         if len(gt_set) == 0:
             exit(puts_err("No genotypes"))
         gt_set = np.vstack(gt_set)
-        seqs = list(zip(v.samples, np.transpose(gt_set)))
+        seqs = [''.join(gt_set[:, i]) for i in range(gt_set.shape[1])]
+
         if args["fasta"]:
-            for sample, seq in seqs:
-                print((">" + sample))
-                print((''.join(seq)))
+            for sample, seq in list(zip(v.samples, seqs)):
+                print(">" + sample)
+                print(seq)
 
         elif args["tree"]:
             """
             Generate a phylogenetic tree using an aligned fasta with muscle.
             """
 
-            ########
-            # Consider biotite to replace deprecated --maketree function in muscle 
-            ########
-
-            # Check for muscle dependency
-            check_program_exists("muscle")
-            fasta = ""
             with indent(4):
                 puts_err(colored.blue("\nGenerating Fasta\n"))
-            for sample, seq in seqs:
-                fasta += ">" + sample + "\n" + ''.join(seq) + "\n"
-            tree_type = "upgma"  # default is upgma
+            trace = align.Alignment.trace_from_strings(seqs)
+            aligned_seqs = align.Alignment([NucleotideSequence(seq) for seq in seqs], trace)
+            distances = 1 - align.get_pairwise_sequence_identity(aligned_seqs, mode="all")
+
             if args["nj"]:
-                tree_type = "neighborjoining"
-            with indent(4):
-                puts_err(colored.blue("\nGenerating " + tree_type + " Tree\n"))
-            comm = ["muscle", "-maketree", "-in", "-", "-cluster", tree_type]
-            tree, err = Popen(comm, stdin=PIPE, stdout=PIPE).communicate(input=fasta.encode())
-            
+                tree = phylo.neighbor_joining(distances)
+            else:
+                tree = phylo.upgma(distances)
+
             # output tree
-            print(tree.decode("utf-8"))
+            print(tree.to_newick(labels=v.samples))
 
             if args["--plot"]:
+                # graph = tree.as_graph().to_undirected()
+                # fig = plt.figure(figsize=(8.0, 8.0))
+                # ax = fig.gca()
+                # ax.axis("off")
+                # # Calculate position of nodes in the plot
+                # pos = nx.kamada_kawai_layout(graph)
+                # # Assign the gene names to the nodes that represent a reference index
+                # node_labels = {i: name for i, name in enumerate(v.samples)}
+                # nx.draw_networkx_edges(graph, pos, ax=ax)
+                # nx.draw_networkx_labels(
+                #     graph,
+                #     pos,
+                #     ax=ax,
+                #     labels=node_labels,
+                #     font_size=7,
+                #     # Draw a white background behind the labeled nodes
+                #     # for better readability
+                #     bbox=dict(pad=0, color="white"),
+                # )
+                # fig.tight_layout()
+
+                # plt.show()
+
                 prefix = os.path.dirname(os.path.abspath(sys.modules['vcfkit'].__file__)) + "/static"
                 template = open(prefix + "/tree.html",'r').read()
                 tree_template = Template(template)
                 html_out = tempfile.NamedTemporaryFile(suffix=".html", delete=False)
                 with html_out as f:
-                    tree = tree.replace("\n", "")
+                    tree = tree.to_newick(labels=v.samples)
                     sample_len = len(v.samples)
-                    f.write(tree_template.render(**locals()))
+                    f.write(tree_template.render(tree=tree, prefix=prefix, sample_len=sample_len).encode())
                     webbrowser.open("file://" + html_out.name)
 
 if __name__ == '__main__':