1. Seurat对象查看当前的Assay

pbmc <- Read10X('./filtered_gene_bc_matrices/hg19/')
pbmc <- CreateSeuratObject(pbmc,project = 'pbmc3k',min.cells = 3,min.features = 200)
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc) %>% FindVariableFeatures(nfeatures = 2000) %>% ScaleData(vars.to.regress = "percent.mt")
pbmc <- SCTransform(pbmc, vars.to.regress = "percent.mt")
DefaultAssay(pbmc)
#[1] "SCT"

在进行了SCTransform操作后，矩阵默认会变成SCT矩阵，如果不加设置，后续的PCA等操作都是基于SCT矩阵。
修改DefaultAssay：

DefaultAssay(pbmc) <- 'RNA'
DefaultAssay(pbmc) 
#[1] "RNA"

2. Seurat使用FindVariables找到高变基因后增删高变基因

var.gene <- pbmc@assays$RNA@var.features #提取高变基因
pbmc@assays$RNA@var.features <- c(var.gene,c('CD4','CD8'))%>%unique() #增加高变基因
#unique()是为了保证增加的基因不和原有的重复

3. 不同运行步骤中的文件来源和储存位置⚠️

• PCA输入的是@scale.data❗️
  因为做PCA的前提是：1. 主成分分析认为主元之间彼此正交，样本呈高斯分布；2. 主成分分析假设源信号间彼此非正交。
  NormalizeData得到的@data矩阵消除了测序文库差异（对于每个细胞，将每个基因的count除以总数，然后乘以一个scale.factor, 之后以自然对数进行转换），但得到的矩阵仍然是非高斯分布的。而ScaleData对矩阵进行了中心化，得到的@scale.data矩阵结果接近于高斯分布。

• RunTSNE和RunUMAP输入的都是PCA的数据（PCA@cell.embedding），输出结果保存在各自的数据槽 。

• FindNeighbors输入的也是PCA的数据（根据PCA降维结果判断哪个细胞和哪个系细胞距离更近）。这个函数构建SNN矩阵，结果保存在pbmc@graph下面（输入的是SCT的PCA得到的是SCT_nn和SCT_snn，输入RNA的PCA得到RNA_nn和RNA_snn）。

• Harmony读入的是PCA（要注意是RNA的pca还是SCT的pca）的数据，进行调整，在和pca一起保存在reduction下面。其基因数和PC数都和PCA一样。

⚠️：PCA的值是可以被覆盖的，使用三步法对矩阵进行标准化后进行PCA后再使用SCT矩阵进行标准化，PCA的矩阵变成了SCT的PCA矩阵，原有的PCA矩阵不会保留。后续的TSNE和UMAP降维图也和三步法不一样。

• FindClusters输入的是FindNeighbors的结果，其运行结果保存在metadata，使用不同的resolution运行时，每运行一次，meta.data中就会多出一列，记录不同的resolution的分群结果。meta.data中的seurat_cluster记录的是最后一次运行聚类的结果。

• 细胞周期评分用的是@data的数据❗️（经log转换后的矩阵Normalize Data），运行的结果储存在meta.data中

• FindAllMarkers和FindMarkers输入的是@data的数据❗️

• AverageExpression默认输入的也是@data的数据❗️，但可以通过slot参数来改，比如设置slot=‘scale.data’
  延伸阅读：https://mp.weixin.qq.com/s/s25DLc-tj0lPAcsPurn89Q

• 差异分析和数据整合使用的也是@data的数据❗️，虽然官方更推荐使用@scale.data的数据，但目前Seurat还不支持。

• FindVariableFeatures()输入的是@counts的数据（❗️不是Normalized-Data也不是Scaled_Data）
  参考：https://www.biostars.org/p/406388/

FindVariableFeatures()就是在细胞与细胞间进行比较，选择表达量差别最大的基因，默认情况下，会返回2,000个高可变基因。高可变基因的数目用于做PCA。
Q：PCA输入的是scale.data的数据，scale.data默认做的是FindVariableFeatures()找到的高变基因的scale，但是scale.data如果scale了全部基因，做PCA也就是用的全部基因，那FindVariableFeatures()还有什么用？
A：虽然很多教程都提示内存足够的话建议scale所有的基因，但PCA的时候还是要选择高变基因，否则会引入噪声。低丰度，变化低的基因，很多都是噪音，如果用全部基因去跑PCA，这种低丰度基因的噪音是没有办法通过选择pc数去除的。仔细看Seurat官网，虽然它推荐scale全部基因，但是做PCA的时候还是限定了使用高变基因。

