单细胞分析过程中的稀疏矩阵删减

.
├── main.py
├── data
│   ├── barcodes.tsv
│   ├── features.tsv
│   ├── matrix.mtx.gz
│   ├── selected.tsv
├── results
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   └── matrix.mtx.gz

from scipy.io import mmread
import pandas as pd
import numpy as np

# 读取表达矩阵
_index = pd.read_csv("./data/features.tsv", index_col=1, sep="\t", header=None)
_index.index.name = None  # 把索引列的列名去掉

_col = pd.read_csv("data/barcodes.tsv", sep="\t", header=None)
_col.index.name = None  # 把列名向量的名去掉

_data = mmread("data/matrix.mtx.gz").todense()

# 处理表达矩阵
## 挑选需求 col
_col.index
_selected = pd.read_csv("./data/selected.tsv", sep="\t", header=None)
_selected.index.name = None  # 把索引列的列名去掉
filtered_index = list(set(_col[0]).intersection(_selected[0]))
print(filtered_index.__len__())

# 加行名列名
rna_count = pd.DataFrame(
    data=_data, index=_index.index, columns=_col.iloc[:, 0].tolist()
)
filtered_columns = rna_count[filtered_index]
rna_count = filtered_columns
print(rna_count.iloc[0:3, 0:2])
print("gene_ID_len : " + str(rna_count.shape[0]))  ### 获取表达矩阵基因数
print("cell_ID_len : " + str(rna_count.shape[1]))  ### 获取表达矩阵细胞数

# 重新写出 DataFrame 为 10X 格式的 sparse matrix 等相关文件
import os
import shutil
import gzip
import scipy
import time

fmt = "%Y-%m-%d %a %H:%M:%S"
Date = time.strftime(fmt, time.localtime(time.time()))
outdir = "results"
os.makedirs(outdir, exist_ok=True)

# Save matrix.mtx.gz
reAnno_count_sparse_mtx = scipy.sparse.coo_matrix(rna_count.values)
scipy.io.mmwrite(
    os.path.join(outdir, "matrix.mtx"),
    reAnno_count_sparse_mtx,
    comment="Generate DateTime::" + str(Date),
)
with open(os.path.join(outdir, "matrix.mtx"), "rb") as mtx_in:
    with gzip.open(
        os.path.join(outdir, "matrix.mtx") + ".gz", "wb"
    ) as mtx_gz:  # 创建一个读写文件'matrix.mtx.gz'，用以将 matrix.mtx 拷贝过去
        shutil.copyfileobj(mtx_in, mtx_gz)
os.remove(os.path.join(outdir, "matrix.mtx"))

# Save barcodes.tsv.gz
barcodesFile = pd.DataFrame(rna_count.columns)
barcodesFile.to_csv(
    os.path.join(outdir, "barcodes.tsv.gz"), sep="\t", header=False, index=False
)

# Save features.tsv.gz
featuresFile = _index
featuresFile.to_csv(
    os.path.join(outdir, "features.tsv.gz"), sep="\t", header=False, index=False
)

numpy==1.24.3
pandas==2.0.1
scipy==1.11.4