首先打开rstudio

打开后用getwd()查看当前工作路径

getwd()

如果路径跟上次的不一样，重新设置一下路径

setwd("X:/xxx/xxxx")
setwd("F:/Rstudio_data/001_singlecell_code/raw/BC21/")

重新加载R包

library(Seurat)
library(tidyverse)
library(patchwork)
library(dplyr)

加载数据

scRNA <-load("scRNA1.Rdata")

官方推荐是2000个高变基因，很多文章也有设置30000的，这个因自己的实验项目决定

scRNA1 <- FindVariableFeatures(scRNA1, selection.method = "vst", nfeatures = 3000) 

Identify the 10 most highly variable genes，把top10的高变基因挑选出来，目的是为了作图

top10 <- head(VariableFeatures(scRNA1), 10) 

plot variable features with and without labels 画出来不带标签的高变基因图

plot1 <- VariableFeaturePlot(scRNA1) 

把top10的基因加到图中

plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE, size=2.5) 
plot <- CombinePlots(plots = list(plot1, plot2),legend="bottom") 
plot 

数据标准化（中心化）

如果内存足够最好对所有基因进行中心化

scale.genes <-  rownames(scRNA1)
scRNA1 <- ScaleData(scRNA1, features = scale.genes)

如果内存不够，可以只对高变基因进行标准化

scale.genes <-  VariableFeatures(scRNA)
scRNA <- ScaleData(scRNA, features = scale.genes)

scRNA对象中原始表达矩阵经过标准化和中心化之后，已经产生了三套基因表达数据，可以通过以下命令获得

原始表达矩阵

GetAssayData(scRNA,slot="counts",assay="RNA") 

标准化之后的表达矩阵

GetAssayData(scRNA,slot="data",assay="RNA")

中心化之后的表达矩阵

GetAssayData(scRNA,slot="scale.data",assay="RNA") 

细胞周期回归：上一步找到的高变基因，常常会包含一些细胞周期相关基因。

它们会导致细胞聚类发生一定的偏移，即相同类型的细胞在聚类时会因为细胞周期的不同而分开。

cc.genes
CaseMatch(c(cc.genes$s.genes,cc.genes$g2m.genes),VariableFeatures(scRNA1))

细胞周期评分

g2m_genes = cc.genes$g2m.genes
g2m_genes = CaseMatch(search = g2m_genes, match = rownames(scRNA1))
s_genes = cc.genes$s.genes
s_genes = CaseMatch(search = s_genes, match = rownames(scRNA1))
scRNA1 <- CellCycleScoring(object=scRNA1,  g2m.features=g2m_genes,  s.features=s_genes)

查看细胞周期基因对细胞聚类的影响

scRNAa <- RunPCA(scRNA1, features = c(s_genes, g2m_genes))
p <- DimPlot(scRNAa, reduction = "pca", group.by = "Phase")
p

VlnPlot(scRNAa, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.HB","G2M.Score","S.Score"), ncol = 6)
ggsave("cellcycle_pca.png", p, width = 8, height = 6)

如果需要消除细胞周期的影响

scRNAb <- ScaleData(scRNA1, vars.to.regress = c("S.Score", "G2M.Score"), features = rownames(scRNA1))

PCA降维并提取主成分

PCA降维

scRNA1 <- RunPCA(scRNA1, features = VariableFeatures(scRNA1)) 
plot1 <- DimPlot(scRNA1, reduction = "pca", group.by="orig.ident") 
### 画图
plot1 
####确定数据的维度 Determine the ‘dimensionality’ of the dataset 
###ElbowPlot() 可以快速的检查降维的效果
plot2 <- ElbowPlot(scRNA1, ndims=20, reduction="pca") 
##画图
plot2
###我们一般选择拐点作为降维的度数。
plotc <- plot1+plot2
ggsave("pca.pdf", plot = plotc, width = 8, height = 4) 
ggsave("pca.png", plot = plotc, width = 8, height = 4)

后续分析要根据右图选择提取的pc轴数量，一般选择斜率平滑的点之前的所有pc轴，此图我的建议是选择前13个pc轴。

可以看出大概在PC为13的时候，每个轴是有区分意义的。

pc.num=1:13

细胞聚类

Identify clusters of cells by a shared nearest neighbor (SNN) modularity optimization based clustering algorithm. First calculate k-nearest neighbors and construct the SNN graph. Then optimize the modularity function to determine clusters. For a full description of the algorithms, see Waltman and van Eck (2013) The European Physical Journal B. Thanks to Nigel Delaney (evolvedmicrobe@github) for the rewrite of the Java modularity optimizer code in Rcpp!

scRNA1 <- FindNeighbors(scRNA1, dims = pc.num) 
scRNA1 <- FindClusters(scRNA1, resolution = 0.5)

这个resolution(分辨率)是可以自定义的，当我们的样本细胞数较大时候resolution 要高一些，一般情况2万细胞以上都是大于1.0的
查看每一类有多少个细胞

table(scRNA1@meta.data$seurat_clusters)
metadata <- scRNA@meta.data
cell_cluster <- data.frame(cell_ID=rownames(metadata), cluster_ID=metadata$seurat_clusters)
write.csv(cell_cluster,'cluster/cell_cluster.csv',row.names = F)

可视化降维有两个方法tSNE和UMAP

非线性降维——这个目的是为了可视化，而不是特征提取（PCA），虽然它也可以用来做特征提取。

tSNE

scRNA1 = RunTSNE(scRNA1, dims = pc.num)
embed_tsne <- Embeddings(scRNA1, 'tsne')
write.csv(embed_tsne,'embed_tsne.csv')
plot1 = DimPlot(scRNA1, reduction = "tsne") 
##画图
plot1
###label = TRUE把注释展示在图中
DimPlot(scRNA1, reduction = "tsne",label = TRUE) 
###你会发现cluster都标了图中
ggsave("tSNE.pdf", plot = plot1, width = 8, height = 7)
##把图片保存一下

UMAP---第二种可视化降维

scRNA1 <- RunUMAP(scRNA1, dims = pc.num)
embed_umap <- Embeddings(scRNA1, 'umap')
write.csv(embed_umap,'embed_umap.csv') 
plot2 = DimPlot(scRNA1, reduction = "umap") 
plot2
ggsave("UMAP.pdf", plot = plot2, width = 8, height = 7)

合并tSNE与UMAP

plotc <- plot1+plot2+ plot_layout(guides = 'collect')
plotc
ggsave("tSNE_UMAP.pdf", plot = plotc, width = 10, height = 5)

保存数据这节课的数据

saveRDS(scRNA1, file="scRNA1.rds")