https://satijalab.org/seurat/articles/pbmc3k_tutorial.html

此外如果使用了SCTransform，默认是得到3000个高变基因，Scale的也是高变基因。再下一步直接做RunPCA()，也是默认使用高变基因。

4. @data标准化矩阵和@scale.data 归一化矩阵的区别
单细胞RNA 测序数据中，文库之间测序覆盖率的系统差异通常是由细胞间的cDNA 捕获或PCR 扩增效率方面的技术差异引起的，这归因于用最少的起始材料难以实现一致的文库制备。标准化旨在消除这些差异（例如长度，GC 含量），以使它们不干扰细胞之间表达谱的比较。这样可以确保在细胞群体中观察到的任何异质性或差异表达都是由生物学而不是技术偏倚引起的（批次矫正仅在批次之间发生，并且必须同时考虑技术偏差和生物学差异，标准化只需考虑技术差异）。
软件Seurat 提供了三种标准化的方法，分别为LogNormalize、CLR、RC，通常情况下我们采用LogNormalize 的方式进行标准化，计算公式为：log1p((Feature counts/total counts) ∗ scale. factor)

归一化的目的则是使特征具有相同的度量尺度

参考：Seurat的normalization和scaling

5. 关于有些细胞属于同一个cluster但是在umap或者tsne图上相聚较远的问题：UMAP和TSNE是各自的算法在PCA降维的基础上再进行非线形降维，在二维图上把其各自算法认为相近的细胞聚在一起。但是FindClusters输入的不是UMAP或TSNE降维的数据，而是FindNeighbors的数据，而FindNeighbors输入的数据是PCA降维数据，是用另外一种算法计算的细胞之间的距离。因此会出现有些细胞被认为是同一个cluster，但是在umap或者tsne图上相聚较远（尤其是一些散在的，脱离主群的细胞）

6. marker基因鉴定，查看marker基因的表达是使用RNA矩阵还是sct矩阵？
这是一个争议性问题，两个都可以，目前建议最好使用RNA矩阵。
sct的到的count并不是真实的基因表达值，而是通过scaledata倒推出来的，它是一个回归，运算之后的残差。

7. 关于FindAllMarks找到的基因
如下图，先看cluster0的Marker基因：cluster0的差异基因是cluster0的细胞和剩下的所有的cluster合在一起的细胞做对比得到的。pct.1是这个基因在cluster0中的表达比例（S100A8在cluster0的细胞中的表达比例是100%），pct.2是这个基因在除了cluster0以外的所有细胞中的表达比例（S100A8在除cluster0以外的细胞中的表达比例是51.2%）。avg_log2FC是表达差异倍数，p_val_adj是校正后的p值。

8. 在提取Marker基因时比较好的办法：因为单细胞矩阵算出来的结果，p_val_adj==0的有很多，所以可以先把p_val_adj==0先提出来。再把p_val_adj<0.01的按差异倍数靠前的20/30/500...（按需要）个基因提出来，然后把两个矩阵合在一起（取交集）用来做细胞鉴定。（结合p值和fc来做筛选效果更好）

top = 30   #可根据需要调整
TopMarkers1 <- ClusterMarker %>% filter(p_val_adj == 0) %>% group_by(cluster) %>% 
  top_n(n = top, wt = avg_log2FC)
TopMarkers2 <- ClusterMarker %>% filter(p_val_adj < 0.01) %>% group_by(cluster) %>%
  top_n(n = top, wt = avg_log2FC)
TopMarkers <- rbind(TopMarkers1, TopMarkers2) %>% unique.matrix() %>% arrange(cluster)
write.csv(TopMarkers,'CellType/TopMarkers.csv', row.names=F)

⚠️：提取没有核糖体和线粒体的marker基因更好。（这些基因对鉴定没有帮助）

