MATLAB高维转录分析代码片段
在 MATLAB 中进行高维转录分析(例如单细胞 RNA 测序 scRNA-seq 分析)可以使用以下步骤,包括数据导入、预处理、降维、聚类和差异表达分析。 MATLAB 的矩阵操作能力和工具箱(如统计和机器学习工具箱)提供了高效的分析手段。
以下是高维转录分析的步骤:
1. 数据导入
使用 readmatrix
导入基因表达数据,通常数据表包含基因名称和细胞的表达值。
matlab">% 读取数据
data = readmatrix('gene_expression.csv'); % 数据格式:行 = 基因, 列 = 细胞
genes = data(:,1); % 假设第一列为基因名称
counts_matrix = data(:,2:end); % 提取基因表达矩阵
2. 质量控制和预处理
过滤表达较低的细胞和基因以提高分析质量。
matlab">% 计算每个细胞和基因的计数
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. 归一化和对数转换
使用每个细胞的总表达值进行归一化,然后取对数转换。
matlab">% 归一化每个细胞的总计数
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 降维并可视化。
matlab">% 中心化数据
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)。
matlab">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 更适合高维数据的可视化。
matlab">% 执行 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. 差异基因表达分析
找到每个簇的标志基因,计算每个基因在不同簇中的平均表达。
matlab">% 计算每个簇的平均表达
mean_expression = zeros(size(log_normalized_counts, 1), k);
for i = 1:kmean_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. 导入必要的库
python">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
对象。
python"># 加载数据
adata = sc.read_csv('gene_expression.csv')
adata.var_names_make_unique() # 确保基因名称唯一
3. 数据预处理
包括过滤低质量的细胞和基因,进行归一化和对数转换等步骤。
python"># 基因和细胞过滤
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 图,例如每个细胞的基因数和表达水平的分布。
python"># 计算一些常用的 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),这些基因是后续降维分析的重点。
python"># 选择高变基因
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 分析,并可以选择绘制散点图。
python"># PCA 降维
sc.tl.pca(adata, svd_solver='arpack')
sc.pl.pca(adata, color='total_counts') # 使用细胞总计数上色
7. t-SNE 和 UMAP 降维
使用 t-SNE 和 UMAP 进行可视化,将数据映射到二维平面。
python"># 计算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
算法进行聚类,以区分细胞群体。
python"># 聚类分析
sc.tl.leiden(adata, resolution=0.5)
sc.pl.umap(adata, color='leiden') # 使用聚类结果上色
9. 差异基因表达分析
可以使用聚类标签进行差异表达分析,找出每个簇的特征基因。
python"># 差异表达分析
sc.tl.rank_genes_groups(adata, 'leiden', method='t-test')
sc.pl.rank_genes_groups(adata, n_genes=10, sharey=False)
10. 保存和导出分析结果
分析完成后,可以将 AnnData
对象或分析结果保存为文件。
python"># 保存数据
adata.write('processed_gene_expression.h5ad')
总结
此流程利用了 scanpy
的多种方法来进行高维转录分析,包括预处理、降维、聚类和差异表达分析等步骤。这些步骤可以帮助识别细胞类型、分群以及特征基因。