寻找高变基因

##寻找高变基因
library(Seurat)
library(tidyverse)
library(patchwork)
library(dplyr)
rm(list=ls())
#加载数据
scRNA <-load("scRNA1.Rdata")
###Identification of highly variable features (feature selection)，找到高变基因
###官方推荐是2000个高变基因，很多文章也有设置30000的，这个因自己的实验项目决定
scRNA1 <- FindVariableFeatures(scRNA1, selection.method = "vst", nfeatures = 2000) 
# 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 

横坐标是某基因在所有细胞中的平均表达值，纵坐标是此基因的方差。红色的点是被选中的高变基因，黑色的点是未被选中的基因，变异程度最高的10个基因在如图中标注了基因名称。

在进行PCA降维之前还有要对数据进行中心化（必选）和细胞周期回归分析（可选），它们可以使用ScaleData()函数一步完成。

数据中心化

#数据中心化
##如果内存足够最好对所有基因进行中心化
scale.genes <-  rownames(scRNA1)
scRNA1 <- ScaleData(scRNA1, features = scale.genes)
##如果内存不够，可以只对高变基因进行标准化
scale.genes <-  VariableFeatures(scRNA1)
scRNA1 <- ScaleData(scRNA1, features = scale.genes)

##scRNA对象中原始表达矩阵经过标准化和中心化之后，已经产生了三套基因表达数据，可以通过以下命令获得
#原始表达矩阵
GetAssayData(scRNA1,slot="counts",assay="RNA") 
#标准化之后的表达矩阵              
GetAssayData(scRNA1,slot="data",assay="RNA")
#中心化之后的表达矩阵 
GetAssayData(scRNA1,slot="scale.data",assay="RNA") 

细胞周期回归

上一步找到的高变基因，常常会包含一些细胞周期相关基因。它们会导致细胞聚类发生一定的偏移，即相同类型的细胞在聚类时会因为细胞周期的不同而分开。如在Spatially and functionally distinct subclasses of breast cancer-associated fibroblasts revealed by single cell RNA sequencing 中 vCAFs和cCAFs中只有cell cycle genes的表达差异。
查看我们选择的高变基因中有哪些细胞周期相关基因：

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)
colnames(scRNA1@meta.data)
#[1] "orig.ident"   "nCount_RNA"   "nFeature_RNA" "percent.mt"  
#[5] "percent.HB"   "S.Score"      "G2M.Score"    "Phase"       

以上代码运行之后会在scRNA1@meta.data中添加S.Score、G2M.Score和Phase三列有关细胞周期的信息。

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

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

##如果需要消除细胞周期的影响
#scRNAb <- ScaleData(scRNA1, vars.to.regress = c("S.Score", "G2M.Score"), features = rownames(scRNA1))

PCA降维并提取主成分

PCA降维

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

左图是根据主成分1和2的值将细胞在平面上展示出来，右图展示前20个主成分的解释量（pca中等同于标准差）。后续分析要根据右图选择提取的pc轴数量，一般选择斜率平滑的点之前的所有pc轴，此图作者的建议是选择前18个pc轴。

pc.num=1:18

获取PCA结果

#此部分代码为分析非必须代码，不建议运行！！！
#获取基因在pc轴上的投射值
Loadings(object = scRNA1[["pca"]])
#获取各个细胞的pc值
Embeddings(object = scRNA1[["pca"]])
#获取各pc轴解释量方差
Stdev(scRNA1)
#查看决定pc值的top10基因， 此例查看pc1-pc5轴
print(scRNA1[["pca"]], dims = 1:5, nfeatures = 10) 
#查看决定pc值的top10基因在500个细胞中的热图，此例查看pc1-pc9轴
DimHeatmap(scRNA1, dims = 1:9, nfeatures=10, cells = 500, balanced = TRUE)

细胞聚类

此步利用 细胞-PC值 矩阵计算细胞之间的距离，然后利用距离矩阵来聚类。其中有两个参数需要人工选择，第一个是FindNeighbors()函数中的dims参数，需要指定哪些pc轴用于分析，选择依据是之前介绍的cluster/pca.png文件中的右图。第二个是FindClusters()函数中的resolution参数，需要指定0.1-0.9之间的一个数值，用于决定clusters的相对数量，数值越大cluters越多。

#细胞聚类
###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) 
###这个分辨率是可以自定义的，当我们的样本细胞数较大时候resolution 要高一些，一般情况2万细胞以上都是大于1.0的
scRNA1 <- FindClusters(scRNA1, resolution = 0.5)
## 查看每一类有多少个细胞
table(scRNA1@meta.data$seurat_clusters)
 # 0   1   2   3   4   5   6   7   8   9 
#297 181 121 112  70  63  53  50  44  22 
metadata <- scRNA1@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
scRNA2= scRNA1
scRNA2 = RunTSNE(scRNA2, dims = pc.num)
embed_tsne <- Embeddings(scRNA2, 'tsne')
write.csv(embed_tsne,'embed_tsne.csv')
plot1 = DimPlot(scRNA2, reduction = "tsne") 
##画图
plot1
###label = TRUE把注释展示在图中
DimPlot(scRNA1, reduction = "tsne",label = TRUE) 
#运行这一步的时候出错了 。。暂时还不知道怎么解决。。解决了哈哈
plot1=DimPlot(scRNA1, reduction = "tsne",label = TRUE) 
###你会发现cluster都标了图中
ggsave("cluster/tSNE.pdf", plot = plot1, width = 8, height = 7)
##把图片保存一下

#UMAP---第二种可视化降维
scRNA3=scRNA1
scRNA3 <- RunUMAP(scRNA3, dims = pc.num)
embed_umap <- Embeddings(scRNA3, 'umap')
write.csv(embed_umap,'embed_umap.csv') 
plot2 = DimPlot(scRNA3, reduction = "umap") 
plot2
plot2 = DimPlot(scRNA3, reduction = "umap",label = TRUE)
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")