人的核糖体基因是RPS/RPL开头
人的线粒体基因是MT-开头
小鼠的核糖体基因是Rps/Rpl开头
小鼠的线粒体基因是mt-开头

##人
### 提取没有线粒体、核糖体的Markers
ClusterMarker_noRibo <- ClusterMarker[!grepl("^RP[SL]", ClusterMarker$gene, ignore.case = F),]
ClusterMarker_noRibo_noMito <- ClusterMarker_noRibo[!grepl("^MT-", ClusterMarker_noRibo$gene, ignore.case = F),]
top = 100   #可根据需要调整
TopMarkers1 <- ClusterMarker_noRibo_noMito %>% filter(p_val_adj == 0) %>% 
  group_by(cluster) %>% top_n(n = top, wt = avg_log2FC)
TopMarkers2 <- ClusterMarker_noRibo_noMito %>% filter(p_val_adj < 0.01) %>% 
  group_by(cluster) %>% top_n(n = top, wt = avg_log2FC)
ClusterMarker_noRibo_noMito <- rbind(TopMarkers1, TopMarkers2) %>% unique.matrix() %>% arrange(cluster)
write.csv(ClusterMarker_noRibo_noMito,'CellType/TopMarkers_noRibo_noMito.csv', row.names=F)

##小鼠
### 提取没有线粒体、核糖体的Markers
ClusterMarker_noRibo <- ClusterMarker[!grepl("^Rp[sl]", ClusterMarker$gene, ignore.case = F),]
ClusterMarker_noRibo_noMito <- ClusterMarker_noRibo[!grepl("^mt-", ClusterMarker_noRibo$gene, ignore.case = F),]

有些基因比如Foxp3，对细胞鉴定很重要，但常常在筛选Marker基因的时候筛选不出来。非负矩阵分解可能更好。
参考：过滤线粒体核糖体基因

9. 提取亚群

Idents(pbmc) <- "celltype" #先把celltype设成Idents（active ident)
CD4T <- subset(pbmc, idents = "CD4_T") #通过idents来提取
tmp <- c("orig.ident", "percent.mt","percent.rb","percent.HB","S.Score",
         "G2M.Score","Phase","DF.classifications")
CD4T <- CreateSeuratObject(CD4T@assays$RNA@counts, meta.data = CD4T@meta.data[,tmp])
#metadata保留了，降维聚类结果没有保留

⚠️新提取的亚群需要重新进行降维聚类（和大群相比，标记基因发生了变化），并重新寻找marker基因，重新分群，注释。❗️subset提取子集后，不同样本间批次校正的信息也被去除，需要重新进行批次矫正

参考：Seurat取子集时会用到的函数和方法⚠️⚠️⚠️

10. 取子集后如何去除空子集（还存在这个level，但里面包含的细胞为0，如何去除）as.factor(as.character())

Idents(CD4T) <- "seurat_clusters"
CD4T.sub <- subset(CD4T, idents = c(0,1,2,3,4,5,6))
CD4T.sub$seurat_clusters <- as.factor(as.character(CD4T.sub$seurat_clusters))
CD4T.sub$celltype <- as.factor(as.character(CD4T.sub$celltype))
saveRDS(CD4T.sub, "CD4T/CD4T.classified.rds")

11. 双细胞的预测和去除如DoubletFinder建议单样本进行，不建议双样本一起预测。除此之外，其他步骤都可以多样本一起做，质控的时候也可以多样本一起做，但是建议每个样本都单独看一看。

12. 单细胞多样本整合：merge()；多样本拆分：SplitObject()

scRNA <- merge(scRNAlist[[1]], scRNAlist[2:length(scRNAlist)])
scRNAlist <- SplitObject(scRNA, split.by = "orig.ident")

13. 在做多组数据整合，每个组又有多个样本的时候，最好把单独的每个样本当成批次，而不是把不同的组当成批次。

14. 多核运算

library(future)
options(future.globals.maxSize = 20 * 1024^3)
plan(multisession, workers = 8) #开启多核运算 (8个核)

plan('sequential') #终止多核运算

参考：单细胞数据分析中future包的使用

15. pbmc3k.final@commands$FindClusters 可以查看FindClusters运行时间和记录。Seurat是记录其分析过程的，也可查看command下其他操作

