使用Python进行DNA序列状态转移矩阵分析和PCA降维可视化

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

# 对标准化后的矩阵进行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

# 绘制PAC图
def plot_pac(matrix):
    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.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)

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

        # 绘制PCA散点图
        label = os.path.splitext(os.path.basename(fasta_file))[0]  # 提取文件名作为label
        plt.scatter(pca_result[:, 0], pca_result[:, 1], label=label)

        # 绘制PAC图
        plot_pac(standardized_matrix)

    plt.legend()
    plt.show()

该代码可以从指定文件夹中读取多个FASTA文件，分别计算它们的DNA序列状态转移矩阵，进行标准化，然后进行PCA降维，最后绘制PCA散点图和PAC图。

该代码在绘制PCA散点图时，通过提取文件名作为label，解决了No handles with labels found to put in legend.错误。

你可以根据自己的需求修改代码，例如修改文件路径、修改绘制图形的样式等。