Transform dna alignment into numpy array using biopython
<p>I have several DNA sequences that have been aligned and I would like to keep  only the bases that are variable at a specific position. </p>

<p>This maybe could be done if we first transform the alignment into an array. I tried using the code in the Biopython tutorial but it gives an error.</p>

<pre><code>import numpy as np
from Bio import AlignIO
alignment = AlignIO.parse("ma-all-mito.fa", "fasta")
align_array = np.array([list(rec) for rec in alignment], np.character)
print("Array shape %i by %i" % align_array.shape)
</code></pre>

<p>The error I get:</p>

<pre><code>Traceback (most recent call last):

File "C:/select-snps.py", line 8, in &lt;module&gt;
    print("Array shape %i by %i" % align_array.shape)
TypeError: not all arguments converted during string formatting
</code></pre>

<p>I'm answering to your problem instead of fixing your code. If you want to keep only certain positions, you want to use <code>AlignIO</code>:</p>

<p>FASTA sample <code>al.fas</code>:</p>

<pre><code>&gt;seq1
CATCGATCAGCATCGACATGCGGCA-ACG
&gt;seq2
CATCGATCAG---CGACATGCGGCATACG
&gt;seq3
CATC-ATCAGCATCGACATGCGGCATACG
&gt;seq4
CATCGATCAGCATCGACAAACGGCATACG
</code></pre>

<p>Now suppose you want to keep only certain positions. <a href="http://biopython.org/DIST/docs/api/Bio.Align.MultipleSeqAlignment-class.html" rel="nofollow" title="MultipleSeqAlignment">MultipleSeqAlignment</a> allows you to <em>query</em> the alignment like a numpy array:</p>

<pre><code>from Bio import AlignIO


al = AlignIO.read("al.fas", "fasta")

# Print the 11th column
print(al[:, 10])

# Print the 12-15 columns
print(al[:, 11:14])
</code></pre>

<p>If you want to know the shape of the alignment, use <code>len</code> and <code>get_alignment_length</code>:</p>

<pre><code>&gt;&gt;&gt; print(len(al), al.get_alignment_length())
4 29
</code></pre>

<hr>

<p>When you use <code>AlignIO.parse()</code> to load an alignment, it assumes the file to be parsed could contain more than one alignment (PHYLIP does this). Thus the parser returns an iterator over each alignment and not over records as your code implies. But your FASTA file only contain one alignment per file and <code>parse()</code> yields only one <code>MultipleSeqAlignment</code>. So the fix to your code is:</p>

<pre><code>alignment = AlignIO.read("ma-all-mito.fa", "fasta")
align_array = np.array(alignment, np.character)
print("Array shape %i by %i" % align_array.shape)
</code></pre>

python numpy biopython