Seurat应用笔记 数据为Seurat官网示范数据
library(Seurat)
library(dplyr)
library(magrittr)
每次重新进入R都要重新library包
每次得到什么中间结果可以
save.image("保存路径")
或者直接点environment左上方的保存按钮
（当然可以一开始用setwd设定好工作位置，然后调用save.image的时候就可以直接用“文件名+.RData”的形式），就可以对每次的数据分开保存。
load的作用就是把一个RData类型文件载入到environment中，这样变量就全部在当前环境下，可以直接使用
setwd("运行环境位置")
单细胞数据文件分类：
由单细胞得的10x数据，怎么处理？
code1：
matrix <- Read10X_h5("filtered_gene_bc_matrices_h5.h5")
该函数直接读取10x的h5文件，并且得到一个matrix矩阵
code2：
matrix <- Read10X("data/matrix_dir")
默认读取运行环境中的 data/matrix_dir目录下包含文件，如果要读取运行环境之外的文件，就要从头什么盘开始输入位置。
包括matrix.mtx, genes.tsv, and barcodes.tsv 把他们都放到matrix中
导入文件后用Seurat包处理
pbmc <-CreateSeuratObject(raw.data = matrix,min.cells = 3,
min.genes = 200, project = "NAME")
创建一个叫做"pbmc"的东西，是SeuratObject。
其中命名raw.data为以上matrix（有的代码中命名不同，或者一开始导入的时候不以matrix形式导入，而以counts形式，总之这里raw.data = matrix可写“你要创建的文件的名字”=“你前面导入东西的名字”）；基因至少在min.cells个细胞中表达；每个细胞中至少表达min.genes个基因，可以自己设定；project的命名自己定义就好。
之后Seurat将数据保存在不同的slot中，如pbmc@raw.data, pbmc@data, pbmc@meta.data, pbmc@ident，其中raw.data存放的是每个细胞中每个gene的原始UMI数据，data存放的是gene的表达量，meta.data存放的是每个细胞的统计数据如UMI数目，gene数目等，ident此时存放的是project信息。
由于技术原因，一个GEM（单细胞测试的油包水滴）中可能会包含2个或多个细胞，也可能不包含细胞，这时候可以通过观察每个barcode中的基因数目或UMI数目来判断。
处理筛掉异常GEM的方法：
- 可以计算每个barcode中的线粒体基因含量等，从而更加仔细的观察数据的质量。
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^mt-")
mt-开头的是线粒体基因，这里将其进行标记并统计其分布频率
VlnPlot(object = pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")+ NoLegend()
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")+ NoLegend()
CombinePlots(plots = list(plot1, plot2))
以上对 pbmc 对象做小提琴图，分别为基因数，细胞数和线粒体占比
以上根据图片中基因数和线粒体数，分别设置过滤参数，这里基因数200-2500，线粒体百分比为小于 5%（常用的过滤参数）
作图，直观地看一下数据的质量。要输入代码plot或者Combineplot，才可以在右边出现图
标准化
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
鉴定高变基因，之后就关注这些基因进行研究
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
top10 <- head(VariableFeatures(pbmc), 10)
使用默认参数，用vst方法选取2000个高变基因。
然后选取top10变化的
plot3 <- VariableFeaturePlot(pbmc)
plot4 <- LabelPoints(plot = plot3, points = top10, repel = TRUE,xnudge=0,ynudge=0)
plot3+plot4
得到图，看到变化最大的前2000个和10个的图
接下来进行数据归一化
1.线性姜维
有一堆细胞，但是你需要将这些细胞进行分群，但是你又不知道这些细胞如果进行分群的话是按照什么依据进行分群，就要先使用PCA 降维找到细胞分群的主要特点，这也就是我们常说的主成分分析
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc,features = all.genes,vars.to.regress = "percent.mt")
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
数据归一化，对所有基因进行标准化，默认只是标准化高变基因（ 2000 ），速度更快，不影响 PCA 和分群，但影响热图的绘制。
线性降维（PCA），默认用前面得到的2000个高变基因集，但也可通过 features 参数自己指定 需要挺长时间的，从1%开始加载
定义可视化细胞和功能的几种有用的方式PCA，包括VizDimReduction，DimPlot，和DimHeatmap 得到一堆positive negative基因
VizDimLoadings(object = pbmc, dims = 1:2, reduction = "pca")
DimPlot(pbmc, reduction = "pca")+ NoLegend()
DimHeatmap(pbmc, dims = 1:2, cells = 500, balanced = TRUE)
以上三行得到三个图，定义可视化细胞和功能，还不知道有啥用？
鉴定数据集的可用维度，确定哪些主成分所代表的基因可以进入下游分析用于后续细胞分类。
这里可以使用JackStraw做重抽样分析。可以用JackStrawPlot可视化看看哪些主成分可以进行下游分析。
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot(pbmc, dims = 1:15)
ElbowPlot(object = pbmc)
第一步需要一定计算时间最后得到图JackStrawPlot和Elbowplot，后者比较直观。
第二行，得到JackStrawPlot。是彩色的，虚线以上的为可用维度，dim是展示多少个。可以调节dim，画出不同数量的 pca 查看。
也就是：基于每个主成分对方差解释率的排名。建议尝试选择多个主成分个数做下游分析，对整体影响不大；在选择此参数时，建议选择偏高的数字（ “宁滥勿缺”，为了获取更多的稀有分群）；有些亚群很罕见，如果没有先验知识，很难将这种大小的数据集与背景噪声区分开来。
第三行，得到Elbowplot，点图。比较直观地看到，前面的几个点分得比较开。
之后就可以决定选取多少主成分用于细胞分类。从上到下，从p-value最小的开始选择top few。
尊重官网的建议，这里选取top 10
- 非线性降维（ UMAP/tSNE)
TSNE降维与PCA降维并不相同，TSNE降维主要是非线性降维。
TSNE优势：使得相似的对象有更高的概率被选择，而不相似的对象有较低的概率被选择。它将多维数据映射到适合于人类观察的两个或多个维度，广泛应用于图像处理。
一般情况下，在进行细胞分群时：先进行PCA主成分分析，再进行TSNE降维分析进行细胞分群，这样的分群结果可靠度更高。
Seurat 提供了几种非线性降维的方法进行数据可视化（在低维空间把相似的细胞聚在一起），比如 UMAP 和 t-SNE，运行 UMAP 需要先安装'umap-learn'包，两种方法都可以使用，但不要混用，这样，后面的结算结果会将先前的聚类覆盖掉，只能保留一个。
2.1 umap
umap优于tsne：
UMAP 与 t-SNE 均是非线性降维，多用于数据可视化
UMAP 结构与t-SNE一致
UMAP 计算更快
UMAP 能更好地反映高纬结构，比t-SNE有着更好的连续性
这种连续性反映到单细胞分析中就是能更好滴可视化细胞的分化状态（UMAP better represents the multi-branched continuous trajectory of hematopoietic development）
umap安装包如下：
install(packages = 'umap-learn')
基于PCA 空间中的欧氏距离计算 nearest neighbor graph，优化任意两个细胞间的距离权重（输入上一步得到的 PC 维数）
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)
resolution 参数决定下游聚类分析得到的分群数，对于 3K 左右的细胞，设为 0.4-1.2 能得到较好的结果(官方说明)；如果数据量增大，该参数也应该适当增大；增加的值会导致更多的群集。
使用 Idents（）函数可查看不同细胞的分群；查看每一类有多少个细胞
head(Idents(pbmc), 8)
table(pbmc@active.ident)
得到前8个细胞分群的table
pbmc <- RunUMAP(object = pbmc, dims = 1:10)
DimPlot(object = pbmc, reduction = "umap")
umapplot<-UMAPPlot(pbmc,label = TRUE, pt.size = 1.5)+ NoLegend()
这里用 DimPlot()函数绘制散点图。如果只做tsne，参数reduction = "tsne"；如果不设定reduction，默认先从搜索 umap， 然后 tsne, 再然后 pca；
也可以直接使用这 3 个函数 PCAPlot()、 TSNEPlot()、UMAPPlot()； cols， pt.size 分别调整分组颜色和点的大小；
2.2 TSNE
pbmc <- RunTSNE(pbmc, dims = 1:10)
DimPlot(object = pbmc, reduction = "tsne")
tsneplot<-TSNEPlot(pbmc,label = TRUE, pt.size = 1.5)+ NoLegend()
做完以上降维处理保存数据
saveRDS(pbmc, file = "pbmc_tutorial.rds")
save(pbmc,file="pbmc.RData")
load(file = "pbmc.RData")
保存在了工作路径下文件夹中，并且下次用的时候直接load就好
寻找差异表达的特征（聚类生物标志物）
寻找某个聚类和其他所有聚类显著表达的基因
cluster1.markers <- FindMarkers(object = pbmc, ident.1 = 1, min.pct = 0.25)
head(x = cluster1.markers, n = 5)
find all markers distinguishing cluster 5 from clusters 0 and 3
cluster5.markers <- FindMarkers(object = pbmc, ident.1 = 5, ident.2 = c(0,3), min.pct = 0.25)
write.table(as.data.frame(cluster5.markers), file="cluster1vs2.xls", sep="\t", quote = F)
find markers for every cluster compared to all remaining cells, report only the positive ones
pbmc.markers <- FindAllMarkers(object = pbmc, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25)
pbmc.markers <- FindAllMarkers(object = pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25, return.thresh = 0.01)
library(dplyr)
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(object = pbmc, features = top10$gene) + NoLegend()
画图
VlnPlot(object = pbmc, features = c("MS4A1", "CD79A"))
FeaturePlot(object = pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
you can plot raw UMI counts as well
VlnPlot(object = pbmc, features.plot = c("NKG7", "PF4"), use.raw = TRUE, y.log = TRUE)
FeaturePlot(object = pbmc, features = c("hbaa1","hbaa2","ighv1-4","cd79a","cd79b","cd4-1","cd8a","cd8b","itga2b"), cols = c("grey", "blue"), reduction = "tsne")
将单元类型标识分配给集群
new.cluster.ids <- c("1", "3", "4", "2", "5", "6", "7", "8")
names(x = new.cluster.ids) <- levels(x = pbmc)
pbmc <- RenameIdents(object = pbmc, new.cluster.ids)
DimPlot(object = pbmc, label = TRUE, pt.size = 1.5) + NoLegend()
气泡图
markers.to.plot <- c("cd74a", "cd74b")
DotPlot(pbmc, features = markers.to.plot) + RotatedAxis()
