DNA序列分析：利用PCA和PAC分析序列的碱基转移模式

import os
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

def read_fasta_file(file_path):
    sequences = {}
    with open(file_path) as f:
        sequence_name = ''
        sequence = ''
        for line in f:
            if line.startswith('>'):
                if sequence_name:
                    sequences[sequence_name] = sequence
                sequence_name = line.strip()[1:]
                sequence = ''
            else:
                sequence += line.strip()
        if sequence_name:
            sequences[sequence_name] = sequence
    return sequences

def get_transition_matrix(sequences):
    # 初始化状态转移频率矩阵
    transition_matrix = np.zeros((4, 4))
    # 遍历所有序列
    for sequence in sequences.values():
        # 遍历序列中的所有碱基对
        for i in range(len(sequence)-1):
            # 获取当前碱基和下一个碱基的状态
            current_state = get_base_state(sequence[i])
            next_state = get_base_state(sequence[i+1])
            # 更新状态转移频率矩阵
            transition_matrix[current_state, next_state] += 1
    # 将状态转移频率矩阵转换为状态转移概率矩阵
    transition_matrix = transition_matrix / np.sum(transition_matrix, axis=1, keepdims=True)
    return transition_matrix

def get_base_state(base):
    if base == 'A':
        return 0
    elif base == 'C':
        return 1
    elif base == 'G':
        return 2
    elif base == 'T':
        return 3
    else:
        return -1

def count_adjacent_transitions(sequences):
    # 初始化相邻状态转移计数器
    transition_count = {
        'AA': 0,
        'AC': 0,
        'AG': 0,
        'AT': 0,
        'CA': 0,
        'CC': 0,
        'CG': 0,
        'CT': 0,
        'GA': 0,
        'GC': 0,
        'GG': 0,
        'GT': 0,
        'TA': 0,
        'TC': 0,
        'TG': 0,
        'TT': 0,
    }
    # 遍历所有序列
    for sequence in sequences.values():
        # 遍历序列中的所有碱基对
        for i in range(len(sequence)-1):
            # 获取当前碱基和下一个碱基的状态
            current_state = get_base_state(sequence[i])
            next_state = get_base_state(sequence[i+1])
            # 判断当前和下一个碱基是否都合法
            if current_state != -1 and next_state != -1:
                # 统计相邻状态转移
                transition_key = '{}{}'.format(get_state_base(current_state), get_state_base(next_state))
                transition_count[transition_key] += 1
    return transition_count

def get_state_base(state):
    if state == 0:
        return 'A'
    elif state == 1:
        return 'C'
    elif state == 2:
        return 'G'
    elif state == 3:
        return 'T'
    else:
        raise ValueError('Invalid state: {}'.format(state))

def plot_pac(transition_matrix, file_name):
    # 计算PAC
    pac = np.zeros((4, 4))
    for i in range(4):
        for j in range(4):
            pac[i, j] = transition_matrix[i, j] * transition_matrix[j, i] / (transition_matrix[i, i] * transition_matrix[j, j])
    # 绘制PAC图
    fig, ax = plt.subplots()
    im = ax.imshow(pac, cmap='coolwarm', vmin=0, vmax=1)
    ax.set_xticks(np.arange(4))
    ax.set_yticks(np.arange(4))
    ax.set_xticklabels(['A', 'C', 'G', 'T'])
    ax.set_yticklabels(['A', 'C', 'G', 'T'])
    ax.set_title('PAC for {}'.format(file_name))
    cbar = ax.figure.colorbar(im, ax=ax)
    cbar.ax.set_ylabel('PAC value', rotation=-90, va="bottom")
    plt.show()

# 遍历文件夹中的所有fasta文件
folder_path = 'FASTA 文件'
for file_name in os.listdir(folder_path):
    if file_name.endswith('.fasta'):
        file_path = os.path.join(folder_path, file_name)
        print('Processing file: {}'.format(file_path))
        # 读取FASTA文件
        sequences = read_fasta_file(file_path)
        # 计算状态转移频率矩阵
        transition_matrix = get_transition_matrix(sequences)
        print('Transition Matrix:')
        print(transition_matrix)
        # 统计两个相邻状态转移的次数
        transition_count = count_adjacent_transitions(sequences)
        print('Adjacent Transition Count:')
        print(transition_count)
        # 标准化状态转移频率矩阵
        standardized_matrix = (transition_matrix - np.mean(transition_matrix)) / np.std(transition_matrix)
        # 进行PCA分析
        pca = PCA(n_components=2)
        pca.fit(standardized_matrix)
        transformed_matrix = pca.transform(standardized_matrix)
        # 绘制散点图
        fig, ax = plt.subplots()
        ax.scatter(transformed_matrix[:, 0], transformed_matrix[:, 1])
        for i, txt in enumerate(['A', 'C', 'G', 'T']):
            ax.annotate(txt, (transformed_matrix[i, 0], transformed_matrix[i, 1]))
        ax.set_xlabel('PC1 ({:.2f}%)'.format(pca.explained_variance_ratio_[0] * 100))
        ax.set_ylabel('PC2 ({:.2f}%)'.format(pca.explained_variance_ratio_[1] * 100))
        ax.set_title('PCA for {}'.format(file_name))
        plt.show()
        # 绘制PAC图
        plot_pac(transition_matrix, file_name)