基于状态转移矩阵和PCA的序列分析

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

# 读取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


# 测试代码
if __name__ == '__main__':
    # 读取fasta文件，计算状态转移矩阵
    matrix1 = calc_transition_matrix('WBKG00000000.fasta')
    matrix2 = calc_transition_matrix('WHPN00000000.fasta')
    matrix3 = calc_transition_matrix('WSRO00000000.fasta')
    # 对矩阵进行标准化
    standardized_matrix1 = standardize(matrix1)
    standardized_matrix2 = standardize(matrix2)
    standardized_matrix3 = standardize(matrix3)
    # 对标准化后的矩阵进行PCA分析
    pca_result1 = pca(standardized_matrix1)
    pca_result2 = pca(standardized_matrix2)
    pca_result3 = pca(standardized_matrix3)
    # 绘制PCA散点图
    plt.scatter(pca_result1[:, 0], pca_result1[:, 1], label='WBKG00000000')
    plt.scatter(pca_result2[:, 0], pca_result2[:, 1], label='WHPN00000000')
    plt.scatter(pca_result3[:, 0], pca_result3[:, 1], label='WSRO00000000')
    plt.legend()
    plt.show()

该代码使用Biopython库读取FASTA文件，并计算状态转移矩阵。然后，使用NumPy库对矩阵进行标准化，并使用PCA方法对标准化后的矩阵进行降维。最后，使用matplotlib库绘制PCA散点图，以可视化不同序列之间的差异。

步骤：

1. **导入库：**导入必要的库，包括Bio.SeqIO用于读取FASTA文件，numpy用于矩阵操作，matplotlib.pyplot用于绘制图像。
2. 计算状态转移矩阵：calc_transition_matrix函数读取FASTA文件，计算所有序列的字母组合，构建状态转移矩阵。
3. 标准化矩阵：standardize函数使用z-score标准化方法对状态转移矩阵进行标准化。
4. PCA分析：pca函数对标准化后的矩阵进行PCA分析，并返回投影结果。
5. **绘制PCA散点图：**测试代码读取三个FASTA文件，计算并标准化状态转移矩阵，然后进行PCA分析，最后使用plt.scatter函数绘制PCA散点图。

该代码可以用于分析不同序列之间的差异，并可视化它们的差异模式。通过修改FASTA文件名和参数设置，可以对其他序列进行分析。