数据分析思维之分而治之

在前面的笔记中： 扎克伯格背刺基于R语言的Seurat单细胞生态 我提到了一个比较大数据量（接近100万个细胞啦）的单细胞转录组项目，是 929,686 cells derived from 156 fresh clinical samples obtained from 41 HGSOC patients 。

这样的话，它是 31815个基因 ，是 929690 个细胞，所以数值会很恐怖，是29578087350，但是因为是稀疏矩阵，所以这个单细胞表达量矩阵里面绝大部分都是0值，真正有数值的地方是2178171554，差不多是7.36%，这个是单细胞转录组的特性：drop-out。如下所示：

> 2178171554/29578087350
[1] 0.07364139

既然数值的地方是2178171554，那么在R编程语言里面读取sparse_matrix.mtx文件会 报错：

Error in scan(file, nmax = 1, what = what, quiet = TRUE, ...) : 
  scan() expected 'an integer', got '2178171554'
Calls: readMM -> scan1 -> scan
Execution halted

==> inputs/sparse_matrix.mtx <==
%%MatrixMarket matrix coordinate real general
%
31815 929690 2178171554

因为上面的数值 2178171554 确实超出了 R 中 int 类型的范围。在 32 位系统中，R 的 int 类型通常是一个 32 位整数，其范围是 -2^31 到 2^31-1，即 -2,147,483,648 到 2,147,483,647。

但是经费是有限的，不可能说是无限制添加内存条，所以我提出来了： 在R编程环节有所限制未必不是好事 。我找人工智能大模型帮我写了一个代码把上面的单细胞表达量矩阵拆分成为两个：

import scanpy as sc
import numpy as np
import scipy.io as sio

# 读取单细胞数据
all_data = sc.read_h5ad("./fdff1375-aafc-48ef-b5b1-b1ff8466d92c.h5ad")

# 随机拆分成两个部分
np.random.seed(42)  # 设置随机种子以确保可重复性
indices = np.random.permutation(all_data.shape[0])
split_index = indices[int(0.5 * len(indices))]  # 找到中间索引进行拆分

# 创建两个子集
subset1 = all_data[indices[:split_index], :]
subset2 = all_data[indices[split_index:], :]

# 为每个子集导出数据
def export_subset(subset, subset_name):
    # 导出细胞信息
    subset.obs.to_csv(f"{subset_name}_cellinfo.csv")
    # 导出基因信息
    subset.var.to_csv(f"{subset_name}_geneinfo.csv")
    # 导出表达矩阵
    raw_data = subset.raw.X
    # 转置矩阵，因为通常细胞是行，基因是列
    sparse_matrix = raw_data.T 
    sio.mmwrite(f"{subset_name}_sparse_matrix.mtx", sparse_matrix)

# 导出第一个子集的数据
export_subset(subset1, 'subset1')

# 导出第二个子集的数据
export_subset(subset2, 'subset2')

关键的Python代码是使用了 numpy 库来随机拆分一个矩阵（或数组）成两个子集。下面是代码的逐行解释：

1. np.random.seed(42) ：这行代码设置了随机数生成器的种子为42。设置随机种子可以确保每次运行代码时生成的随机数序列是相同的，这有助于结果的可重复性。

2. indices = np.random.permutation(all_data.shape[0]) ：这行代码生成了一个随机排列的索引数组。 all_data.shape[0] 给出了 all_data 矩阵的行数， np.random.permutation 函数对这个数字进行随机排列，生成一个从0到 all_data 行数减1的随机排列数组。

3. split_index = indices[int(0.5 * len(indices))] ：这行代码计算了拆分索引，即在随机排列的索引数组中间位置的值。 int(0.5 * len(indices)) 计算了数组长度的一半（向下取整），然后使用这个值从 indices 数组中取出一个索引，这个索引将用来将 all_data 矩阵拆分成两个大致相等的部分。

4. subset1 = all_data[indices[:split_index], :] ：这行代码创建了第一个子集 subset1 。它使用切片 indices[:split_index] 来选择 all_data 矩阵的前半部分行（根据随机排列的索引）， : 表示选择所有的列。因此， subset1 包含了 all_data 矩阵中随机选择的前半部分行。

5. subset2 = all_data[indices[split_index:], :] ：这行代码创建了第二个子集 subset2 。它使用切片 indices[split_index:] 来选择 all_data 矩阵从 split_index 开始的后半部分行，同样 : 表示选择所有的列。因此， subset2 包含了 all_data 矩阵中随机选择的后半部分行。

总的来说，这段代码将 all_data 矩阵随机拆分成两个子集 subset1 和 subset2 ，每个子集包含 all_data 的一半行，且行的选择是基于随机排列的索引。这种方法常用于机器学习中的数据集拆分，例如将数据集拆分成训练集和测试集。

之前没有拆分的时候得到的单细胞表达量矩阵会肉眼看起来就很多：

ls -lh inputs/
total 55G
-rw-rw-r-- 1 t180559 t180559 537M 11月  9 23:44 cellinfo.csv
-rw-rw-r-- 1 t180559 t180559 2.1M 11月  9 23:44 geneinfo.csv
-rw-rw-r-- 1 t180559 t180559  54G 11月 10 00:52 sparse_matrix.mtx

现在拆分成为了两个：

$ ls -lh subset*
-rw-rw-r-- 1 t180559 t180559 332M 11月 10 18:00 subset1_cellinfo.csv
-rw-rw-r-- 1 t180559 t180559 2.1M 11月 10 18:00 subset1_geneinfo.csv
-rw-rw-r-- 1 t180559 t180559  34G 11月 10 18:30 subset1_sparse_matrix.mtx
-rw-rw-r-- 1 t180559 t180559 205M 11月 10 18:31 subset2_cellinfo.csv
-rw-rw-r-- 1 t180559 t180559 2.1M 11月 10 18:31 subset2_geneinfo.csv
-rw-rw-r-- 1 t180559 t180559  21G 11月 10 18:48 subset2_sparse_matrix.mtx

==> subset1_sparse_matrix.mtx <==
%%MatrixMarket matrix coordinate real general
%
31815 575647 1348277546 
==> subset2_sparse_matrix.mtx <==
%%MatrixMarket matrix coordinate real general
%
31815 354043 829894008

可以看到，之前的929690细胞是同一个单细胞转录组数据集，然后被拆分成为了 575647 和 354043个细胞数量的两个项目，不知道为什么不是平均分配。

拆分成为了两个项目就可以各自独立的跑单细胞转录组数据分析的降维聚类分群流程啦，而且也可以很容易给出来生物学名字，如下所示：

但是呢，这个是因为是跑了两次单细胞转录组流程，所以后面的分析比如如果是要把里面的内皮细胞子集提取出来做细分的时候，就可以把两个对象里面的内皮细胞合并，然后继续走单细胞转录组流程。 也就是说， 除了第一层次降维聚类分群是需要这个项目里面的全部的单细胞之外，后面的分析都可以细化到具体的细胞亚群里面去。