文章目錄
- MATLAB 基因表達數據處理與可視化全流程案例
- 一、數據準備與預處理
- 1. 模擬/加載基因表達數據
- 二、高維基因表達數據降維可視化
- 三、差異表達基因分析與可視化
- 四、基因表達時序數據可視化（擴展）
- 五、關鍵技巧與注意事項
- 六、真實數據加載示例
MATLAB 基因表達數據處理與可視化全流程案例
基因表達數據（如RNA-seq、微陣列數據）通常是高維稀疏矩陣（樣本數×基因數，維度可達上萬），核心需求包括：數據預處理、降維可視化、差異表達分析、聚類分析等。以下是完整可運行的實戰案例，涵蓋從數據加載到可視化的全流程。
matlab數據分析教程包含環境配置（matlab下載安裝配置教程）:
一、數據準備與預處理
1. 模擬/加載基因表達數據
首先生成模擬的基因表達數據（或加載真實數據，如GEO數據庫下載的.txt/.csv文件）：
% ===================== 步驟1：生成模擬基因表達數據 =====================
rng(1); % 固定隨機種子
n_samples = 50; % 50個樣本（如患者/對照組）
n_genes = 1000; % 1000個基因（高維特徵）
% 分組：20個對照組（Control）、30個處理組（Treatment）
group = categorical([repmat('Control',n_samples*0.4,1); repmat('Treatment',n_samples*0.6,1)]);
% 生成表達矩陣：處理組部分基因上調/下調，模擬差異表達
expr_matrix = randn(n_samples, n_genes); % 基礎表達（正態分佈）
% 前50個基因在處理組中上調（倍數2）
expr_matrix(group=='Treatment', 1:50) = expr_matrix(group=='Treatment', 1:50) + 2;
% 51-100個基因在處理組中下調（倍數2）
expr_matrix(group=='Treatment', 51:100) = expr_matrix(group=='Treatment', 51:100) - 2;
% ===================== 步驟2：數據預處理（關鍵！） =====================
% 1. 過濾低表達基因（去除表達量接近0的基因，減少噪聲）
expr_filtered = expr_matrix(:, mean(expr_matrix)>0.1); % 保留均值>0.1的基因
% 2. 標準化（z-score，消除基因間量綱差異）
expr_norm = zscore(expr_filtered);
% 3. 處理缺失值（如存在NaN，用均值填充）
expr_norm(isnan(expr_norm)) = 0; % 簡化處理，實際可⽤knnimpute二、高維基因表達數據降維可視化
基因表達數據維度極高（上千/上萬基因），需通過降維映射到2D/3D空間，展示樣本聚類：
% ===================== 步驟3：降維（UMAP/PCA/t-SNE） =====================
% 方法1：PCA（快速看全局分佈）
[coeff, score, ~, ~, explained] = pca(expr_norm);
X_pca = score(:,1:2); % 前2主成分
% 方法2：UMAP（兼顧全局+局部，推薦基因數據）
X_umap = umap(expr_norm, 'NumDimensions', 2, 'NumNeighbors', 10);
% 方法3：t-SNE（聚焦局部聚類）
X_tsne = tsne(expr_norm, 'Perplexity', 20, 'Algorithm','barnes-hut');
% ===================== 步驟4：可視化對比 =====================
figure('Color','w','Position',[100,100,1200,400]);
tiledlayout(1,3,'Padding','compact');
% 子圖1：PCA結果
nexttile;
gscatter(X_pca(:,1), X_pca(:,2), group, [0.2 0.6 0.8; 0.8 0.2 0.4], 'os', 8);
xlabel(['PC1 (解釋方差: ', num2str(explained(1),2), '%)']);
ylabel(['PC2 (解釋方差: ', num2str(explained(2),2), '%)']);
title('PCA降維：基因表達樣本分佈');
grid on; legend('Location','best');
% 子圖2：UMAP結果
nexttile;
gscatter(X_umap(:,1), X_umap(:,2), group, [0.2 0.6 0.8; 0.8 0.2 0.4], 'os', 8);
xlabel('UMAP維度1'); ylabel('UMAP維度2');
title('UMAP降維：基因表達樣本分佈');
grid on; legend('Location','best');
% 子圖3：t-SNE結果
nexttile;
gscatter(X_tsne(:,1), X_tsne(:,2), group, [0.2 0.6 0.8; 0.8 0.2 0.4], 'os', 8);
xlabel('t-SNE維度1'); ylabel('t-SNE維度2');
title('t-SNE降維：基因表達樣本分佈');
grid on; legend('Location','best');
sgtitle('基因表達高維數據降維可視化','FontSize',14,'FontWeight','bold');結果分析：
- PCA 能看到兩組樣本的大致分離，但因基因表達的非線性特徵，可能存在重疊；
- UMAP 完全分離對照組/處理組，且保留樣本間的全局相對關係（推薦基因數據）；
- t-SNE 局部聚類更緊湊，但全局結構可能失真。
三、差異表達基因分析與可視化
識別兩組間表達量顯著差異的基因，並可視化其表達模式：
% ===================== 步驟5：差異表達基因篩選（t檢驗） =====================
% 分組提取表達數據
expr_control = expr_norm(group=='Control', :);
expr_treatment = expr_norm(group=='Treatment', :);
% 雙樣本t檢驗（計算p值）
[p_val, ~] = ttest2(expr_control, expr_treatment);
% 校正p值（多重檢驗校正，Benjamini-Hochberg）
p_adjust = mafdr(p_val, 'BHFDR', true);
% 篩選差異基因：校正p值<0.05 且 表達量差異>1.5倍
fold_change = mean(expr_treatment) - mean(expr_control); % 差異倍數
diff_genes_idx = (p_adjust < 0.05) & (abs(fold_change) > 1.5);
diff_genes_expr = expr_norm(:, diff_genes_idx); % 差異基因表達矩陣
fprintf('篩選出 %d 個差異表達基因\n', sum(diff_genes_idx));
% ===================== 步驟6：差異基因可視化 =====================
% 1. 火山圖（展示p值 vs 差異倍數）
figure('Color','w');
scatter(fold_change, -log10(p_adjust), 30, 'k', 'filled');
hold on;
% 標記顯著差異基因（紅色：上調，綠色：下調）
scatter(fold_change(diff_genes_idx), -log10(p_adjust(diff_genes_idx)), ...
50, fold_change(diff_genes_idx)>0, 'filled');
colormap([0.8 0.2 0.2; 0.2 0.8 0.2]); % 紅=上調，綠=下調
% 添加閾值線
plot([-1.5, -1.5], [0, max(-log10(p_adjust))], 'k--');
plot([1.5, 1.5], [0, max(-log10(p_adjust))], 'k--');
plot([min(fold_change), max(fold_change)], [-log10(0.05), -log10(0.05)], 'k--');
xlabel('差異倍數（Treatment - Control）');
ylabel('-log10(校正p值)');
title('基因表達火山圖');
grid on; colorbar('Ticks',[0,1],'TickLabels',{'下調','上調'});
% 2. 熱圖（展示差異基因的表達模式聚類）
figure('Color','w');
% 對差異基因表達矩陣歸一化（每行歸一化，突出樣本間模式）
diff_genes_norm = zscore(diff_genes_expr, 0, 1);
% 繪製熱圖
heatmap(diff_genes_norm, 'RowLabels',group, 'Colormap',parula, ...
'Title','差異表達基因熱圖（樣本×基因）');結果分析：
- 火山圖中，紅色點為上調基因、綠色點為下調基因，虛線為篩選閾值；
- 熱圖可清晰看到：對照組/處理組的差異基因表達模式完全分離，驗證篩選有效性。
四、基因表達時序數據可視化（擴展）
若為時序基因表達數據（如不同時間點的表達量），可可視化基因表達趨勢：
% ===================== 模擬時序基因表達數據 =====================
time_points = 0:2:10; % 0/2/4/6/8/10小時
n_time = length(time_points);
n_genes = 5; % 選5個核心差異基因
% 模擬5個基因的時序表達趨勢
expr_time = zeros(n_time, n_genes);
expr_time(:,1) = 1 + sin(time_points/2); % 正弦趨勢
expr_time(:,2) = linspace(0, 3, n_time); % 線性上調
expr_time(:,3) = linspace(3, 0, n_time); % 線性下調
expr_time(:,4) = 2 + randn(size(time_points))*0.1; % 穩定表達
expr_time(:,5) = 0.5*cos(time_points) + 1; % 餘弦趨勢
% ===================== 時序趨勢可視化 =====================
figure('Color','w');
plot(time_points, expr_time, 'o-', 'LineWidth',1.5, 'MarkerSize',8);
xlabel('時間（小時）'); ylabel('標準化表達量');
title('基因表達時序趨勢');
legend({'基因1','基因2','基因3','基因4','基因5'}, 'Location','best');
grid on;五、關鍵技巧與注意事項
- 數據預處理核心：
- 過濾低表達基因：避免噪聲基因干擾降維和差異分析；
- 標準化：基因表達量範圍差異大，z-score是基礎操作；
- 缺失值處理：推薦用
knnimpute（K近鄰填充）替代簡單均值填充。
- 降維方法選型：
- 基因數據優先選 UMAP（兼顧全局/局部，速度快）；
- 快速探索用PCA，聚焦局部聚類用t-SNE；
- 大數據（>1000樣本）：先PCA降維到50維，再用UMAP/t-SNE（減少計算量）。
- 差異分析優化：
- 小樣本用非參數檢驗（如
ranksum）替代t檢驗； - 批量數據用
edgeR/DESeq2（MATLAB可調用R包，需安裝Statistics and Machine Learning Toolbox）。
- 可視化優化：
- 熱圖行/列聚類：
heatmap可添加'ClusterRows',true,'ClusterColumns',true； - 火山圖標註關鍵基因：用
text函數標註高差異倍數的基因名。
六、真實數據加載示例
若從GEO數據庫下載了基因表達數據（.txt格式），加載代碼：
% 加載真實基因表達數據（樣本×基因，第一列是分組標籤）
data = readtable('GSExxxx_expr.txt', 'TextType','string');
group = categorical(data(:,1)); % 分組標籤
expr_matrix = table2array(data(:,2:end)); % 表達矩陣
% 後續預處理/降維/可視化同上述案例如需針對特定基因表達分析場景（如單細胞RNA-seq、加權基因共表達網絡分析WGCNA）的定製化代碼，可補充説明，我會進一步細化！
本文章為轉載內容，我們尊重原作者對文章享有的著作權。如有內容錯誤或侵權問題，歡迎原作者聯繫我們進行內容更正或刪除文章。