16. 关于质控标准：同一组织的最好用同样的标准，不同组织的可以不一样。不同组织线粒体含量等可能存在差异。

17. 可视化的方法总结
参考：//www.greatytc.com/p/0d1e2c7d21a4

18. circos图绘制

19. 单细胞数据思维导图，有利于查看单细胞数据格式。
//www.greatytc.com/p/7560f4fd0d77

20. 对于旧版本Seurat对象的更新
scRNA <- UpdateSeuratObject(scRNA)
UpdateSeuratObject {SeuratObject} ：Update old Seurat object to accommodate new features

21. 对有些操作需要用到python设置的情况

Sys.setenv(RETICULATE_PYTHON="/usr/bin/python3")
# ⚠️要根据自己python3的路径来设置，可以在终端使用which python3来查看 

22. 单细胞数据做pooling的好处：可以尽量的降低dropout的问题。（dropout就是矩阵中的zero，这些zero实际上并不是0，而是每个液滴里面起始反应量太低了。而一般的反转录效率只能到30%左右，70%的转录本实际上在反转录那一步是被丢掉的，这是单细胞测序一个比较大的问题）。
但是一旦做了pooling，你必须要证明pooling对结果是没有影响的（下图的右面三个图）。

这个文章用了pooling

23. Seurat的VlnPlot中的combine参数，在如下画三个基因的情况下，设置成T就画一张图，设置成False，会将三个基因各画一张图。

VlnPlot(seuratObj, features = c("Il15", "Cxcl10","Cxcl16"), split.by = "aggregate", pt.size = 0, combine = T)

24. rev()这一步是将横坐标的基因反过来排序

DotPlot(scRNA, features = nichenet_output$top_ligands, cols = "RdYlBu") + RotatedAxis()
DotPlot(seuratObj, features = nichenet_output$top_ligands %>% rev(), cols = "RdYlBu") + RotatedAxis()

这两个画出来的图横坐标基因的顺序是相反的（见NicheNet）

25. 堆叠小提琴图的绘制
完成这个需求有以下几种实现方法：
1. Seurat包直接就可以实现（stack = T）
2. 通过Seuart->scanpy来实现，第一张是Seurat包VlnPlot函数画的图，第二张是scanpy中stacked_violin函数画的图，那么现在问题就变成为Seurat对象到scanpy对象的转换
3. 用R原生函数实现StackedVlnPlot的方法
4. 使用基于scanpy包衍生的scanyuan包来说实现
5. 使用R包MySeuratWrappers来实现

最简单的方法1如下：

pbmc = readRDS('pbmc.rds')
markers <- c('CD14','CD68','CD4','CD19','MS4A1','CD79A','GNLY','NKG7','GZMB','CD8A','CD8B','FCGR3A',
             'FOXP3','PPBP','FCER1A','CD34')
markers <- CaseMatch(markers, rownames(pbmc))
markers <- as.character(markers)
VlnPlot(pbmc, features = markers, pt.size = 0, group.by = 'cell_type',stack = T)+NoLegend()

26. 当绘图需要设置横坐标顺序时，先把Ident设置为需要绘图的变量，再使用factor进行设置。

Idents(scRNA) <- 'orig.ident'
My_levels <- c( "con1","con2", "con3", "case1","case2", "case3")
Idents(scRNA) <- factor(Idents(scRNA), levels= My_levels)

如果不设置level，会按字母顺序排列，case会自动排在con前面。

27. 在绘制多个质控参数的小提琴图时，可以先生成一个空list，将每张图都存成list的一个对象，再使用patchwork包的wrap_plots函数画成一张图。

plot.featrures = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.HB")
plots = list()
for(i in seq_along(plot.featrures)){
    plots[[i]] = VlnPlot(scRNA, pt.size = 0,
                         features = plot.featrures[i],log = T) + theme.set2 + NoLegend()}
wrap_plots(plots = plots, nrow=2) 
ggsave(filename = 'beforeQC9log.pdf', width = 8, height = 6)

28. 在进行细胞注释时，可以使用recode函数

