脑脊液单细胞转录组数据处理
前些天我们公众号弄了一个活动，详见：春节期间单细胞转录组数据分析全免费，收到了上百个需求， 本来呢我们自己就算是春节前后14天不吃不喝不眠不休也不可能完成这么多单细胞数据处理。好在我灵机一动，想起来了前面两个月培养的一百多个在线实习生，毕竟教了他们R语言，转录组，以及单细胞转录组。
所以我写了一个还算是比较自动化的单细胞转录组数据处理代码，如果是我自己的，可以在十几分钟就完成复现文章的第一层次降维聚类分群图，比如数据集：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE200164
来源于文章：Cerebrospinal fluid immune dysregulation during healthy brain aging and cognitive impairment. Cell 2022 Dec 22;185(26):5028-5039.e13. PMID: 36516855
主要是比较两个分组的CSF immune systems：
- 45 cognitively typical subjects ranging from 54-82 years old.
- 14 subjects with mild cognitive impairment or Alzheimer’s disease
虽然有59个样品但是单细胞总数才7万多，因为是比较免疫细胞，所以第一层次降维聚类分群就是跟前面一直举例的PBMC单细胞数据集类似：
免疫细胞分群
文章的主图质量很高，更值得一提的是全部的绘图代码都是公开的哦，在：Github: https://github.com/gatelabnw/csf_aging，只需要你会降维聚类分群，后续的各种可视化基本上就是复制粘贴调试代码的小活了：
Fig1D-F_Fig4H_SuppFig1C_1G_celltype_annotation.R
Fig1G-H_linear_modeling.R
Fig2A-B_SuppFig2_loess.R
Fig2C_mast_agebin.R
Fig2D_nonclassical_mono_loess.R
Fig2F_2H_de-swan.R
Fig2G_de-swan_lm_comparison.R
Fig2I_APOE_APOC1_PLTP_loess_diagnosis.R
Fig2J_mast_diagnosis_advanced.R
Fig3A-E_cellchat.R
Fig3F_Fig4E_CXCL16_CXCR6_expression.R
Fig4A_levenshtein_network.R
Fig4B_ci_hc_connections.R
Fig4C-D_mast_diagnosis_clonality.R
Fig4F_SuppFig5G_CXCR6_supervised_celltypes.R
Fig4G_SuppFig5D_supervised_celltype_annotation.R
Fig4H_cd8tem_CXCR6.R
Fig4J_olink_quanterix_partial_correlation.R
Fig4K_somamer_partial_correlation.R
Fig5A_SuppFig7A_ad_risk_gene.R
SuppFig1E_soupx_visual.R
SuppFig1F-H_general_qc.R
SuppFig3A-B_mast_age_bin_sex.R
SuppFig4A-B_mast_diagnosis.R
SuppFig4C_treg_umap.R
SuppFig5B-C_CXCR4_CXCR6.R
SuppFig5F_supervised_mast_diagnosis.R
SuppFig6A-B_olink_somamer_cxcl16_over_diagnosis.R
首先读取作者提供的表达量矩阵文件
因为作者在GEO上面给的一个表达量矩阵文件（GSE200164_counts.csv.gz），所以很简单的读取方式：
ct=fread( 'GSE200164_counts.csv.gz',data.table = F)
ct[1:4,1:10]
length(unique(ct$V1))
ct[nrow(ct),1:4]
rownames(ct)=ct[,1]
ct=ct[,7:ncol(ct)]
ct[1:4,1:4]
sce.all = CreateSeuratObject(counts = ct )
as.data.frame(sce.all@assays$RNA@counts[1:10, 1:2])
head(sce.all@meta.data, 10)
table(sce.all$orig.ident)
rm(ct)
值得注意的是虽然这个脑脊液单细胞转录组数据有五十多个样品，但是总细胞数量并不大，所以上面的代码很容易运行，而且内存要求并不高，但是如果细胞数量超过了20万，上面的代码在 CreateSeuratObject(counts = ct ) 会报错，其实有更好的办法绕过它。
然后质量控制并且降维聚类分群即可
质量控制其实每个数据集不一样的，取决于单细胞转录组来源的技术，特殊情况下需要个性化的修改 scRNA_scripts/qc.R 里面的代码：
###### step2:QC质控 ######
dir.create("./1-QC")
setwd("./1-QC")
# 如果过滤的太狠,就需要去修改这个过滤代码
source('../scRNA_scripts/qc.R')
sce.all.filt = basic_qc(sce.all)
print(dim(sce.all))
print(dim(sce.all.filt))
setwd('../')
但是多样品整合，以及后续的降维聚类分群这个代码在任意单细胞转录组数据集都是大概率不需要修改的：
sp='human'
###### step3: harmony整合多个单细胞样品 ######
dir.create("2-harmony")
getwd()
setwd("2-harmony")
source('../scRNA_scripts/harmony.R')
# 默认 ScaleData 没有添加"nCount_RNA", "nFeature_RNA"
# 默认的
sce.all.int = run_harmony(sce.all.filt)
setwd('../')
###### step4: 降维聚类分群和看标记基因库 ######
# 原则上分辨率是需要自己肉眼判断,取决于个人经验
# 为了省力,我们直接看 0.1和0.8即可
table(Idents(sce.all.int))
table(sce.all.int$seurat_clusters)
table(sce.all.int$RNA_snn_res.0.1)
table(sce.all.int$RNA_snn_res.0.8)
getwd()
dir.create('check-by-0.1')
setwd('check-by-0.1')
sel.clust = "RNA_snn_res.0.1"
sce.all.int <- SetIdent(sce.all.int, value = sel.clust)
table(sce.all.int@active.ident)
source('../scRNA_scripts/check-all-markers.R')
setwd('../')
getwd()
dir.create('check-by-0.8')
setwd('check-by-0.8')
sel.clust = "RNA_snn_res.0.8"
sce.all.int <- SetIdent(sce.all.int, value = sel.clust)
table(sce.all.int@active.ident)
source('../scRNA_scripts/check-all-markers.R')
setwd('../')
getwd()
然后就是合理的生物学命名了
我给大家准备的基因列表主要是针对肿瘤单细胞的第一层次降维聚类分群 ， 是：
- immune (CD45+,PTPRC),
- epithelial/cancer (EpCAM+,EPCAM),
- stromal (CD10+,MME,fibo or CD31+,PECAM1,endo)
参考我前面介绍过 CNS图表复现08—肿瘤单细胞数据第一次分群通用规则，这3大单细胞亚群构成了肿瘤免疫微环境的复杂。绝大部分文章都是抓住免疫细胞亚群进行细分，包括淋巴系（T,B,NK细胞）和髓系（单核，树突，巨噬，粒细胞）的两大类作为第二次细分亚群。但是也有不少文章是抓住stromal 里面的fibo 和endo进行细分，并且编造生物学故事的。
celltype=data.frame(ClusterID=0:17 ,
celltype= 0:17)
#定义细胞亚群
celltype[celltype$ClusterID %in% c( 10, 1,8),2]='CD8'
celltype[celltype$ClusterID %in% c(0,2,3,6,7,9,12),2]='CD4'
celltype[celltype$ClusterID %in% c( 15 ),2]='Bcells'
celltype[celltype$ClusterID %in% c( 11 ),2]='NK'
celltype[celltype$ClusterID %in% c( 16 ),2]='cycle'
celltype[celltype$ClusterID %in% c( 4,17 ),2]='MAC'
celltype[celltype$ClusterID %in% c( 13 ),2]='mono'
celltype[celltype$ClusterID %in% c( 5 ),2]='cDC'
celltype[celltype$ClusterID %in% c( 14 ),2]='pDC'
可以看到我们的命名还算是跟原文差不多了：
跟原文差不多