一、矩阵处理
- 情况一：常规矩阵
  - 例子1：常规处理
  - 补充例子1：最常用的转换获取TPM的操作
- 情况二：稀疏矩阵
二、cellAnnotion处理
- 例子1: cellID刚好能够和title对应
- 例子2: 存在某个线索可以实现cellID和sampleID的转化
三、快速填写marker基因
- 情况1： markerGene文件以一个表格的形式存放
- 情况2： markerGene文件以多重Excel表格的形式存放,不同的sheet用对应的clusterName表示

所有代码只做参考，大家需要根据自己情况修改

一、矩阵处理

情况一：常规矩阵

例子1：常规处理

最常见的矩阵形式是行为基因，列为细胞，其中一般第一列为geneName。当然也有前多列为gene信息的情况, 如下：

处理思路：

去除矩阵中非数据列（通常是geneName）

矩阵转置

依次添加cellID和normalizationMethod

当然也可以先转置后去除非数据行，效果是一样的。

#读取矩阵
df_downloaded = pd.read_csv('../downloaded_data/GSE_supplementary_data/GSE130812_FPKM_C1_table.txt.gz',sep='\t')
df_downloaded.head()

##单独提取基因映射文件（前两列为基因信息列，列名为'gene_id'和'gene_name'），用于后面矩阵转置后的列重命名；
##小建议：一般建议使用ensemblID（形如ENSMUSG00000000001）来标识基因，这样不会出现基因名字重复的问题。
geneMap = df_downloaded.iloc[:,[0,1]] 
df_downloaded = df_downloaded.drop(['gene_id','gene_name'],1) ##去除矩阵中的非数据列；便于转置
###有时候读取的文件中gene不在第一列，而是作为行名出现，此时可以忽略此步操作。 当然，你也可以使用gene=df_downloaded.index.tolist()提取基因

##提取cellID
cellID = df_downloaded.columns.tolist()
df_norm = df_downloaded.T ##最常见的矩阵形式是行为基因，列为细胞，所以必须转置
df_norm.rename(columns=geneMap.iloc[:,0].str.replace('(.[0-9]+$)',""), inplace = True)##用ensemblID重命名列（即第一步提取的geneMap文件）
###ensemblID有时候会在后面带上版本号（如ENSMUSG00000000001.4中的.4）,这些版本号必须去除。这里使用正则匹配的方式直接将末尾的版本号替换成了空
##很多pandas中类似的字符提取的操作建议使用正则表达式来完成，速度较快；
df_norm.insert(0,'cellID',cellID)#插入cellID
df_norm.insert(0,'normalizationMethod','FPKM') ##插入normalizationMethod
df_norm.head()

补充例子1：最常用的转换获取TPM的操作

FPKM 转 TPM

df_tpm = df_norm ##直接将df_norm赋给df_tpm. (注意：运行完这步之后千万别再保存df_norm到文件中，因为后续对df_tpm的操作同时也会修改df_norm, 保存的话会将TPM矩阵覆盖至原本已经保存的norm文件中)
df_tpm.iloc[:,2:] = df_tpm.iloc[:,2:]/np.sum(df_tpm.iloc[:,2:],axis=1).values.reshape(-1,1)*1e6
df_tpm['normalizationMethod']='TPM from FPKM data'
df_tpm.head()

log2(TPM+1) 转 TPM

df_tpm = df_norm ##直接将df_norm赋给df_tpm. (注意：运行完这步之后千万别再保存df_norm到文件中，因为后续对df_tpm的操作同时也会修改df_norm, 保存的话会将TPM矩阵覆盖至原本已经保存的norm文件中)
df_tpm.iloc[:,2:] = df_tpm.iloc[:,2:].apply(np.exp2)-1 ##如果是ln(TPM+1) --> TPM,则直接将.apply(np.exp2)-1改为.apply(np.expm1)
df_tpm['normalizationMethod']='TPM from log2(TPM+1) data'
df_tpm.head()

log2(FPKM+1) 转 TPM

df_tpm = df_norm##直接将df_norm赋给df_tpm. (注意：运行完这步之后千万别再保存df_norm到文件中，因为后续对df_tpm的操作同时也会修改df_norm, 保存的话会将TPM矩阵覆盖至原本已经保存的norm文件中)
df_tpm.iloc[:,2:] = df_tpm.iloc[:,2:].apply(np.exp2)-1
df_tpm.iloc[:,2:] = df_tpm.iloc[:,2:]/np.sum(df_tpm.iloc[:,2:],axis=1).values.reshape(-1,1)*1e6
df_tpm['normalizationMethod']='TPM from log2(FPKM+1) data'

