Msa类处理多序列比对数据

同源搜索，多序列比对等都是常用的方式，但是有很多的软件可以实现这些同源搜索和多序列比对，但是不同的软件输出的文件格式却是不完全一致，有熟悉的FASTA格式的，也有A2M, A3M,stockholm等格式。

详细介绍：

https://github.com/soedinglab/hh-suite/wiki#multiple-sequence-alignment-formats

A3M格式文件示例：

1. 每个序列都以 > 开头的行开始，并包含序列的标识信息。
2. 在序列标识行之后，是与该序列相关的比对信息，通常使用字母来表示氨基酸或核酸。‘-’表示缺失，小写字母表示插入。

Stockholm格式文件示例：

import dataclasses
from typing import Sequence, Tuple
import string
import collections# Sequence 表示序列类型，内部的 Sequence[int] 表示整数序列。
# DeletionMatrix 表示一个由整数组成的二维数组。
DeletionMatrix = Sequence[Sequence[int]]### 1. 定义Msa类
# Python中，dataclass 是一个装饰器（Decorator），用于创建称为数据类（data class）的类。
# dataclass 装饰器自动生成一些特殊方法，如 __init__、__repr__、__eq__ 等，
# 减少了编写这些方法的样板代码。
@dataclasses.dataclass(frozen=True)
class Msa:"""Class representing a parsed MSA file."""## 初始化参数sequences: Sequence[str]deletion_matrix: DeletionMatrixdescriptions: Sequence[str]# __post_init__ 是Python数据类（data class）中的特殊方法，# 用于在创建数据类的实例之后进行进一步的初始化操作def __post_init__(self):if not (len(self.sequences) ==len(self.deletion_matrix) ==len(self.descriptions)):raise ValueError('All fields for an MSA must have the same length. 'f'Got {len(self.sequences)} sequences, 'f'{len(self.deletion_matrix)} rows in the deletion matrix and 'f'{len(self.descriptions)} descriptions.')def __len__(self):return len(self.sequences)def truncate(self, max_seqs: int):return Msa(sequences=self.sequences[:max_seqs],deletion_matrix=self.deletion_matrix[:max_seqs],descriptions=self.descriptions[:max_seqs])m_seq = ["AAALLL","AT-LAL","S-ALLI"] # 多序列比对后的数据m_del_matrix = [[0,0,0,0,0,0],[0,0,0,0,0,0],[0,0,0,0,0,0]]m_descriptions = ["seq1","seq2","seq3"]# 实例化
test_msa = Msa(m_seq, m_del_matrix, m_descriptions)
print(test_msa)
print(len(test_msa))
# 去除msa第三条序列
print(test_msa.truncate(2))### 2. 定义函数，解析fasta格式字符串
def parse_fasta(fasta_string: str) -> Tuple[Sequence[str], Sequence[str]]:"""Parses FASTA string and returns list of strings with amino-acid sequences.Arguments:fasta_string: The string contents of a FASTA file.Returns:A tuple of two lists:* A list of sequences.* A list of sequence descriptions taken from the comment lines. In thesame order as the sequences."""sequences = []descriptions = []index = -1for line in fasta_string.splitlines():line = line.strip()if line.startswith('>'):index += 1descriptions.append(line[1:])  # Remove the '>' at the beginning.sequences.append('')continueelif not line:continue  # Skip blank lines.  sequences[index] += linereturn sequences, descriptionswith open("test_aln.a3m") as f:a3m_string = f.read()
sequences, description = parse_fasta(a3m_string)print(sequences)
print(description)## 多序列比对a3m格式：
## 1. 每个序列都以 > 开头的行开始，并包含序列的标识信息。
## 2.在序列标识行之后，是与该序列相关的比对信息，通常使用字母来表示氨基酸或核酸。
##   ‘-’表示缺失，小写字母表示插入。### 3.定义函数，解析a3m格式的msa字符串，生成Msa实例，该函数调用parse_fasta函数
def parse_a3m(a3m_string: str) -> Msa:"""Parses sequences and deletion matrix from a3m format alignment.Args:a3m_string: The string contents of a a3m file. The first sequence in thefile should be the query sequence.Returns:A tuple of:* A list of sequences that have been aligned to the query. Thesemight contain duplicates.* The deletion matrix for the alignment as a list of lists. The elementat `deletion_matrix[i][j]` is the number of residues deleted fromthe aligned sequence i at residue position j.* A list of descriptions, one per sequence, from the a3m file."""sequences, descriptions = parse_fasta(a3m_string)deletion_matrix = []for msa_sequence in sequences:deletion_vec = []deletion_count = 0for j in msa_sequence:if j.islower():deletion_count += 1else:deletion_vec.append(deletion_count)deletion_count = 0deletion_matrix.append(deletion_vec)# Make the MSA matrix out of aligned (deletion-free) sequences.# string.ascii_lowercase, string模块提供的字符串常量，包含了所有小写字母的 ASCII 字符# str.maketrans 是 Python 字符串方法，用于创建一个字符映射表（translation table），# ''换成''并删除string.ascii_lowercasedeletion_table = str.maketrans('', '', string.ascii_lowercase)# str.translate 使用映射表执行字符转换（删除小写字母）aligned_sequences = [s.translate(deletion_table) for s in sequences]return Msa(sequences=aligned_sequences,deletion_matrix=deletion_matrix,descriptions=descriptions)with open("test_aln.a3m") as f:a3m_string = f.read()msa1 = parse_a3m(a3m_string)
print(msa1)### 4.定义函数， 解析stockholm格式的msa字符串，生成Msa实例
def parse_stockholm(stockholm_string: str) -> Msa:"""Parses sequences and deletion matrix from stockholm format alignment.Args:stockholm_string: The string contents of a stockholm file. The firstsequence in the file should be the query sequence.Returns:A tuple of:* A list of sequences that have been aligned to the query. Thesemight contain duplicates.* The deletion matrix for the alignment as a list of lists. The elementat `deletion_matrix[i][j]` is the number of residues deleted fromthe aligned sequence i at residue position j.* The names of the targets matched, including the jackhmmer subsequencesuffix."""## 有序字典，保持多序列比对中的序列顺序name_to_sequence = collections.OrderedDict()for line in stockholm_string.splitlines():line = line.strip()# 去除空行和注释行if not line or line.startswith(('#', '//')):continuename, sequence = line.split()if name not in name_to_sequence:name_to_sequence[name] = ''name_to_sequence[name] += sequencemsa = []deletion_matrix = []query = ''keep_columns = []for seq_index, sequence in enumerate(name_to_sequence.values()):## 第一行为query序列if seq_index == 0:# Gather the columns with gaps from the queryquery = sequencekeep_columns = [i for i, res in enumerate(query) if res != '-']# Remove the columns with gaps in the query from all sequences.aligned_sequence = ''.join([sequence[c] for c in keep_columns])msa.append(aligned_sequence)# Count the number of deletions w.r.t. query.deletion_vec = []deletion_count = 0# query序列相对于每一个同源序列，氨基酸位置的缺失情况，累加连续缺失for seq_res, query_res in zip(sequence, query):                    if seq_res != '-' or query_res != '-':                if query_res == '-':deletion_count += 1else:deletion_vec.append(deletion_count)deletion_count = 0deletion_matrix.append(deletion_vec)return Msa(sequences=msa,deletion_matrix=deletion_matrix,descriptions=list(name_to_sequence.keys()))with open("test_aln.stockholm") as f:stockholm_string = f.read()
print(stockholm_string)msa2 = parse_stockholm(stockholm_string) 
print(msa2)## 注：parse_stockholm 和 parse_a3m 函数生成Msa对象中，
##    deletion_matrix中在查询序列deletion位置填上缺失的个数，
##    下一个氨基酸位置的0跳过，所以总长度相等
##    如函数输入msa中第一条序列（query序列）为：“A--CE-H”， 则函数输出的第一条序列为：“ACEH”， 
##    deletion_matrix的第一个元素为：[0,2,0,1]