from Bio import SeqIO import numpy as np import matplotlib.pyplot as plt import os

读取fasta文件，计算状态转移矩阵

def calc_transition_matrix(fasta_file): sequences = list(SeqIO.parse(fasta_file, 'fasta')) alphabet = list(set(''.join([str(seq.seq) for seq in sequences])))
states = len(alphabet) matrix = np.zeros((states, states)) for seq in sequences: seq_str = str(seq.seq) for i in range(len(seq_str) - 1): from_state = alphabet.index(seq_str[i]) to_state = alphabet.index(seq_str[i + 1]) matrix[from_state, to_state] += 1 return matrix

对矩阵进行z-score标准化

def standardize(matrix): standardized_matrix = (matrix - np.mean(matrix, axis=0)) / np.std(matrix, axis=0) return standardized_matrix

绘制PAC图

def plot_pac(matrix, fig_title): cov = np.cov(matrix.T) pac = np.zeros_like(cov) for i in range(cov.shape[0]): for j in range(cov.shape[1]): pac[i, j] = cov[i, j] / np.sqrt(cov[i, i] * cov[j, j]) plt.imshow(pac, cmap='coolwarm') plt.colorbar() plt.title(fig_title) plt.show()

测试代码

if name == 'main': # 获取文件夹中的所有fasta文件 fasta_folder = './FASTA 文件' fasta_files = [os.path.join(fasta_folder, f) for f in os.listdir(fasta_folder) if f.endswith('.fasta')]

# 对每个fasta文件进行状态转移矩阵的计算、标准化、PCA分析和PAC图绘制
for fasta_file in fasta_files:
    # 读取fasta文件，计算状态转移矩阵
    matrix = calc_transition_matrix(fasta_file)
    print(matrix)

    # 对矩阵进行标准化
    standardized_matrix = standardize(matrix)
    print(standardized_matrix)

    # 绘制PAC图
    fig_title = os.path.splitext(os.path.basename(fasta_file))[0]
    plot_pac(standardized_matrix, fig_title)