scRNA$celltype <- scRNA$SCT_snn_res.0.3
scRNA$celltype <- recode(scRNA$celltype,
                         "0" = "CD14_Monocyte",
                         "1" = "CD14_Monocyte",
                         "2" = "NK_cell",
                         "3" = "B_cell",
                         "4" = "T_cell",
                         "5" = "T_cell",
                         "6" = "CD16_Monocyte",
                         "7" = "Megakaryocyte",
                         "8" = "NK_cell",
                         "9" = "T_cell",
                         "10" = "T_cell",
                         "11" = "T_cell",
                         "12" = "T_cell",
                         "13" = "T_cell",
                         "14" = "B_cell",
                         "15" = "DC",
                         "16" = "Neutrophil",
                         "17" = "T_cell",
                         "18" = "Proliferating_lymphocyte",
                         "19" = "CD14_Monocyte",
                         "20" = "B_cell")

29. 对某一metadata进行重命名

#查看meta.data中有哪些内容
colnames(scRNA@meta.data) 
#meta.data中第13个进行重命名
colnames(scRNA@meta.data)[13] <- 'Cytokine_Score' 

使用Seurat的RenameIdents函数也可以

scRNA <- RenameIdents(scRNA, 'oldIdent'='newIdent')

30. 进行marker基因绘图时，使用CaseMatch函数去除不存在的基因。

markers <- c('CD79A','CD79B','MS4A1','CD19','MZB1','GNLY','NKG7','NCAM1','CD3D','CD3E')
markers <- CaseMatch(markers, rownames(scRNA))
markers <- as.character(markers)

31. prop.table函数，生成细胞比例表格，用于绘制细胞比例条形图

prop.table(x, margin = NULL)

x: table
margin: a vector giving the margins to split by. E.g., for a matrix 1 indicates rows, 2 indicates columns, c(1, 2) indicates rows and columns. When x has named dimnames, it can be a character vector selecting dimension names.

m <- matrix(1:4, 2)
m
#      [,1] [,2]
# [1,]    1    3
# [2,]    2    4
proportions(m, 1)
#           [,1]      [,2]
# [1,] 0.2500000 0.7500000
# [2,] 0.3333333 0.6666667

#根据orig.ident来计算细胞百分比
cell.prop<-as.data.frame(prop.table(table(Idents(scRNA), scRNA$orig.ident)))
colnames(cell.prop)<-c("cluster","Group","proportion")
write.csv(cell.prop,row.names = FALSE,file = 'cell_prop_all.csv')
cell.prop$Group <- factor(cell.prop$Group,ordered=TRUE,levels=c("HC1","HC2","HC3","Case1","Case2","Case3"))
ggplot(cell.prop,aes(Group,proportion,fill=cluster))+
geom_bar(stat="identity",position="fill")+
  ggtitle("")+
  theme_bw()+
  theme(axis.ticks.length=unit(0.5,'cm'))+
  guides(fill=guide_legend(title=NULL))
ggsave(
  'NK/satck_prop_NK.pdf',
  plot = last_plot(),
  device = NULL,
  path = NULL,
  scale = 1,
  width = 10,
  height = 8,
  dpi = 300,
  limitsize = TRUE,
)

32. 在从基因表达矩阵（行为基因，列为样本）直接生成Seurat对象时，如下：

HC_1 <- CreateSeuratObject(counts = HC_1)

得到的HC_1样本的orig.ident默认是样本名中第一个_号的前一部分。所以要保证矩阵的列名是样本名_细胞barcode这样的格式。
如果有多个分组，例如两个样本矩阵中细胞分别命名为HC_1_barcode，HC_2_barcode，在直接通过如下方法得到两个Seurat对象，再对其进行merge之后，两个样本会被合并成一个。也就是样本信息只保留了第一个_号之前的HC，没有保留_号之后的1和2。

HC_1 <- CreateSeuratObject(counts = HC_1)
HC_2 <- CreateSeuratObject(counts = HC_2)
scRNA <- merge(HC_1,HC_2)
table(scRNA$orig.ident)
# HC
#3000

为了避免这种情况，可以在构建Seurat对象时通过参数进行设置

HC_1 <- CreateSeuratObject(counts = HC_1,names.field = 1:2)
HC_2 <- CreateSeuratObject(counts = HC_2,names.field = 1:2)
scRNA <- merge(HC_1,HC_2)
table(scRNA$orig.ident)
# HC_1  HC_2 
# 1500  1500

