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

定义碱基序列

alphabet = ['AA', 'AC', 'AG', 'AT', 'CA', 'CC', 'CG', 'CT', 'GA', 'GC', 'GG', 'GT', 'TA', 'TC', 'TG', 'TT']

读取fasta文件，计算相邻状态转移计数矩阵和状态转移矩阵

def calc_transition_matrix(fasta_file): sequences = list(SeqIO.parse(fasta_file, 'fasta')) matrix = np.zeros((len(alphabet), len(alphabet))) for seq in sequences: seq_str = str(seq.seq) for i in range(len(seq_str) - 1): from_state = alphabet.index(seq_str[i:i+2]) to_state = alphabet.index(seq_str[i+1:i+3]) matrix[from_state, to_state] += 1 # 计算状态转移矩阵 transition_matrix = matrix / np.sum(matrix, axis=1, keepdims=True) return matrix, transition_matrix

对矩阵进行z-score标准化

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

对标准化后的矩阵进行PCA分析

def pca(matrix): cov = np.cov(matrix.T) eig_vals, eig_vecs = np.linalg.eig(cov) idx = np.argsort(eig_vals)[::-1] eig_vecs = eig_vecs[:, idx] projection = np.dot(matrix, eig_vecs) return projection

测试代码

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

# 对每个fasta文件进行相邻状态转移计数矩阵和状态转移矩阵的计算、标准化和PCA分析
for fasta_file in fasta_files:
    # 读取fasta文件，计算相邻状态转移计数矩阵和状态转移矩阵
    matrix, transition_matrix = calc_transition_matrix(fasta_file)
    print('相邻状态转移计数矩阵：')
    print(matrix)
    print('状态转移矩阵：')
    print(transition_matrix)

    # 对矩阵进行标准化
    standardized_matrix = standardize(matrix)
    print('标准化后的相邻状态转移计数矩阵：')
    print(standardized_matrix)

    # 对标准化后的矩阵进行PCA分析
    pca_result = pca(standardized_matrix)

    # 绘制PCA散点图
    plt.scatter(pca_result[:, 0], pca_result[:, 1], label=fasta_file)

plt.legend()
plt.show()