MATLAB和R及Python亚群差异表达分析

在进行亚群差异表达分析以检测和分析选择性扰动效应时，MATLAB 可以通过一系列步骤来实现这一目的。以下是如何使用 MATLAB 进行这样的分析的简要步骤和示例。

MATLAB片段

1. 数据导入和预处理

确保数据已经格式化为适合 MATLAB 的结构，例如矩阵或表格格式。通常，这涉及将基因表达矩阵（如单细胞 RNA-seq 数据）导入 MATLAB，并将亚群标签和扰动信息相匹配。

% 导入数据，假设数据保存在 CSV 文件中
expressionData = readmatrix('expression_data.csv'); % 行表示基因，列表示样本
groupLabels = readmatrix('group_labels.csv');       % 样本的亚群标签
perturbationLabels = readmatrix('perturbation_labels.csv'); % 样本的扰动标签

2. 数据标准化和归一化

为了使数据可比，通常需要对其进行归一化和标准化处理。

% 数据归一化处理，如 log 变换或 z-score 标准化
expressionData = log2(expressionData + 1); % 避免 log(0)
zScoredData = zscore(expressionData, 0, 2); % 对每个基因进行标准化

3. 差异表达分析

使用 t 检验或其他统计方法来比较不同亚群之间的表达差异。

% 初始化结果存储
pValues = zeros(size(expressionData, 1), 1);
foldChange = zeros(size(expressionData, 1), 1);

% 比较两个条件下的表达水平，例如：扰动组 vs 对照组
for i = 1:size(expressionData, 1)
    group1 = expressionData(i, perturbationLabels == 1); % 扰动组
    group2 = expressionData(i, perturbationLabels == 0); % 对照组
    
    % 使用 t 检验进行差异分析
    [~, pValues(i)] = ttest2(group1, group2);
    
    % 计算 fold change
    foldChange(i) = mean(group1) - mean(group2);
end

% 多重比较校正（Benjamini-Hochberg 法）
[~, ~, adjPValues] = fdr_bh(pValues);

4. 结果可视化

绘制火山图和热图来展示差异表达的结果。

% 火山图
figure;
scatter(foldChange, -log10(adjPValues));
xlabel('Fold Change');
ylabel('-log10(Adjusted P-Value)');
title('Volcano Plot');
grid on;

% 热图绘制显著差异表达的基因
significantGenes = adjPValues < 0.05; % 选择显著性基因
figure;
heatmap(zScoredData(significantGenes, :));
title('Heatmap of Differentially Expressed Genes');

5. 进一步分析

可以对结果进行更多下游分析，例如基因富集分析（使用 MATLAB 的 gsea 工具箱或与 R 工具结合）来解释扰动的生物学含义。

% 假设你有 MATLAB Bioinformatics Toolbox
% 可以使用 clustergram 进行聚类分析
clustergram(zScoredData(significantGenes, :), 'Standardize', 2);

6. 解释和报告结果

对识别出的显著基因进行注释，并将其映射到生物学通路以获得更深入的理解。

这个流程概述了如何使用 MATLAB 进行差异表达分析，以检测和分析亚群中的选择性扰动效应。更多高级分析可以通过集成其他 MATLAB 工具箱或外部工具（如 R 中的 DESeq2 或 Seurat）来实现。

R片段

在 R 中进行亚群差异表达分析，以检测和分析选择性扰动效应时，可以使用各种工具和包，如 Seurat、DESeq2、edgeR 等。下面是一个典型的分析流程示例：

1. 数据导入和预处理

通常情况下，单细胞 RNA 测序数据会使用 Seurat 来处理。首先，导入数据并创建一个 Seurat 对象。

library(Seurat)

# 假设你有一个表达矩阵，行是基因，列是样本
expression_matrix <- read.csv("expression_data.csv", row.names = 1)
meta_data <- read.csv("metadata.csv", row.names = 1)

# 创建 Seurat 对象
seurat_obj <- CreateSeuratObject(counts = expression_matrix, meta.data = meta_data)

# 查看数据
head(seurat_obj@meta.data)

2. 数据标准化和标记

使用 Seurat 进行数据归一化和标记。

# 归一化数据
seurat_obj <- NormalizeData(seurat_obj)

# 找到高变基因
seurat_obj <- FindVariableFeatures(seurat_obj)

# 数据缩放
seurat_obj <- ScaleData(seurat_obj, features = rownames(seurat_obj))

3. 亚群识别和注释

如果你还没有亚群信息，你可以使用 Seurat 的聚类方法来识别亚群。

# 运行 PCA
seurat_obj <- RunPCA(seurat_obj, features = VariableFeatures(object = seurat_obj))

# 聚类并标记亚群
seurat_obj <- FindNeighbors(seurat_obj, dims = 1:10)
seurat_obj <- FindClusters(seurat_obj, resolution = 0.5)

# 运行 UMAP/T-SNE 进行可视化
seurat_obj <- RunUMAP(seurat_obj, dims = 1:10)
DimPlot(seurat_obj, reduction = "umap", group.by = "seurat_clusters")

4. 差异表达分析

进行不同条件（例如扰动和对照组）之间的差异表达分析。

# 识别扰动效应
Idents(seurat_obj) <- "perturbation_condition"  # 假设元数据中有条件标签

# 寻找差异表达基因
diff_expr_results <- FindMarkers(seurat_obj, ident.1 = "treated", ident.2 = "control")

# 查看结果
head(diff_expr_results)

5. 多重比较校正和过滤显著基因

应用多重比较校正来过滤显著的基因。

# 将 P 值进行调整
diff_expr_results$adjusted_p_val <- p.adjust(diff_expr_results$p_val, method = "BH")

# 筛选显著基因
significant_genes <- subset(diff_expr_results, adjusted_p_val < 0.05)

6. 结果可视化

绘制火山图和热图来展示差异表达分析结果。

library(ggplot2)

# 火山图
ggplot(diff_expr_results, aes(x = avg_log2FC, y = -log10(adjusted_p_val))) +
  geom_point(alpha = 0.4) +
  theme_minimal() +
  xlab("Log2 Fold Change") +
  ylab("-log10 Adjusted P-Value") +
  ggtitle("Volcano Plot of Differential Expression Analysis")

# 热图
library(pheatmap)

top_genes <- rownames(significant_genes)[1:20]  # 选择前 20 个显著基因
pheatmap(seurat_obj@assays$RNA@scale.data[top_genes, ], cluster_rows = TRUE, cluster_cols = TRUE)

7. 下游分析

执行基因集富集分析来解释选择性扰动效应的生物学意义。可以使用 clusterProfiler 或其他工具。

library(clusterProfiler)
library(org.Hs.eg.db)

# 转换基因名称
gene_list <- rownames(significant_genes)
gene_list <- bitr(gene_list, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)

# 运行基因集富集分析
gse_result <- enrichGO(gene = gene_list$ENTREZID, OrgDb = org.Hs.eg.db, ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05)

# 查看结果
head(gse_result)

8. 报告和解释

将所有结果整合到报告中，并注释显著基因及其功能路径，以便获得更深入的生物学理解。

这个流程为你提供了使用 R 进行亚群差异表达分析来检测和分析选择性扰动效应的完整步骤。通过 R 中的包（如 Seurat 和 clusterProfiler）可以高效地进行数据处理、可视化和功能分析。

👉更新：亚图跨际