33. 在做umap和tsne的时候有很多聚不在一起的小群，可以考虑下面几种情况
    (1) 双细胞的问题
    双细胞倾向于独立成群，或位于大群的边缘或线状连接部分。
    (2) 高变基因数目选择问题
    FindVariableFeatures()时nfeatures选的越多，图越分散
    (3) 降维时pc数目的选择
    做tsne和umap降维时选择的pc数越多，图越分散。
    参考：如何确定细胞聚类的PC数

⚠️PC数的选择：Seurat官网提供的三种方法只能给出PC数的粗略范围，选择不同PC数目，细胞聚类效果差别较大，因此，需要一个更具体的PC数目。作者提出一个确定PC阈值的三个标准：

# Determine percent of variation associated with each PC
pct <- pbmc [["pca"]]@stdev / sum( pbmc [["pca"]]@stdev) * 100
# Calculate cumulative percents for each PC
cumu <- cumsum(pct)
# Determine which PC exhibits cumulative percent greater than 90% and % variation associated with the PC as less than 5
co1 <- which(cumu > 90 & pct < 5)[1]
co1
# Determine the difference between variation of PC and subsequent PC
co2 <- sort(which((pct[1:length(pct) - 1] - pct[2:length(pct)]) > 0.1), decreasing = T)[1] + 1
# last point where change of % of variation is more than 0.1%.
co2
# Minimum of the two calculation
pcs <- min(co1, co2)
pcs
# Create a dataframe with values
plot_df <- data.frame(pct = pct,   cumu = cumu,   rank = 1:length(pct))
# Elbow plot to visualize 
ggplot(plot_df, aes(cumu, pct, label = rank, color = rank > pcs)) + 
  geom_text() + 
  geom_vline(xintercept = 90, color = "grey") + 
  geom_hline(yintercept = min(pct[pct > 5]), color = "grey") +
  theme_bw()

34. 亚群注释注意事项

一般先选默认分辨率(0.8)，大概可能会分出十几个群。因为最终都是要注释到每一个barcode，所以首先可以看大类marker的分布(不受分辨率影响)，可以根据marker基因的分布来调整分辨率。是否需要精细的分群得看精细的分群对研究有没有决定作用，还有很重要的一点是看分出的各个cluster在Findallmarkers给出的结果中marker的热图是不是能明显分开。精细划分的细胞本来就很类似，如果有部分小群的热图明显分不开或者非常类似，就可以考虑把分辨率调小。

比如像这种情况，3和4可能就是同一类细胞

35. 跑umap/tsne跟跑FindNeighbors输入的都是PCA的矩阵，但是用的pc数可以是不一样的。
    在一些教程里面，我们常常看到这样的代码，在看了一下ElbowPlot之后选定了PC数就直接RunTSNE，RunUMAP和FindNeighbors都用一样的pc数。

pbmc <- RunPCA(pbmc, verbose = F)
ElbowPlot(pbmc, ndims = 50)
pc.num=1:20
pbmc <- pbmc %>% RunTSNE(dims=pc.num) %>% RunUMAP(dims=pc.num) %>%
  FindNeighbors(dims = pc.num) %>% FindClusters()

这实际上是没有必要的必须保持一致的。下游的都是用pca之后的，pca是为了压缩数据。
umap和tsne是为了可视化(仅仅是可视化)，但是FindNeighbor是计算细胞间距离矩阵。找类群数目和可视化可以说没有关系。

36. UMAP图的分群可视化

library(dplyr)
library(Seurat)
library(patchwork) 
library(tidyverse)
library(cowplot)
library(clustree)

# Load the PBMC dataset
pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
# Initialize the Seurat object with the raw (non-normalized data).
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc,dims=1:20,resolution = seq(from=0,by=.2,length=10))
pbmc <- RunUMAP(pbmc, dims = 1:10)
clustree(pbmc)

Idents(pbmc) <- "RNA_snn_res.1"
p<- map(c(levels(Idents(pbmc))),function(x){DimPlot(pbmc, cells.highlight = CellsByIdentities(object = pbmc, idents = x))})
wrap_plots(plots = p, nrow=3) 
## 对比27，不必先设置p=list()就可以直接使用wrap_plots是因为map生成的对象就是个list
# 或 plot_grid(plotlist = p)

