DNA 序列状态转移矩阵 PCA 分析：使用 Python 从 FASTA 文件中提取特征

这是一个使用 Python 脚本从 FASTA 文件中提取 DNA 序列特征并进行 PCA 分析的示例。该脚本使用 Biopython 库解析 FASTA 文件，并使用 NumPy 和 Matplotlib 计算状态转移矩阵，进行 Z-score 标准化和 PCA 分析，最终生成 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


# 测试代码
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 = calc_transition_matrix(fasta_file)
        print(matrix)

        # 对矩阵进行标准化
        standardized_matrix = standardize(matrix)
        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()

该脚本首先定义了三个函数：

1. calc_transition_matrix：该函数读取 FASTA 文件，计算状态转移矩阵。
2. standardize：该函数对状态转移矩阵进行 Z-score 标准化。
3. pca：该函数对标准化后的矩阵进行 PCA 分析。

然后，主程序部分遍历文件夹中的所有 FASTA 文件，对每个文件进行状态转移矩阵的计算、标准化和 PCA 分析，并最后绘制 PCA 散点图。

警告信息

如果你在运行该脚本时遇到 RuntimeWarning: Glyph 25991 missing from current font. ... RuntimeWarning: Glyph 20214 missing from current font. 的警告信息，这是因为当前的字体不支持这些字符。这通常意味着当前的字体不支持这些字符，将使用备用字体来显示它们。如果你不想看到这些警告，你可以尝试更改 matplotlib 的字体设置，或者使用支持这些字符的字体。

使用方法

1. 将该脚本保存为 Python 文件，例如 dna_analysis.py。
2. 将你的 FASTA 文件放入一个名为 fasta_files 的文件夹中。
3. 在终端中运行该脚本：python dna_analysis.py。

该脚本将生成一个 PCA 散点图，显示不同 DNA 序列之间的差异。

注意：

• 该脚本只考虑了 DNA 的四个碱基：A、C、G、T。
• 该脚本仅供参考，你可能需要根据你的具体需求进行修改。