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><code>AlignIO</code> doesn't seem to be the tool you want for this job. You have a file presumably with many sequences, not with many multiple sequence alignments, so you probably want to use <code>SeqIO</code>, not <code>AlignIO</code> (<a href="http://biopython.org/wiki/AlignIO" rel="nofollow">source</a>).  This is why the shape of your array is (1, 99, 16926), because you have 1 alignment of 99 sequences of length 16926.</p>

<p>If you just want an array of the sequences (which it appears you do from the <code>np.character</code> dtype supplied to <code>np.array</code>), then do the following:</p>

<pre><code>import numpy as np
from Bio import SeqIO
records = SeqIO.parse("ma-all-mito.fa", "fasta")
align_array = np.array([record.seq for record in records], np.character)
print("Array shape %i by %i" % align_array.shape)
# expect to be (99, 16926)
</code></pre>

<p>Note above that technically each element of <code>records</code> is also a BioPython <code>SeqRecord</code> which includes the sequence in addition to metadata.  <code>list(record)</code> is a shortcut for getting the sequence, the other way being <code>record.seq</code>.  Either should work, but I chose using the attribute way since it is more explicit.</p>

python numpy biopython