|
1
|
|
|
from Bio import pairwise2 |
|
2
|
|
|
from Bio.Seq import Seq |
|
3
|
|
|
|
|
4
|
|
|
# Define two sequences |
|
5
|
|
|
from Bio import SeqIO |
|
6
|
|
|
from Bio import AlignIO |
|
7
|
|
|
# Open and read the sequence |
|
8
|
|
|
|
|
9
|
|
|
|
|
10
|
|
|
# Access the sequence and other information |
|
11
|
|
|
##print(f"Sequence ID: {seq1.id}") |
|
12
|
|
|
#print(f"Sequence: {seq1.seq}") |
|
13
|
|
|
#print(f"Description: {seq1.description}") |
|
14
|
|
|
|
|
15
|
|
|
|
|
16
|
|
|
file_path = "output.sto" |
|
17
|
|
|
|
|
18
|
|
|
# Read the alignment |
|
19
|
|
|
alignment = AlignIO.read(file_path, "stockholm") |
|
20
|
|
|
|
|
21
|
|
|
# Loop through the alignment and extract sequences |
|
22
|
|
|
for seq1 in alignment: |
|
23
|
|
|
for seq2 in alignment: |
|
24
|
|
|
print(f"Sequence ID: {seq1.id}") |
|
25
|
|
|
#print(f"Sequence: {record.seq}") |
|
26
|
|
|
#print(f"Description: {record.description}") |
|
27
|
|
|
#print() # Blank line between sequences |
|
28
|
|
|
|
|
29
|
|
|
# Align the sequences using global alignment |
|
30
|
|
|
#seq2 = Seq(seq1.seq.replace('-', '')) |
|
31
|
|
|
|
|
32
|
|
|
#print(seq1.seq) |
|
33
|
|
|
#print(seq2) |
|
34
|
|
|
|
|
35
|
|
|
alignments = pairwise2.align.globalxx(seq1.seq, seq2.seq) |
|
36
|
|
|
best_alignment = alignments[0] |
|
37
|
|
|
from icecream import ic;import sys;ic.configureOutput(outputFunction=lambda *a: print(*a, file=sys.stderr), includeContext=True,contextAbsPath=True);ic.configureOutput(prefix='') |
|
38
|
|
|
ic(best_alignment) |
|
39
|
|
|
# Extract the aligned sequences |
|
40
|
|
|
aligned_seq1, aligned_seq2 = best_alignment[0], best_alignment[1] |
|
41
|
|
|
|
|
42
|
|
|
# Calculate similarity (number of matches / length of alignment) |
|
43
|
|
|
matches = sum(1 for a, b in zip(aligned_seq1, aligned_seq2) if a == b) |
|
44
|
|
|
alignment_length = len(aligned_seq1) |
|
45
|
|
|
similarity = (matches / alignment_length) * 100 |
|
46
|
|
|
|
|
47
|
|
|
|
|
48
|
|
|
# Print the alignment and similarity score |
|
49
|
|
|
#print(f"Alignment:\n{aligned_seq1}\n{aligned_seq2}") |
|
50
|
|
|
print(f"{seq1.id.strip()} Similarity: {similarity:.2f}%") |
|
51
|
|
|
|
|
52
|
|
|
#exit() |
|
53
|
|
|
# Take the best alignment |
|
54
|
|
|
|
mmagnus /
rna-tools