情况二：稀疏矩阵

关于何为稀疏矩阵可以参考以下博文：python scipy 稀疏矩阵详解

一般细胞数在1万以上就需要存成稀疏矩阵，这样方便后续的操作，因为后面绝大多数的自动函数都需要读取TPM矩阵，保存常规矩阵会运行得特别慢。按照公司的要求低于1万的数据尽量保存成常规矩阵，除非数据很大严重影响到了你后续的操作（有的时候细胞数稍低也有数据偏大的情况）

情况2.1: 单个稀疏矩阵

建议直接读取稀疏矩阵后就保存。也可以直接对稀疏矩阵进行操作（矩阵切片提取部分细胞、FPKM转TPM等）后再保存

例子1: 数据不需要做任何操作，直接读取然后保存

from scipy import sparse
from scipy.sparse import coo_matrix
from scipy.io import mmread, mmwrite, mminfo
from scipy.sparse import coo_matrix, hstack,vstack
##读取稀疏矩阵三个关键文件
cellID_raw= pd.read_csv('../downloaded_data/E-MTAB-8483.aggregated_filtered_counts.mtx_cols',header=None,names=['cellID'])
genes_raw  = pd.read_csv('../downloaded_data/E-MTAB-8483.aggregated_filtered_counts.mtx_rows',sep='\t',header=None,names=['ENSMUG1','ENSMUG2'])  ##这里需要先读取基因文件查看文件格式，然后视情况自行对names进行修改
coo_mtx_raw =mmread('../downloaded_data/E-MTAB-8483.aggregated_filtered_counts.mtx')

##单独提取cellID和gene
genes = genes_raw['ENSMUG1'].tolist()
cellID = cellID_raw['cellID'].tolist()

#保存结果。 读取的矩阵一般都是行为基因，列为样本。咱们保存的正好相反。所以需要转置
my_builder.save_template_mtx(mtx=coo_mtx_raw.T, mtx_name='expressionMatrix_rawCounts.mtx', genes=genes, cellID=cellID)

例子2: 需要提取部分细胞（直接对稀疏矩阵进行细胞筛选）

from scipy import sparse
from scipy.sparse import coo_matrix
from scipy.io import mmread, mmwrite, mminfo
from scipy.sparse import coo_matrix, hstack,vstack
##读取稀疏矩阵三个关键文件
cellID_raw= pd.read_csv('../downloaded_data/E-MTAB-8483.aggregated_filtered_counts.mtx_cols',header=None,names=['cellID'])
genes_raw  = pd.read_csv('../downloaded_data/E-MTAB-8483.aggregated_filtered_counts.mtx_rows',sep='\t',header=None,names=['ENSMUG1','ENSMUG2'])  ##这里需要先读取基因文件查看文件格式，然后视情况自行对names进行修改
coo_mtx_raw =mmread('../downloaded_data/E-MTAB-8483.aggregated_filtered_counts.mtx')
genes = genes_raw['ENSMUG1'].tolist()

####直接对稀疏矩阵进行细胞筛选（例如我们只提取第55列到61列进行保存）
cellID = cellID_raw['cellID'].tolist()
filter_list = cellID[56:60]
###提权被筛选细胞的索引
cellID_index = [cellID.index(i) for i in cellID if i in filter_list]
##获得筛选后的cellID
cellID_filter = [cellID[i] for i in cellID_index]	
##基于cellID的索引对提取稀疏矩阵中的子列
coo_mtx_raw_filter = coo_mtx_raw.tocsr()[:,cellID_index]	
##提取的稀疏矩阵，type 为scipy.sparse.csr.csr_matrix，需要需要转化一些方便储存
coo_mtx_raw_filter = sparse.coo_matrix(coo_mtx_raw_filter)

##保存结果。 读取的矩阵一般都是行为基因，列为样本。咱们保存的正好相反。所以需要转置
my_builder.save_template_mtx(mtx=coo_mtx_raw_filter.T, mtx_name='expressionMatrix_rawCounts.mtx', genes=genes, cellID=cellID_filter)

情况2.2: 多个稀疏矩阵直接进行拼接再保存