map函数：
R语言循环第三境界：purrr包map函数！
浅析R语言中map(映射)与reduce(规约)

37. RunUMAP()不同参数对分群的影响

38. 批量读入矩阵并创建seurat对象
    ⚠️assign()的用法

# Data loading and QC

### 读入sample消息
samples <- read_excel("../data/metadata/patients_metadata.xlsx", range = cell_cols("A:A")) %>% .$sample_id
samples
#  [1] "p018t" "p018n" "p019t" "p019n" "p023t" "p024t" "p027t" "p027n" "p028n"
# [10] "p029n" "p030t" "p030n" "p031t" "p031n" "p032t" "p032n" "p033t" "p033n"
# [19] "p034t" "p034n"

### import cellranger files from different data sets
for (i in seq_along(samples)){
  assign(paste0("scs_data", i), Read10X(data.dir = paste0("../data/cellranger/", samples[i], "/filtered_feature_bc_matrix")))
}
# 读入的每一个文件都是一个对象，命名为scs_data1-20

### create seurat objects from cellranger files
for (i in seq_along(samples)){
  assign(paste0("seu_obj", i), CreateSeuratObject(counts = eval(parse(text = paste0("scs_data", i))), project = samples[i], min.cells = 3))
}
# 对每一个scs_data创建一个Seurat对象，命名为seu_obj1-20

### merge data sets
seu_obj <- merge(seu_obj1, y = c(seu_obj2, seu_obj3, seu_obj4, seu_obj5, seu_obj6, seu_obj7, seu_obj8, seu_obj9, seu_obj10, seu_obj11, seu_obj12, seu_obj13, seu_obj14, seu_obj15, seu_obj16, seu_obj17, seu_obj18, seu_obj19, seu_obj20), add.cell.ids = samples, project = "lung")

### calculate mitochondrial, hemoglobin and ribosomal gene counts
seu_obj <- PercentageFeatureSet(seu_obj, pattern = "^MT-", col.name = "pMT")
seu_obj <- PercentageFeatureSet(seu_obj, pattern = "^HBA|^HBB", col.name = "pHB")
seu_obj <- PercentageFeatureSet(seu_obj, pattern = "^RPS|^RPL", col.name = "pRP")

qcparams <- c("nFeature_RNA", "nCount_RNA", "pMT", "pHB", "pRP")
for (i in seq_along(qcparams)){
  print(VlnPlot(object = seu_obj, features = qcparams[i], group.by = "orig.ident", pt.size = 0))
}
for (i in seq_along(qcparams)){
  print(RidgePlot(object = seu_obj, features = qcparams[i], group.by = "orig.ident"))
}

VlnPlot(seu_obj, features = c("nFeature_RNA", "nCount_RNA", "pMT"), pt.size = 0, group.by = "orig.ident", ncol = 1, log = T)
ggsave2("SuppFig1B.pdf", path = "../results", width = 30, height = 20, units = "cm")

39. 绘图时颜色的传参

ClusterName_color_panel <- c(
  "Naive CD4 T" = "#DC143C", "Memory CD4 T" = "#0000FF", "CD14+ Mono" = "#20B2AA",
  "B" = "#FFA500", "CD8 T" = "#9370DB", "FCGR3A+ Mono" = "#98FB98",
  "NK" = "#F08080", "DC" = "#0000FF", "Platelet" = "#20B2AA"
)
ggplot(df, aes(Pseudotime, colour = cell_type, fill=cell_type)) +
  geom_density(bw=0.5,size=1,alpha = 0.5)+theme_classic2()+ scale_fill_manual(name = "", values = ClusterName_color_panel)+scale_color_manual(name = "", values = ClusterName_color_panel)

参考：monocle2

40. 绘图气泡图/小提琴图的时候翻转纵坐标
    使用scale_y_discrete+rev

p1=DotPlot(pbmc,features = c('CCR2','CD3D'),group.by = 'cell_type')+theme_bw()
p2=DotPlot(pbmc,features = c('CCR2','CD3D'),group.by = 'cell_type')+theme_bw()+scale_y_discrete(limits = rev(levels(pbmc$cell_type)))
p1|p2

