跟学单细胞周更(一)

数据链接 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE175687

获取数据：

这里给大家介绍一下我个人下载时发现的一个小技巧

我想要下载这个数据集下四个矩阵文件

复制打包好的tar文件的http链接 发现是一个跳转链接，无法直接服务器中下载(如果点击下载下到自己电脑上，再传给服务器，有点麻烦)

检索资料，发现有人说使用curl下载 -L参数可以跟随重定向进行下载

但经curl -L -O测试使用，发现下载到的文件和前面直接wget是一样的

于是查看这个下载下来的文件

发现是一个html的格式

并且内部包含了真正的ftp下载链接

这时便可以直接下载了

可以发现这里是h5文件，使用 Read10X_h5函数读取

初探单细胞下游 可以记住这个10X技术输出文件目录和格式，以后使用 Read10X函数读取Seurat对象时可以检查一下

project = gsub('_filtered_feature_bc_matrix.h5','',gsub('^GSM[0-9]*_','',pro) )

从字符串 pro中去除以"GSM"开头、以"_"结尾的部分，以及去除字符串中的"_filtered_feature_bc_matrix.h5"部分

可以看到效果：

这对我们区分原始样本分组，了解不同实验分组的数据分布很有帮助

特别是在一开始决定需不需要去除批次效应，这一点我们在后面的推文也会谈到【flag】

初探单细胞下游

Step2-质控

过滤指标

• 最少表达基因数的细胞&最少表达细胞数的基因
• 线粒体/核糖体基因比例
• 过滤特定基因

Step3-降维聚类分群

• 标准化
• 识别高变基因
• 归一化

PCA：sce <- RunPCA(sce, features = VariableFeatures(object = sce))可以发现，并不是拿所有的基因进行PCA，而是前面识别高变基因找到的VariableFeatures，可以类比bulk分析中先看看所有基因PCA分组情况，如果效果不好可以进一步看看MAD前5000或差异表达的基因

关于harmony算法：

library(harmony)
seuratObj <- RunHarmony(sce, "orig.ident")

需要记住的是这个算法是基于我们前面的PCA的：

Harmony应用主成分分析，将转录组表达谱 嵌入到低维空间中，然后应用迭代过程去除数据集特有的影响

#设置不同的分辨率，观察分群效果(可视化后选择合适分辨率)
for (res in c(0.01, 0.05, 0.1, 0.2, 0.3, 0.5,0.8,1)) {
  sce.all=FindClusters(sce.all, #graph.name = "CCA_snn", 
                       resolution = res, algorithm = 1)
}
colnames(sce.all@meta.data)
apply(sce.all@meta.data[,grep("RNA_snn",colnames(sce.all@meta.data))],2,table)

初探单细胞下游

Step4-检查常见分群情况

Step5-细胞亚群生物学命名

这里就是需要我们手动鉴定细胞类型的地方，根据我们自己给的marker的表达分布情况来判断

初探单细胞下游

marker的选择会很大程度上影响我们判断细胞类型

原推文选用的marker：

genes_to_check = c('PTPRC', 'CD3D', 'CD3E', 'CD4','CD8A',
                   'CD19', 'CD79A', 'MS4A1' ,
                   'IGHG1', 'MZB1', 'SDC1',
                   'CD68', 'CD163', 'CD14', 
                   'TPSAB1' , 'TPSB2',  # mast cells,
                   'RCVRN','FPR1' , 'ITGAM' ,
                   'C1QA',  'C1QB',  # mac
                   'S100A9', 'S100A8', 'MMP19',# monocyte
                   'FCGR3A',
                   'LAMP3', 'IDO1','IDO2',## DC3 
                   'CD1E','CD1C', # DC2
                   'KLRB1','NCR1', # NK 
                   'FGF7','MME', 'ACTA2', ## fibo 
                   'DCN', 'LUM',  'GSN' , ## mouse PDAC fibo 
                   'MKI67' , 'TOP2A', 
                   'PECAM1', 'VWF',  ## endo 
                   'EPCAM' , 'KRT19','KRT7', # epi 
                   'FYXD2', 'TM4SF4', 'ANXA4',# cholangiocytes
                   'APOC3', 'FABP1',  'APOA1',  # hepatocytes
                   'Serpina1c',
                   'PROM1', 'ALDH1A1' )

关于marker的选择，参考，我们后面的推文也会提到【flag】

然而初步的marker鉴定只是大致分类，细胞更加细致的亚型，需要我们选取对应的细胞类型数据，进一步重新降维分群

这个时候也需要更加细分的亚型对应的marker

下面我们在本学习专辑中将第一次遇到这种情况

可视化后可以看出来在mac细胞还混入了一部分上皮细胞

进一步细分上皮细胞

Step6-提取指定的上皮细胞单细胞亚群

sce=sce[,sce$celltype=='epi']
#sce.epi=sce[,sce$RNA_snn_res.0.8 %in% c(1,3,4,7,18,22,24,25,26,32,34,35)]用这句代码同样的效果。

这里是直接获取前面分好的细胞类型对应的Seurat对象

后面我们也会介绍，重新根据样本信息重新创建Seurat对象

两样本原代乳腺癌细胞与肺部扩散癌细胞单细胞测序数据分析

需要特别注意的是，因为是重新降维分群，所以counts应该输入的是原始表达矩阵，而不是处理过的矩阵，相关标准化、识别 高变基因、归一化都需要再做

在上一期Seurat对象内部结构评论区有人提到了这个问题，这确实是需要注意的

• counts

保存的是未经处理的原始数据，适合存放稀疏矩阵

• data

原始数据经过标准化后，会存放在@data中，和counts 一样也是一个特殊的 Matrix 对象

• scale.data

当数据进行scale归一化后，存放在名为scale.data中

然后就是同样的流程了，走到选择新的细分marker进行细胞亚群生物学命名

往期回顾

单细胞多组学刻画EGCs特征

上下游，合体！

单细胞测序最好的教程（七）: 数据整合与批次效应校正

我问汤神的第一个单细胞问题，他的问答竟然是_______

单细胞分析工具||ShinyCell交互式展示单细胞数据