稀疏矩阵越多越有必要使用此种方式处理。常规的转换成常规矩阵再拼接，内存压力太大，速度太慢容易崩溃。（5星推荐）

例子如下：整个过程分三段

循环匹配文件名并对文件名排序

基于提取的文件名循环读取稀疏矩阵获取稀疏矩阵list

合并所有稀疏矩阵

使用这段代码之前建议先简单查看一下各个文件：

查看barcode（判断是否需要添加GSM号用于区分不同文件中的barcode，这个是为了方便后面cellAnnotion中cellID和sample的匹配）

查看gene格式（明确哪一列是geneName，哪一列是ensemblID,以及ensemblID是否带有版本号；更重要的是查看所有样本的gene文件是否相同，不相同没法简单合并，得做基因提取，不过这种情况极少出现）

使用shape查看mtx的行列数，以此与对应barcode和gene文件对比，判断是否需要对mtx进行转置（目前见过的都需要转置）

###循环匹配文件，稀疏矩阵读取
import fnmatch

filetxt_genes = []
filetxt_barcodes = []
filetxt_matrix = []
root_dir = '../downloaded_data/GSE_supplementary_data/GSE111360_RAW/'#此处的路径要做对应的修改
for root,dirs,files in os.walk(root_dir):
    for file in files:
        if fnmatch.fnmatch(file, '*genes.tsv.gz'):###匹配对应的gene文件，看情况修改
            filetxt_genes.append(file)
        if fnmatch.fnmatch(file, '*barcodes.tsv.gz'):###匹配对应的barcodes文件，看情况修改
            filetxt_barcodes.append(file)
        if fnmatch.fnmatch(file, '*matrix.mtx.gz'):###匹配对应的matrix文件，看情况修改
            filetxt_matrix.append(file)

###排序文件，保持文件一致（这段代码非常重要，不进行排序可能直接导致barcode和矩阵数据不对应，到时候错都不知道错在哪）
filetxt_genes.sort()
filetxt_genes

filetxt_barcodes.sort()
filetxt_barcodes

filetxt_matrix.sort()
filetxt_matrix

###循环读取稀疏矩阵，获取稀疏矩阵list  dataframe_list 
from scipy import sparse
from scipy.sparse import coo_matrix
from scipy.io import mmread, mmwrite, mminfo
from scipy.sparse import coo_matrix, hstack,vstack
import re
dataframe_list = []
cellID = []
for i in range(0,len(filetxt_barcodes)):
    print(i)
    path_cellfile= root_dir + filetxt_barcodes[i]
    path_genefile= root_dir + filetxt_genes[i]
    path_matrix= root_dir + filetxt_matrix[i]
    df1 = pd.read_csv(path_cellfile ,header=None,names=['cellID'])
    print("Read cellFile!!")
    barcodes = df1['cellID']
    barcodes = (re.sub('_.+$','_',filetxt_barcodes[i]) + barcodes).tolist()###这里是为了区分样本的cellID,便于后期cellAnnotion的处理，需要看情况修改
    cellID.extend(barcodes)
    df2 = pd.read_csv(path_genefile ,sep="\t",header=None,names=['ensemblID','geneSymbol'])###先查看一下gene文件的格式，如果第一列是ensemblID，第二列是geneSymbol就不用修改
    print("Read geneFile!!")
    gene = df2['ensemblID'] ##根据df2的格式修改
    gene = gene.tolist()
    print(len(gene)) ##根据输出信息，查看所有样本的gene数是否相同，相同的话则基本说明他们的基因顺序是一致的，可以直接合并；不一致则需要做基因提取，参照情况2.1的例子2
    coo_mtx = mmread(path_matrix)
    print("Read matrix!!")
    dataframe_list.append(coo_mtx.T)##千万别忘记转置，目前我们碰到的稀疏矩阵都需要转置

##直接将dataframe_list中的稀疏矩阵合并
mtx_data = vstack(dataframe_list)

#然后用脚本上的代码保存mtx_data、gene、cellID就好了
my_builder.save_template_mtx(mtx=mtx_data, mtx_name='expressionMatrix_rawCounts.mtx', genes=gene, cellID=cellID)

补充例子2：常规矩阵和稀疏矩阵的转换

dataframe 转 coo_matrix

这段代码在多个矩阵合并之后细胞数过多时使用（一般是细胞数超过1万）

