DNA序列状态转移矩阵分析及可视化

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

def read_fasta_file(file_name): sequences = [] with open(file_name, 'r') as f: sequence = '' for line in f: if line.startswith('>'): if sequence != '': sequences.append(sequence) sequence = '' else: sequence += line.strip() if sequence != '': sequences.append(sequence) return sequences

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

def get_transition_counts(sequences): num_states = 4 transition_counts = np.zeros((num_states, num_states)) for sequence in sequences: current_state = get_base_state(sequence[0]) for i in range(1, len(sequence)): next_state = get_base_state(sequence[i]) transition_counts[current_state, next_state] += 1 current_state = next_state return transition_counts

def get_transition_matrix(transition_counts): row_sums = transition_counts.sum(axis=1) transition_matrix = transition_counts / row_sums[:, np.newaxis] return transition_matrix

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

def standardize(matrix): standardized_matrix = (matrix - np.mean(matrix, axis=0)) / np.std(matrix, axis=0) return standardized_matrix

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

def plot_pac(matrix, label=''): 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(label) plt.show()

if name == 'main': fasta_folder = './FASTA 文件' fasta_files = [os.path.join(fasta_folder, f) for f in os.listdir(fasta_folder) if f.endswith('.fasta')]

for fasta_file in fasta_files:
    matrix = calc_transition_matrix(fasta_file)
    print(matrix)

    standardized_matrix = standardize(matrix)
    print(standardized_matrix)

    pca_result = pca(standardized_matrix)

    # 修改后的绘图代码
    plt.plot(pca_result[:, 0], pca_result[:, 1], label=fasta_file)

    plot_pac(standardized_matrix, label=fasta_file)

plt.legend()
plt.show()