41. 添加matadata
    custom是n行1列的矩阵，行名是细胞名，一列是想要添加的信息，比如手动注释结果等。

scRNA <- AddMetaData(scRNA, metadata=custom)

42. 批量绘制细胞群marker基因
    空转也是一样，把FeaturePlot()换成SpatialFeaturePlot()即可

markerlist <- list(
  Tcell = c("CD3D","CD3E"),
  CD4T = c("CD4","CD40LG"),
  CD8T = c("CD8A","CD8B"),
  Treg = c("FOXP3", "IL2RA"),
  Bcell = c("CD79B","MS4A1"),
  Plasma = c("MZB1", "XBP1"),
  Myeloid = c("LYZ","CST3"),
  Monocyte = c("FCN1","S100A9"),
  Macrophage = c("CD163","CD68"),
  DC_CD1 = c("CLEC9A", "BATF3"),
  DC_CD2 = c("CD1C", "FCER1A"),
  DC_LAMP3 = c("LAMP3", "FSCN1"),
  pDC = c("LILRA4", "IL3RA"),
  Neutrophil = c("FCGR3B", "CSF3R"),
  Epithelial = c("EPCAM", "KRT5"),
  Endothelial = c("PECAM1", "VWF"),
  Fibroblasts = c("ACTA2", "COL1A1")
)
for(i in names(markerlist)){
  markers <- markerlist[[i]]
    p <- FeaturePlot(pbmc, features = markers) + plot_layout()&theme(legend.position = "right")
  ggsave(paste0(i, ".pdf"), p, width = 10, height = 8)
}
markers <- do.call("c",markerlist)
p <- FeaturePlot(pbmc, features = markers,ncol = 4) + plot_layout()&theme(legend.position = "right")
ggsave("Markers_PBMC_all.pdf", p, width = 18, height = 30)

43. 查看不同cluster的中位基因

a<-pbmc@meta.data%>%group_by(seurat_clusters)%>%
  summarise(count=n(),nCount_RNA_M =median(nCount_RNA,na.rm = TRUE),nFeature_RNA_M=median(nFeature_RNA,rm=TRUE))
a

查看不同细胞群的中位基因也是一样

a=pbmc@meta.data%>%group_by(cell_type)%>%
  summarise(count=n(),nCount_RNA_M =median(nCount_RNA,na.rm = TRUE),nFeature_RNA_M=median(nFeature_RNA,rm=TRUE))
a

查看不同样品的中位基因也是一样

pbmc@meta.data%>%group_by(sample)%>%
summarise(count=n(),nCount_RNA_M =median(nCount_RNA,na.rm = TRUE),nFeature_RNA_M=median(nFeature_RNA,rm=TRUE))%>%head()

或者也可以

dim(pbmc)
df<-data.frame(pbmc$nFeature_RNA,pbmc$nCount_RNA)
apply(df,2,median)
# pbmc.nFeature_RNA   pbmc.nCount_RNA 
#               819              2213 

44. 计算线粒体基因
    直接使用^(MT|mt|Mt)-，人和小鼠的数据都适用

pbmc[['percent.mito']] <- PercentageFeatureSet(object = pbmc, pattern = '^(MT|mt|Mt)-')

45. 使用ggplot2优化Seurat绘图
    需要注意的点：
    1、坐标提取
    2、颜色对接，尤其是colorRampPalette的使用
    3、viridis的使用
    4、ggplot2的aes和aes_string

46. 使用subset去了单细胞子集之后，想要在原seurat对象上查看哪些细胞被滤掉了的方法

pbmc_new <- subset(pbmc, subset = nFeature_RNA > 500 & nFeature_RNA < 2500 & percent.mt < 5)
cellp <-rownames(pbmc_new@meta.data)
pbmc$filter <-ifelse(rownames(pbmc@meta.data) %in% cellp,"Cell","filtered")
p1 <-DimPlot(pbmc,label = T,repel = T,reduction = "umap",group.by="filter")
p1

47.提取seurat自带颜色

library(scales)
p1 <- DimPlot(scRNA)
x<-ggplot_build(p1)
info = data.frame(colour = x$data[[1]]$colour, group = x$data[[1]]$group)
info <- unique((arrange(info, group)))
cols <- as.character(info$colour)