DNA 序列状态转移频率矩阵分析 - Python 代码实现

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 = ['a', 'c', 'g', 't']
    states = len(alphabet)
    matrix = np.zeros((states, states))
    for seq in sequences:
        seq_str = str(seq.seq).lower()
        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)

该代码可以将 DNA 序列看成一个有 4 个状态 (a, c, g, t) 的随机过程，计算每个碱基对之间状态转移的频率矩阵，并进行 z-score 标准化和 PCA 分析，最终绘制 PAC 图。

代码说明：

• calc_transition_matrix 函数读取 FASTA 文件，计算每个碱基对之间状态转移频率矩阵。
• standardize 函数对矩阵进行 z-score 标准化处理。
• plot_pac 函数绘制 PAC 图。
• 主程序部分遍历文件夹中的所有 FASTA 文件，依次进行状态转移矩阵计算、标准化和 PAC 图绘制。

使用说明：

1. 将 FASTA 文件放在名为 FASTA 文件 的文件夹中。
2. 运行代码，程序会自动读取所有 FASTA 文件并进行分析。
3. 程序会输出每个 FASTA 文件对应的状态转移频率矩阵、标准化矩阵和 PAC 图。

注意：

• 该代码需要安装 Biopython、numpy 和 matplotlib 库。
• 代码中的 FASTA 文件 文件夹路径需要根据实际情况进行修改。
• PAC 图的解释需要参考相关文献和资料。