这个新专辑有以下几点
而我在将自己的学习笔记排版成推文时也会遵循以下行文特点:
前面的学习中我们已经掌握了基本的单细胞上下游分析流程,接下来就是两个基本方向,①加深对基础流程代码的理解,夯实基础;②在基础上拓宽加深。而在学bulk转录组分析时我也是跟着转录组周更走下来,所以接下来本专辑将会开启一部分单细胞周更的跟学,在更加细致地深入、夯实基础代码的同时进行进阶。
本片推文将学习:整合素β1基因敲除前后小鼠肺腺癌上皮细胞水平变化
数据链接 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】
过滤指标
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)
这里就是需要我们手动鉴定细胞类型的地方,根据我们自己给的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细胞还混入了一部分上皮细胞
进一步细分上皮细胞
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对象内部结构评论区有人提到了这个问题,这确实是需要注意的
保存的是未经处理的原始数据,适合存放稀疏矩阵
原始数据经过标准化后,会存放在@data中,和counts 一样也是一个特殊的 Matrix 对象
当数据进行scale归一化后,存放在名为scale.data中
然后就是同样的流程了,走到选择新的细分marker进行细胞亚群生物学命名