MATLAB和Python多语高维转录分析
MATLAB高维转录分析代码片段
在 MATLAB 中进行高维转录分析（例如单细胞 RNA 测序 scRNA-seq 分析）可以使用以下步骤，包括数据导入、预处理、降维、聚类和差异表达分析。 MATLAB 的矩阵操作能力和工具箱（如统计和机器学习工具箱）提供了高效的分析手段。
以下是高维转录分析的步骤：
1. 数据导入
使用 readmatrix 导入基因表达数据，通常数据表包含基因名称和细胞的表达值。
% 读取数据
data = readmatrix('gene_expression.csv'); % 数据格式：行 = 基因, 列 = 细胞
genes = data(:,1); % 假设第一列为基因名称
counts_matrix = data(:,2:end); % 提取基因表达矩阵
2. 质量控制和预处理
过滤表达较低的细胞和基因以提高分析质量。
% 计算每个细胞和基因的计数
cell_counts = sum(counts_matrix, 1);
gene_counts = sum(counts_matrix, 2);
% 设置过滤标准
min_gene_threshold = 200; % 每个细胞最少基因数
min_cell_threshold = 10; % 每个基因最少细胞数
% 筛选满足条件的细胞和基因
filtered_cells = cell_counts > min_gene_threshold;
filtered_genes = gene_counts > min_cell_threshold;
% 过滤表达矩阵
counts_matrix_filtered = counts_matrix(filtered_genes, filtered_cells);
genes_filtered = genes(filtered_genes);
3. 归一化和对数转换
使用每个细胞的总表达值进行归一化，然后取对数转换。
% 归一化每个细胞的总计数
total_counts_per_cell = sum(counts_matrix_filtered, 1);
normalized_counts = counts_matrix_filtered ./ total_counts_per_cell;
% 对数转换
log_normalized_counts = log1p(normalized_counts);
4. 降维分析（主成分分析 PCA）
使用 PCA 降维并可视化。
% 中心化数据
data_centered = log_normalized_counts - mean(log_normalized_counts, 2);
% 执行 PCA
[coeff, score, ~, ~, explained] = pca(data_centered');
% 选择前两个主成分绘制散点图
scatter(score(:,1), score(:,2));
title('PCA of Gene Expression');
xlabel('PC1');
ylabel('PC2');
5. 聚类分析（K-means 聚类）
对 PCA 降维结果进行聚类（如 K-means）。
k = 5; % 假设分为5类
idx = kmeans(score(:,1:10), k); % 使用前10个主成分
% 绘制聚类结果
gscatter(score(:,1), score(:,2), idx);
title('K-means Clustering on PCA-reduced Data');
xlabel('PC1');
ylabel('PC2');
6. t-SNE 或 UMAP 降维（可选）
MATLAB 提供 tsne 函数来进一步降维，t-SNE 更适合高维数据的可视化。
% 执行 t-SNE
Y = tsne(score(:,1:10)); % 使用前10个主成分
% 绘制 t-SNE 结果
gscatter(Y(:,1), Y(:,2), idx);
title('t-SNE of Gene Expression');
xlabel('t-SNE 1');
ylabel('t-SNE 2');
7. 差异基因表达分析
找到每个簇的标志基因，计算每个基因在不同簇中的平均表达。
% 计算每个簇的平均表达
mean_expression = zeros(size(log_normalized_counts, 1), k);
for i = 1:k
mean_expression(:,i) = mean(log_normalized_counts(:, idx == i), 2);
end
% 找到方差最高的基因作为潜在标志基因
[~, top_genes] = sort(var(mean_expression, 0, 2), 'descend');
top_genes = top_genes(1:10); % 选择前10个基因
% 绘制标志基因表达水平
bar(mean_expression(top_genes, :));
set(gca, 'XTickLabel', genes_filtered(top_genes));
xlabel('Genes');
ylabel('Expression Level');
title('Top Differentially Expressed Genes');
legend('Cluster 1', 'Cluster 2', 'Cluster 3', 'Cluster 4', 'Cluster 5');
总结
此 MATLAB 脚本为高维转录分析提供了基础流程，涵盖了数据导入、质量控制、归一化、降维、聚类和差异基因表达分析。这些步骤可为分析单细胞 RNA 测序数据提供支持，并可根据需求进一步扩展。
Python高维转录分析代码片段
在 Python 中进行高维转录分析（如单细胞 RNA 测序数据分析）通常使用 scanpy 和 seaborn 等库。 scanpy 是一个专门用于高维生物学数据的处理与分析库，它提供了从数据加载、预处理到降维和聚类的一整套分析流程。
以下是使用 Python 进行高维转录分析的主要步骤：
1. 导入必要的库
import scanpy as sc
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
2. 加载数据
假设数据是一个 .csv 文件，其中行表示基因，列表示单细胞的基因表达水平。可以使用 scanpy 将其导入为 AnnData 对象。
# 加载数据
adata = sc.read_csv('gene_expression.csv')
adata.var_names_make_unique() # 确保基因名称唯一
3. 数据预处理
包括过滤低质量的细胞和基因，进行归一化和对数转换等步骤。
# 基因和细胞过滤
sc.pp.filter_cells(adata, min_genes=200) # 保留至少200个基因的细胞
sc.pp.filter_genes(adata, min_cells=3) # 保留至少在3个细胞中表达的基因
# 归一化与对数转换
sc.pp.normalize_total(adata, target_sum=1e4) # 使每个细胞的总表达水平为1e4
sc.pp.log1p(adata) # 对数转换
4. 可视化数据的质量控制
可以绘制一些常用的 QC 图，例如每个细胞的基因数和表达水平的分布。
# 计算一些常用的 QC 指标
adata.obs['n_genes'] = (adata.X > 0).sum(axis=1) # 每个细胞的基因数
adata.obs['total_counts'] = adata.X.sum(axis=1) # 每个细胞的总计数
# 绘制 QC 图
sns.histplot(adata.obs['n_genes'], kde=False, bins=50)
plt.xlabel('Number of genes')
plt.ylabel('Number of cells')
plt.title('Gene count distribution per cell')
plt.show()
sns.histplot(adata.obs['total_counts'], kde=False, bins=50)
plt.xlabel('Total counts')
plt.ylabel('Number of cells')
plt.title('Total count distribution per cell')
plt.show()
5. 高度变量基因的识别
识别高变基因（highly variable genes），这些基因是后续降维分析的重点。
# 选择高变基因
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
adata = adata[:, adata.var['highly_variable']]
6. PCA 降维
scanpy 支持 PCA 分析，并可以选择绘制散点图。
# PCA 降维
sc.tl.pca(adata, svd_solver='arpack')
sc.pl.pca(adata, color='total_counts') # 使用细胞总计数上色
7. t-SNE 和 UMAP 降维
使用 t-SNE 和 UMAP 进行可视化，将数据映射到二维平面。
# 计算t-SNE
sc.tl.tsne(adata)
sc.pl.tsne(adata, color='total_counts')
# 计算UMAP
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40) # 构建邻近图
sc.tl.umap(adata)
sc.pl.umap(adata, color='total_counts')
8. 聚类分析
使用 leiden 算法进行聚类，以区分细胞群体。
# 聚类分析
sc.tl.leiden(adata, resolution=0.5)
sc.pl.umap(adata, color='leiden') # 使用聚类结果上色
9. 差异基因表达分析
可以使用聚类标签进行差异表达分析，找出每个簇的特征基因。
# 差异表达分析
sc.tl.rank_genes_groups(adata, 'leiden', method='t-test')
sc.pl.rank_genes_groups(adata, n_genes=10, sharey=False)
10. 保存和导出分析结果
分析完成后，可以将 AnnData 对象或分析结果保存为文件。
# 保存数据
adata.write('processed_gene_expression.h5ad')
总结
此流程利用了 scanpy 的多种方法来进行高维转录分析，包括预处理、降维、聚类和差异表达分析等步骤。这些步骤可以帮助识别细胞类型、分群以及特征基因。