###读取表达谱
df_downloaded = pd.read_csv('../downloaded_data/GSE_supplementary_data/GSE130812_FPKM_C1_table.txt.gz',sep='\t')
df_downloaded.head()

####因为咱们的稀疏矩阵中只保留表达值，cellID和geneName需要单独提出保存在其他文件里，所以得事先将它们提取出来，并删除矩阵中的非数据列
geneMap = df_downloaded.iloc[:,[0,1]]
df_downloaded = df_downloaded.drop(['gene_id','gene_name'],1)
genes=geneMap.iloc[:,0].str.replace('.[0-9]+$',"").tolist()##使用ensemblID表示基因
###ensemblID有时候会在后面带上版本号（如ENSMUSG00000000001.4中的.4）,这些版本号必须去除。这里使用正则匹配的方式直接将末尾的版本号替换成了空

#提取cellID
cellID = df_downloaded.columns.tolist()

df_raw = df_downloaded.T##矩阵中行为基因，列为样本，和我们储存的要求刚好相反，需要转置

##dataframe转coo_matrix
from scipy import sparse, io
mtx_raw = sparse.coo_matrix(df_raw)

#最后使用脚本里的代码保存稀疏矩阵就完事了
my_builder.save_template_mtx(mtx=mtx_raw , mtx_name='expressionMatrix_rawCounts.mtx', genes=genes, cellID=cellID )

coo_matrix 转 dataframe

建议在处理时先查看所有barcodes文件合并后的细胞数是否超过1万，超过则不用转成dataframe，此时可以直接使用情况二中的代码；

from scipy import sparse
from scipy.sparse import coo_matrix
from scipy.io import mmread, mmwrite, mminfo

#读cellID文件
df1 = pd.read_csv('../downloaded_data/GSM3407039_Jones-OM1_barcodes.tsv.gz',header=None,names=['cellID'])
#提取cellID为list
cellID = df1['cellID']
cellID = cellID.tolist()
cellID
#读geneSymbol或者ensemblID文件
df2 = pd.read_csv('../downloaded_data/GSM3407039_Jones-OM1_genes.tsv.gz',sep="\t",header=None,names=['ensemblID','geneSymbol'])
#提取ensemblID为list
gene = df2['ensemblID']

#读稀疏矩阵
coo_mtx = mmread('../downloaded_data/GSM3407039_Jones-OM1_matrix.mtx.gz')
###将稀疏矩阵转化为dataframe(同上，稀疏矩阵通常也是行为基因列为样本，所以需要转置)
coo = coo_mtx.T.toarray()###转置，并转化为array
df_raw = pd.DataFrame(coo)##将array转为dataframe
df_raw.head()
#把列名columns换为gene
df_raw.rename(columns=gene, inplace = True)
#插入cellID列
df_raw.insert(0,'cellID',cellID)
###df_raw.insert(0,'normalizationMethod','FPKM') ##如果这里处理的是标化矩阵和TPM矩阵记得添加normalizationMethod

#这样就转化完毕了

补充例子3：如何处理HDF5(.h5)格式文件

