Dna-Ausrichtung in numpy Array mit Biopython transformieren

Question

Ich habe mehrere DNA-Sequenzen, die ausgerichtet wurden und ich möchte nur die Basen, die an einer bestimmten Position variabel sind.

Dies könnte möglicherweise getan werden, wenn wir zuerst die Ausrichtung in ein Array transformieren. Ich habe versucht, den Code im Biopython-Tutorial zu verwenden, aber es gibt einen Fehler.

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) 

Der Fehler erhalte ich:

Traceback (most recent call last): 

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

Answer 1

AlignIO nicht das Werkzeug, das Sie für diesen Job wollen zu sein scheint. Sie haben vermutlich eine Datei mit vielen Sequenzen, nicht mit vielen mehrfachen Sequenz-Alignments, also möchten Sie wahrscheinlich SeqIO, nicht AlignIO (source) verwenden. Aus diesem Grund ist die Form Ihres Arrays (1, 99, 16926), weil Sie 1 Ausrichtung von 99 Sequenzen der Länge 16926 haben.

Wenn Sie nur ein Array der Sequenzen (die Sie aus dem erscheint np.character dtype zu np.array geliefert), dann wie folgt vorgehen:

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) 

Hinweis darüber, dass technisch jedes Element records ist auch ein biopython SeqRecord die die Sequenz zusätzlich zu Metadaten enthält. list(record) ist eine Abkürzung, um die Sequenz zu erhalten, die andere Art ist record.seq. Entweder sollte funktionieren, aber ich entschied mich für das Attribut way, da es expliziter ist.

Answer 2

Ich beantworte Ihr Problem, anstatt Ihren Code zu reparieren. Wenn Sie nur bestimmte Positionen halten wollen, möchten Sie AlignIO verwenden:

FASTA Probe al.fas:

>seq1 
CATCGATCAGCATCGACATGCGGCA-ACG 
>seq2 
CATCGATCAG---CGACATGCGGCATACG 
>seq3 
CATC-ATCAGCATCGACATGCGGCATACG 
>seq4 
CATCGATCAGCATCGACAAACGGCATACG 

Angenommen, Sie nur bestimmte Positionen halten wollen. MultipleSeqAlignment können Sie Abfrage die Ausrichtung wie ein numpy Array:

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]) 

Wenn Sie die Form der Ausrichtung wissen wollen, verwenden Sie len und get_alignment_length:

>>> print(len(al), al.get_alignment_length()) 
4 29 

---

Wenn Sie AlignIO.parse(), um eine Ausrichtung zu laden, nimmt sie an, dass die zu analysierende Datei mehr als eine Ausrichtung enthalten könnte (PHYLIP tut dies). Daher gibt der Parser einen Iterator über jede Ausrichtung und nicht über Datensätze zurück, wie es Ihr Code impliziert. Aber Ihre FASTA-Datei enthält nur eine Ausrichtung pro Datei und parse() ergibt nur eine MultipleSeqAlignment. Also die Lösung für Ihren Code ist:

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