MATLAB和Python单细胞RNA测序
MATLAB
在MATLAB中处理单细胞RNA测序（scRNA-seq）数据通常涉及数据读取、预处理、降维和聚类分析等步骤。MATLAB拥有生物信息学工具箱（Bioinformatics Toolbox）和其他数据处理工具箱，可以用于单细胞RNA测序分析。以下是一些基本的步骤和示例代码，帮助您使用MATLAB进行单细胞RNA测序分析。
1. 数据读取与导入
单细胞RNA测序数据通常以不同格式保存，例如10X Genomics格式、CSV文件或HDF5文件。在MATLAB中，可以使用内置函数读取这些格式的数据。
% 示例：读取CSV格式的scRNA-seq数据
data = readtable('scRNA_data.csv');% 如果是10X Genomics的HDF5文件，可以使用h5read函数
expressionData = h5read('filtered_feature_bc_matrix.h5', '/matrix');
2. 数据预处理
预处理是单细胞分析的重要步骤，包括过滤低质量细胞、标准化和对数变换等。
2.1 过滤低质量细胞和基因
可以设置阈值去除低表达的基因和低质量细胞。
% 删除低表达基因
minCounts = 3;
minCells = 3;
filteredData = data(sum(data{:,:} > minCounts, 2) > minCells, :);% 删除低质量细胞
minGenes = 200;
filteredData = filteredData(:, sum(filteredData{:,:} > 0, 1) > minGenes);
2.2 对数标准化
对数变换通常用于降低极端值的影响，使数据更适合后续的降维和聚类分析。
% 对数据进行对数标准化处理
logNormalizedData = log1p(filteredData{:,:}); % log1p是log(1 + x)的计算方法
3. 降维处理
常见的降维方法包括PCA、tSNE和UMAP。MATLAB可以直接使用PCA和tSNE进行降维。
3.1 PCA降维
主成分分析可以将高维数据压缩为少数几个主成分，保留数据的主要信息。
% 使用PCA降维
[coeff, score, ~] = pca(logNormalizedData');
% score变量包含降维后的数据，前几个主成分可以用于后续分析
reducedData = score(:, 1:10); % 选择前10个主成分
3.2 t-SNE降维
t-SNE用于非线性降维，能够很好地展示高维数据中的簇结构。
% 使用t-SNE降维
rng('default'); % 设置随机数种子以获得一致结果
tsneData = tsne(logNormalizedData', 'NumDimensions', 2);
scatter(tsneData(:, 1), tsneData(:, 2));
title('t-SNE Visualization of scRNA-seq Data');
4. 聚类分析
聚类分析用于识别具有相似表达特征的细胞群，常见方法包括K均值聚类和层次聚类。
4.1 K均值聚类
K均值聚类是常用的聚类方法之一，可以将细胞分成若干个簇。
% 使用K均值聚类
numClusters = 5; % 设定簇的数量
clusterIdx = kmeans(reducedData, numClusters);% 可视化聚类结果
gscatter(tsneData(:, 1), tsneData(:, 2), clusterIdx);
title('K-means Clustering on t-SNE Data');
4.2 层次聚类
层次聚类是一种基于距离的聚类方法，适用于更细致的分群。
% 使用层次聚类
z = linkage(reducedData, 'ward');
dendrogram(z);
title('Hierarchical Clustering Dendrogram');
5. 差异表达分析
在细胞聚类后，可以进行差异表达分析，识别不同群体之间的差异表达基因。
% 假设已经有了聚类标签clusterIdx
geneNames = filteredData.Properties.RowNames;
clusterA = logNormalizedData(:, clusterIdx == 1); % 群体A的数据
clusterB = logNormalizedData(:, clusterIdx == 2); % 群体B的数据% t检验找到差异表达基因
pValues = zeros(size(logNormalizedData, 1), 1);
for i = 1:size(logNormalizedData, 1)[~, p] = ttest2(clusterA(i, :), clusterB(i, :));pValues(i) = p;
end% 根据p值筛选显著基因
significantGenes = geneNames(pValues < 0.05);
6. 可视化
可视化是scRNA-seq分析中不可或缺的一部分，通过不同颜色和形状来展示细胞的聚类结果和基因表达情况。
% t-SNE结果上展示聚类
figure;
gscatter(tsneData(:,1), tsneData(:,2), clusterIdx);
xlabel('t-SNE 1');
ylabel('t-SNE 2');
title('t-SNE with Cluster Labels');% 热图展示差异表达基因
significantData = logNormalizedData(ismember(geneNames, significantGenes), :);
heatmap(significantData, 'ColorLimits', [0 3]);
title('Differentially Expressed Genes Heatmap');
总结
MATLAB可以帮助分析和处理单细胞RNA测序数据，从数据导入、预处理、降维、聚类到差异表达分析和可视化。
Python
在Python中，单细胞RNA测序（scRNA-seq）数据分析可以使用多种专门的工具和库，如Scanpy、Seurat (for R, interoperable with Python) 和 Scvi-tools。其中，Scanpy是专门为Python设计的单细胞数据分析包，功能强大，特别适合处理大型scRNA-seq数据。下面将详细介绍如何在Python中使用这些工具进行单细胞RNA测序的标准分析流程。
1. 环境设置与数据导入
1.1 安装依赖库
可以通过以下命令安装主要的分析工具包：
pip install scanpy anndata matplotlib seaborn scikit-learn
1.2 导入数据
scRNA-seq数据通常存储在AnnData格式中，适合存储基因表达矩阵及其关联的元数据。如果使用10X Genomics数据，可以直接读取h5文件。
import scanpy as sc# 从10X数据中读取
adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")
# 或者从CSV/TSV文件读取
adata = sc.read_csv("scRNA_data.csv")
2. 数据预处理
预处理步骤包括筛选细胞和基因、标准化、对数变换以及高变基因选择等。
2.1 筛选低质量的细胞和基因
通常过滤低表达基因和低质量细胞以提高数据质量。
# 保留表达超过200个基因的细胞，保留在3个以上细胞中表达的基因
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
2.2 计算线粒体基因含量并过滤
高线粒体基因比例可能指示低质量细胞或凋亡细胞。
adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True)# 过滤线粒体基因含量高的细胞
adata = adata[adata.obs.pct_counts_mt < 5, :]
2.3 标准化和对数变换
将每个细胞的基因表达值归一化为总表达量的1万分之一，并进行对数变换。
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
2.4 识别高变基因
识别高变基因有助于降低分析维度并专注于生物学上有意义的基因。
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
adata = adata[:, adata.var.highly_variable]
3. 降维分析
降维分析帮助在低维空间中表示细胞间的相似性。常用的降维方法包括PCA、t-SNE和UMAP。
3.1 PCA降维
先使用PCA降低数据维度。
sc.tl.pca(adata, svd_solver="arpack")
sc.pl.pca_variance_ratio(adata, log=True)
3.2 邻近图构建
基于PCA结果构建KNN图，常用于后续的聚类和UMAP/t-SNE降维。
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)
3.3 UMAP或t-SNE降维
UMAP通常比t-SNE更快速，且能更好保留数据全局结构。
sc.tl.umap(adata)
sc.pl.umap(adata, color=["sample", "cell_type"])
4. 聚类分析
聚类分析用于识别具有相似表达谱的细胞群。Scanpy中使用基于图的Louvain或Leiden聚类算法。
sc.tl.leiden(adata, resolution=0.5) # 分辨率参数影响簇的数量
sc.pl.umap(adata, color=["leiden"])
5. 差异表达分析
通过差异表达分析找到特定簇中的标志性基因，以进一步了解每个细胞群的生物学特征。
sc.tl.rank_genes_groups(adata, "leiden", method="t-test")
sc.pl.rank_genes_groups(adata, n_genes=20, sharey=False)
6. 可视化
可视化是scRNA-seq分析中重要的一部分。可以通过UMAP/t-SNE可视化聚类结果、基因表达情况等。
# 查看特定基因在UMAP图中的表达分布
sc.pl.umap(adata, color=["GeneA", "GeneB"])# 绘制热图展示不同簇的差异表达基因
sc.pl.rank_genes_groups_heatmap(adata, n_genes=10, groupby="leiden")
7. 保存和导出数据
可以将分析结果保存为h5ad文件，以便后续分析或共享。
adata.write("processed_scRNA_data.h5ad")
完整代码示例
以下是从导入到聚类和可视化的完整代码示例：
import scanpy as sc# 导入数据
adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")# 数据预处理
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True)
adata = adata[adata.obs.pct_counts_mt < 5, :]
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
adata = adata[:, adata.var.highly_variable]# 降维与聚类
sc.tl.pca(adata, svd_solver="arpack")
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=0.5)# 差异表达分析
sc.tl.rank_genes_groups(adata, "leiden", method="t-test")# 可视化
sc.pl.umap(adata, color=["leiden"])
sc.pl.rank_genes_groups_heatmap(adata, n_genes=10, groupby="leiden")# 保存结果
adata.write("processed_scRNA_data.h5ad")
总结
Python提供了丰富的工具来处理和分析单细胞RNA测序数据，尤其是Scanpy，可以高效地执行数据预处理、降维、聚类和差异表达分析。