h5 文件是以固定的方式将一组文件 (样本测序相关的 metadata、matrix、cluster 信息、坐标，等等）. 打包成为 AnnData 用于 Scanpy 处理流程。AnnData class 包括一个 object adata 储存 data matrix, 即 adata.X，一个 annotation of observations, 即 adata.obs，一个 annotation of variables, 即 adata.var (pd.DataFrame)，以及一个 unstructured annotation, 即 adata.uns (dict)。详情参见

个人理解：h5 文件可以看做是一个包含了多个文件或文件夹的压缩文件，打开文件时同时创建一个嵌套句柄，连接h5中的所有文件；文件夹对应一个group（所以需要用.keys()去获取其中的文件名，意思是需要嵌套的句柄指向该文件夹内的文件）；而文件对应group中的子项（可以直接用.value（注意，现在换成 [()] 的方式取值了）去读取其中的值，即句柄直接指向文件读取）； 同样，也可以理解为一个多层嵌套的文件字典

使用 h5 文件的第一步是将其打开，检查数据内容，并将其中储存成 array 形式的数据整理成矩阵。已打开的文件不可重复打开，需要先 close，如果 close 不成功... 当然也可以直接 Shotdown 重来。详情参见

例子代码如下：

##h5文件的形式多种多样，请根据数据的情况修改，关键是要屡清楚文件中的各个数据长啥样；

##导入相关模块
import h5py  
from scipy import sparse
import numpy as np
from scipy.sparse import csr_matrix
from scipy import sparse, io
import pandas as pd


# 读入hdf5格式文件
f = h5py.File('./miller20-Copy1.processed.h5ad','r+')


# 查看根下的所有group
list(f.keys())  
# Output eg: ['X', 'obs', 'obsm', 'uns', 'var']

#### 一般 X、uns 是 group，obs、obsm、var 是 dataset。
### 同时，X 是 mtx，不带 cellID 和 geneSymbol 的测序结果；obs 是测序结果汇总、cluster 信息和一些 cell metadata；obsm 是 tSNE、UMAP 坐标；uns 是 sample metadata；var 基因相关内容。
##此例子中没有明显的cellID和geneSymbol,但是查看obs和var能够发现，cellID和geneSymbol分别在obs和var中；

###这些内容不是固定的，要根据数据信息自行判定；（把数据自己输出来看看）

###获取cell mate信息
##obs, obsm
cell_mate = f['obs'][()]
cell_mate = pd.DataFrame(cell_mate)
cell_mate.head()

#obsm中存放的tSNE或UMAP信息

##储存形式有点奇怪，不过可以通过以下方式转化
tSNE_array = np.array(f['obsm'][()].tolist())
cell_tSNE = pd.DataFrame({'tSNE1': tSNE_array[:, 0][:,0],'tSNE2': tSNE_array[:, 0][:,1], 'UMAP1': tSNE_array[:, 1][:,0], 'UMAP2': tSNE_array[:, 1][:,1]})
cell_tSNE.insert(0,'cellID',cell_mate['index'].tolist())#插入cellID
cell_tSNE.head()

##提取基因信息
gene_df = f['var'][()]
gene_df = pd.DataFrame(gene_df)
gene_df.head()

###提取表达谱
##表达谱以稀疏矩阵的形式存放
list(f['X'].keys())
##output： ['data', 'indices', 'indptr']
data = f['X']['data']
indices = f['X']['indices']
indptr = f['X']['indptr']

matrix = csr_matrix((data, indices, indptr))
matrix.shape
##output: (23832, 36807)  ##查看矩阵大小，如果细胞数不超过1万，则需要转化为DataFrame；

##barcode
cellID = cell_mate['index'].tolist()

##geneSymbol
genes = gene_df['index'].tolist()

# 并且，根据 matrix的shape来判断是否需要矩阵转置
matrix = matrix.T ##

# 转为coo矩阵
matrix = matrix.tocoo()
##之后使用模板中的代码保存稀疏矩阵即可

###当然也可将它转为DataFrame来保存；视情况而定
matrix.columns= [genes]##列重命名
matrix.insert(0,'cellID',cellID)#插入cellID
matrix.head()

###自己查看uns数据，然后提取其中有用信息
f['uns'].keys()#查看键
f['uns']['Cell_group_categories'][()]##提取值
##也可写个函数将它转成常规字典

注意: 这里的f['uns']中提供的数据是需要和前面f['obs']中提供的cell_meta结合使用的;
例如:cell_meta中的Tissue列中的0,1,2 分别对应f['uns']['Tissue_categories'][()]中的Airway，Distal和Trachea;
同理：Cell_group_categories, Cell_type_categories, phase_categories, sample_categories, orig.ident_categories, 等都可以和cell_meta中的列相对应;
所以，请务必结合两者补全cell_meta中的内容

.h5ad文件处理进阶版（更简单）

前面这些代码可能比较混乱，结构不清晰，但是可以处理多数.h5文件；以下介绍基于scanpy模块的简易处理方式，主要针对.h5ad文件，此方式能够自动将uns中的categories对应到obs中，用起来相对更简单；

import importlib.util
import sys
import os 
import pandas as pd
import numpy as np
import scanpy as sp
import h5py
import scipy
import math
from scipy import sparse

#导入模块,关键是需要导入scanpy

##使用scanpy自带的h5读取函数转化对象
adata=sp.read_h5ad('./GSM4555887_Patient_A-Copy1.h5ad')
adata

#output: #adata中包含以下对象，（表达谱X未显示出来）
#  AnnData object with n_obs × n_vars = 3553 × 18124
#    obs: 'Patient', 'batch', 'n_counts', 'n_genes', 'pheno', 'labels'
#    var: 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'
#    uns: 'Patient_colors', 'dendrogram_labels', 'labels_colors', 'neighbors', 'pca', 'pheno2_colors', 'pheno_colors', 'sample_colors'
#    obsm: 'X_pca', 'X_umap'
#    varm: 'PCs'
#    layers: 'magic'
#    obsp: 'distances', 'connectivities'

##

最后直接使用如下语句直接读取相应的对象。同样需要先观察数据，根据数据的形式应变处理

#提取表达谱
adata.X  ##表达谱，无cellID和geneName
adata.obs_names   ##cellID
adata.var_names   ##geneName

##有的时候会是稀疏矩阵的形式，根据情况自己修改代码

如下代码整理成dataFrame

##先查看数据的行列情况
adata.obs_names.shape,adata.var_names.shape,adata.X.shape
#output
((3553,), (18124,), (3553, 18124))
##行为样本直接转化，列为基因，故columns=var_names，index=obs_name
matrix = pd.DataFrame(columns=adata.var_names.tolist(),index=adata.obs_names.tolist(),data=adata.X)
matrix.head()

#数据比较大的话直接使用 sparse.csc_matrix(matrix)转化成稀疏矩阵即可

##其他相关信息
adata.obs  ##细胞注释矩阵 
adata.obsm  ##细胞相关的补充信息  （这里提供的数据不是固定的，不同的研究有可能有所不同）
adata.obsm['X_pca']  ###每个细胞对应的pca值
adata.obsm['X_umap']  ##每个细胞对应的umap坐标值

adata.var  ##基因注释矩阵
adata.varm  ## 基因相关的补充信息  （这里提供的数据不是固定的，不同的研究有可能有所不同）
adata.varm['PCs']  ##每个基因在每个主成分对应的系数（权重）

adata.uns  ##以dict形式给出所有非结构化的信息，用于对应obs或var中的categories；这个基本上可以忽略了，因为差不多的信息已经自动补全到以上对象中了

二、cellAnnotion处理

此部分处理的关键是找到cellID和sample中某一列的对应关系（一般是title）

例子1: cellID刚好能够和title对应

此时基于两者完成两个dataframe的匹配即可；

匹配不建议使用loc加循环的方式去做，速度特别慢，这里推荐使用pd.merge实现df_cell和sample的匹配

##初始化df_cell矩阵，这里就是模板中的代码
df_cell = my_builder.read_template_tsv(cellAnnotation)
df_tpm = my_builder.read_template_tsv(expressionMatrix_TPM)
df_tpm.head()
df_cell['cellID'] = df_tpm['cellID'].tolist()
df_cell.head()

###读取sample
sample = my_builder.sample_info(GSE = 'GSE130812')
sample.iloc[0] ##使用这种方式去查看sample中的条目，便于匹配，也便于后面条目的填写

关键步骤，使用使用pd.merge实现df_cell和sample的匹配

##重点：使用pd.merge实现df_cell和sample的匹配（这里使用最简单的情况：cellID和sample中的title能够匹配）
df_cell_tmp = pd.merge(df_cell,sample, how = 'left',left_on=['cellID'], right_on=['title'])  
df_cell_tmp.head()###df_cell和sample合并，创建df_cell_tmp。
#df_cell_tmp中包含了df_cell和sample中所有的列，其cellID顺序和df_cell中的一致。

#这时我们只要将df_cell_tmp中的条目直接填入df_cell中即可，如下：
df_cell['meta_sourceName'] = df_cell_tmp['source_name_ch1'].tolist()
df_cell['meta_tissue'] = df_cell_tmp['characteristics_ch1.2.tissue'].tolist()
df_cell['meta_genotype'] = df_cell_tmp['characteristics_ch1.1.genotype'].tolist()

df_cell['clusterID'] = 'notAvailable'
df_cell['clusterName'] = 'notAvailable'
df_cell['sampleID'] = df_cell_tmp['geo_accession'].tolist()
df_cell['cellOntologyName'] = 'pro-T cell'
df_cell['cellOntologyID'] = 'CL:0000827'
df_cell['FACSMarker'] = df_cell_tmp['characteristics_ch1.3.facs purification gate'].tolist()
df_cell['tSNE1'] = ''# 从文章提供的表格或者其他数据来源中查找tSNE信息填入，若无再使用6中的函数生成
df_cell['tSNE2'] = ''

补充说明：

pd.merge中的how参数推荐使用left（df_cell在左边）。原因是pd.merge会根据how的参数对合并的结果进行自动排序。如果how='left'则合并结果会与df_cell的原始顺序一致，不会打乱顺序，也不会丢失结果。即使在df_cell中有个别cellID与sample无法匹配，这参数也能保证cellID位置不变（只不过合并结果的其它列都为NAN）
当cellID和sample中的列无法直接匹配时就得自己想办法处理了，这里千万要仔细。多数情况是能找到匹配线索的，此时可以在sample中添加一列，将线索转化为一列信息然后用于合并就行（请看例子2）
同理，如果有cluster文件，也可以使用此方法将df_cell和cluster文件进行匹配

例子2: 存在某个线索可以实现cellID和sampleID的转化

如下图：已知细胞的barcode（cellID）的尾标和sample是相互对应的，那么我们可直接将cellID转化为sampleID，然后基于sampleID将df_cell和sample合并；

##初始化df_cell矩阵，这里就是模板中的代码
df_cell = my_builder.read_template_tsv(cellAnnotation)
df_tpm = my_builder.read_template_tsv(expressionMatrix_TPM)
df_tpm.head()
df_cell['cellID'] = df_tpm['cellID'].tolist()
df_cell.head()
#读取sample文件
sample = my_builder.sample_info(GSE = 'GSE142541')

使用如下代码直接根据barcode尾标将cellID转化为sampleID

df_cell['sampleID'] = df_cell['cellID'].str.replace('([ATCG]+-1)','GSM4231663') ##这段代码是正则匹配替换，意思是：将cellID中任意多个A/T/C/G四种字符的组合加-1的字符串替换成GSM4231663;并且保存到df_cell['sampleID'] 中；后几句同理；
df_cell['sampleID'] = df_cell['sampleID'].str.replace('([ATCG]+-2)','GSM4231664')##注意这里开始就是对df_cell['sampleID']进行操作了
df_cell['sampleID'] = df_cell['sampleID'].str.replace('([ATCG]+-3)','GSM4231665')
##以上代码也可以简化为一句代码，如下：
#df_cell['sampleID'] = df_cell['cellID'].str.replace('([ATCG]+-1)','GSM4231663').str.replace('([ATCG]+-2)','GSM4231664').str.replace('([ATCG]+-3)','GSM4231665')

df_cell[['cellID','sampleID']].drop_duplicates(['sampleID'])##查看处理效果，确实实现了cellID和sampleID的替换
###这里只处理了part_1，所以只有三个sampleID，part_2道理相同；如果替换的数目过多也可以使用循环的方式处理，这里就看个人能力了；

###有时候cellID上会直接带有sampleID（形如：GSM4231663_AAACCTGGTAAAGTCA-1）,可以用如下代码处理：
#df_cell['sampleID'] = df_cell['cellID'].str.replace('(_.+)','')#即将下划线开始的所有字符去除

如上处理完毕后我们便可根据sampleID直接将df_cell和sample进行合并，代码如下：

df_cell_tmp = pd.merge(df_cell,sample, how = 'left',left_on=['sampleID'], right_on=['geo_accession'])  
df_cell_tmp.head()###df_cell和sample合并，创建df_cell_tmp。
#df_cell_tmp中包含了df_cell和sample中所有的列，其cellID顺序和df_cell中的一致。

#这时我们只要将df_cell_tmp中的条目直接填入df_cell中即可，如下：
df_cell['meta_sourceName'] = df_cell_tmp['source_name_ch1'].tolist()
df_cell['meta_tissue'] = df_cell_tmp['characteristics_ch1.2.tissue'].tolist()
df_cell['meta_genotype'] = df_cell_tmp['characteristics_ch1.1.genotype'].tolist()

df_cell['clusterID'] = 'notAvailable'
df_cell['clusterName'] = 'notAvailable'
df_cell['sampleID'] = df_cell_tmp['geo_accession'].tolist()
df_cell['cellOntologyName'] = 'pro-T cell'
df_cell['cellOntologyID'] = 'CL:0000827'
df_cell['FACSMarker'] = df_cell_tmp['characteristics_ch1.3.facs purification gate'].tolist()

三、快速填写marker基因

情况1： markerGene文件以一个表格的形式存放

如下图：

cellMarker = pd.read_csv('../downloaded_data/GSE_supplementary_data/markerGene_C1.txt',sep='\t')
cellMarker.head()

markerGenes = {}

cellMarker = pd.read_csv('../downloaded_data/GSE_supplementary_data/GSE111360_RAW/cellMarker.txt',sep='\t')
cellMarker.head()

for name,group in cellMarker.groupby('cluster'):
    print(name)
    print(len(group['gene'].tolist()))
    markerGenes[name] = {}#cluster1为clusterName，每个clusterName都要有以下五个keywords
	# 下面等号右边填为List格式，且长度一致，没有的填notAvailable，也需为list格式
    markerGenes[name]['geneSymbol'] = group['gene'].tolist()
    markerGenes[name]['ensemblID'] = my_builder.calculate_ensemblID(group['gene'].tolist())
    markerGenes[name]['pValue'] = group['p_val'].tolist()
    markerGenes[name]['statistics'] = group['avg_logFC'].tolist()# 统计量
    markerGenes[name]['statisticsType'] = ['wilcoxon' for i in range(len(group['gene'].tolist()))]

情况2： markerGene文件以多重Excel表格的形式存放,不同的sheet用对应的clusterName表示

如下图：

##分开的文件保存marker基因	
clusters = set(df_cell['clusterID'])
markerGenes = {}
for i in clusters:
    df_marker = pd.read_excel('../downloaded_data/GSE_supplementary_data/GSE118257_markerGenes.xlsx', sheet_name = i)
    markerGenes[i] = {}
    markerGenes[i]['geneSymbol'] = df_marker['gene'].tolist()
    markerGenes[i]['pValue'] = df_marker['p_val'].tolist()
    markerGenes[i]['statistics'] = df_marker['avg_logFC'].tolist()
    markerGenes[i]['ensemblID'] = my_builder.calculate_ensemblID(markerGenes[i]['geneSymbol'])
    markerGenes[i]['statisticsType'] = ['FScore']*len(markerGenes[i]['geneSymbol'])
    print(i)

补充说明：

绝大多数时候都是情况1，当碰见其他情况是，可以考虑手动将它们处理成上述两种情况再进行处理，速度会快很多；
两种情况都要保证cluster列和cellAnnotation中的clusterName相匹配（除非cellAnnotation中的cluster是用函数自动生成的）

Author: Aimin Xie
UserID: user_34, 有任何疑问可以直接问我；

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

数据处理代码-新手教程.md

数据处理代码-新手教程.md

一、矩阵处理

情况一：常规矩阵

例子1：常规处理

补充例子1：最常用的转换获取TPM的操作

情况二：稀疏矩阵

情况2.1: 单个稀疏矩阵

例子1: 数据不需要做任何操作，直接读取然后保存

例子2: 需要提取部分细胞（直接对稀疏矩阵进行细胞筛选）

情况2.2: 多个稀疏矩阵直接进行拼接再保存

补充例子2：常规矩阵和稀疏矩阵的转换

补充例子3：如何处理HDF5(.h5)格式文件

.h5ad文件处理进阶版（更简单）

二、cellAnnotion处理

例子1: cellID刚好能够和title对应

例子2: 存在某个线索可以实现cellID和sampleID的转化

三、快速填写marker基因

情况1： markerGene文件以一个表格的形式存放

情况2： markerGene文件以多重Excel表格的形式存放,不同的sheet用对应的clusterName表示

Files

数据处理代码-新手教程.md

Latest commit

History

数据处理代码-新手教程.md

File metadata and controls

一、矩阵处理

情况一：常规矩阵

例子1： 常规处理

补充例子1：最常用的转换获取TPM的操作

情况二：稀疏矩阵

情况2.1: 单个稀疏矩阵

例子1: 数据不需要做任何操作，直接读取然后保存

例子2: 需要提取部分细胞（直接对稀疏矩阵进行细胞筛选）

情况2.2: 多个稀疏矩阵直接进行拼接再保存

补充例子2：常规矩阵和稀疏矩阵的转换

补充例子3：如何处理HDF5(.h5)格式文件

.h5ad文件处理进阶版（更简单）

二、cellAnnotion处理

例子1: cellID刚好能够和title对应

例子2: 存在某个线索可以实现cellID和sampleID的转化

三、快速填写marker基因

情况1： markerGene文件以一个表格的形式存放

情况2： markerGene文件以多重Excel表格的形式存放,不同的sheet用对应的clusterName表示

例子1：常规处理