library(dplyr)
library(Seurat)

 # 10X的数据可以使用Read10X这个函数，会返回一个UMI count矩阵，其中的每个值表示每个基因(行)在每个细胞(列)的分子数量
pbmc.data <- Read10X(data.dir = "D:/生信学习/猪年年尾R学习笔记/Seurat/filtered_gene_bc_matrices/hg19")
# class(pbmc.data)
# str(pbmc.data)


# 利用这个UMI count矩阵构建Seurat对象，(使用非标准化的原始数据先构建一个对象);
# Seurat对象是一个容器，里面装了数据(比如表达矩阵)和分析结果(比如PCA、聚类)。关于seurat对象的详细解释，见：    https://github.com/satijalab/seurat/wiki/Seurat#slots
# 在CreateSeuratObject时，会自动计算unique基因数量和总分子数，可以在object meta data中找到。
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc  #结果有个active assay，它存的就是表达矩阵，用pbmc[["RNA"]]@counts就可以访问
str(pbmc)  #可以利用str来查看seurat对象的信息。
pbmc[[]]  #使用[[ ]]是对assay的快速访问。
pbmc@meta.data  #通过object@meta.data也可以访问。
pbmc@version  #如果要查看其他的信息，可以用@
pbmc[["RNA"]]@counts  # 
#关于seurat对象还需要进一步学习！


# ------------------------------------------------------------------------
# 比较稀疏矩阵和矩阵占空间大小（可无！） ----------------------------------------------------------
# # What does data in a count matrix look like?
pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30]  #这个UMI count矩阵长什么样？和我们平常见到的bulk转录组矩阵一样吗？找几个基因，看看前30个细胞的情况
#从结果可以发现它是sparse Matrix，很杂乱。用 . 表示空着，数值为0，即没有检测到该分子。因为单细胞数据和bulk转录组数据还不太一样，scRNA的大部分数值都是0，原因一是由于建库时细胞的起始反应模板浓度就很低；二是测序存在dropout现象(本来基因在这个细胞有表达量，但没有测到)

# Seurat采用了一种聪明的再现方式，比原来用0表示的矩阵大大减小了空间占用，可以对比一下：
dense.size <- object.size(as.matrix(pbmc.data))  #as.matrix是把一个非matrix的变量变成matrix
dense.size
sparse.size <- object.size(pbmc.data)
sparse.size
dense.size/sparse.size  #比较稀疏矩阵跟矩阵之间占内存的大小，发现稀疏矩阵更省空间；
# -----------------------------------------------------------------------


# Standard pre-processing workflow ----------------------------------------
#以下为Seurat中scRNA-seq数据的标准预处理流程。包括：基于QC指标选择和过滤细胞；数据的normalization 和 scaling；寻找高异质性基因。

#   * * * 质控（QC）及筛选细胞 -----------------------------------------------------
# Seurat常用的QC指标：
# 1）* 每个细胞中检测到的unique基因的数量。（因为劣质细胞或空droplets基因数很少；一个droplet中有多个细胞基因数会很多；）
# 2）* Similarly, the total number of molecules detected within a cell (correlates strongly with unique genes) ？？？
# 3）* map到线粒体基因组上的reads比例。 （因为低质量或死亡的细胞会有大量线粒体污染；）
# 这篇文章 https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4758103/ 介绍了一些QC的标准

# Where are QC metrics stored in Seurat? ----------------------------------
# 质控筛选,通过[[ ]]向对象添加列，这里[[ ]]符号在pbmc的metadata这一子集里面添加新列，利用这个方法还可以挑出其他的feature；把表达矩阵里的线粒体QC单独存放在Seurat对象里，等于添加一列到metadata的对象里，

head(pbmc@meta.data, 5) #这时候看一下发现只有三列
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")  #PercentageFeatureSet函数可以计算来自某一个feature的reads count数量。我们在这里需要做的是挑出MT-开头的feature对应的所有基因,因为以 MT- 开头的基因是线粒体基因。然后将线粒体基因数量添加进 pbmc 中，成为单独的一列；
# # 检查下metadata， 查看QC指标存储的地方。
head(pbmc@meta.data, 5)  #如果没有上一步， pbmc[["percent.mt"]]则就不会有percent.mt这一列！meta.data里面多了一列内容。

# 对QC指标可视化，并依QC指标过滤细胞；# unique基因数大于2500或者小于200的都去掉，mitochondrial counts >5%的也去掉。
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

# 还可以用FeatureScatter函数画散点图来可视化两组feature信息的相关性，可以使用Seurat对象里的任何东西，如对象中的列、PC分数等。
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")  # nCount_RNA就是UMI counts数，nFeature_RNA就是指基因数
CombinePlots(plots = list(plot1, plot2))  #这个图的结果说明有表达量的reads和在线粒体的reads数量一般成反比，线粒体reads占比很高的情况一般有表达量的很少；相反，如果是真正表达的基因reads，很少会来自线粒体基因组；(右边????)就是刚才所说的比对结果中的基因数量(feature)一般会和测序得到的reads count值成正比


# 细胞过滤 --------------------------------------------------------------------
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)  # 过滤掉不想要的细胞，

# 数据归一化（normalize） --------------------------------------------------------
# 为何normalize要叫“归一化”？
# 不必纠结名称，原理很好理解：因为我们关心的是差异表达基因，即一个基因在不同样本(不同细胞)的差异，但这种差异往往不能直接比较，因为不同样本的测序深度和文库大小可能不同(比如单细胞建库时每个EP管加的引物，酶之类的可能不同)，那么基因表达量差异就存在一部分因素来自这里。
#为了更真实反映基因差异的生物学意义，就要将数据进行一个转换（例如RPKM、TPM等），可以让同一基因在不同样本中具有可比性。 “归一”归的就是这个起跑线。另外，原始表达量是一个离散程度很高的值，
#有的高表达基因表达量可能成千上万，可有的却只有几十。即使采用了一些归一化的算法消除了一些技术性误差，但真实存在的表达量极差往往会影响之后绘制热图、箱线图的美观性，#于是可以另外采用log 进行一个区间缩放（却不会影响真实值），
#比如原来表达量为1的，用log2(1+1)转换后结果依然为1；原来表达量为1000的，log2(1000+1)转换后约为10。那么原来相差1000倍的变化，现在只差10倍，在不破坏原有数据差异的同时，使数据更加集中。
#关于Seurat归一化，可以看这一篇：https://www.biorxiv.org/content/biorxiv/early/2019/03/18/576827.full.pdf
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)  # 使用pbmc <- NormalizeData(pbmc)的效果是一样的，因为默认是进行一个全局的LogNormalize操作
# 得到的结果存储在：pbmc[["RNA"]]@data


# 鉴定差异基因 ------------------------------------------------------------------
#在细胞与细胞间进行比较，选择表达量差别最大的(即同一个基因在有的细胞中表达量很高，同时在部分细胞中表达很少)，利用FindVariableFeatures函数，会计算一个mean-variance结果，也就是给出表达量均值和方差的关系并且得到top variable features。
# 计算top variable features有三种算法的，分别是：vst(默认)、mean.var.plot、dispersion
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

top10 <- head(VariableFeatures(pbmc), 10)  # Identify the 10 most highly variable genes

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(pbmc)   #绘制不带基因名的火山图
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)   #绘制带基因名的火山图
CombinePlots(plots = list(plot1, plot2))


# 数据标准化（scale） ------------------------------------------------------------
# scale的目的是实现数据的线性转换(scaling)，这也是降维处理(如PCA)之前的标准预处理。主要利用了ScaleData函数。前面的归一化（normalize）是log处理，它是对所有基因的表达量统一对待的，最后放在了一个起跑线上。但是为了真正找到差异，我们还要基于这个起跑线，考虑不同样本对表达量的影响。
# scale使得每个基因在所有细胞中的表达量均值为0，方差为1
# 这一步主要目的就是下一步的降维做准备，使用全部基因是比较慢。如果想提高速度，可以只对鉴定的HVGs(高异质性基因, 之前FindVariableFeatures设置的是2000个）进行操作，其实这个函数默认就是对部分差异基因进行操作。
# 这样操作的结果中降维和聚类不会被影响，但是只对HVGs基因进行标准化，下面画的热图依然会受到全部基因中某些极值的影响。所以为了热图的结果，还是对所有基因进行标准化比较好。
# ScaleData这个函数有两个默认参数：do.scale = TRUE和do.center = TRUE，然后需要输入进行scale/center的基因名（默认是取前面FindVariableFeatures的结果）。center就是对每个基因在不同细胞的表达量都减去均值；scale就是对每个进行center后的值再除以标准差（就是进行了一个z-score的操作）.

# * * * 使用全部基因进行scale ----------------------------------------------------
all.genes <- rownames(pbmc)
length(all.genes)
pbmc <- ScaleData(pbmc, features = all.genes)  #结果存储在 pbmc[["RNA"]]@scale.data 中

# * * * 只对鉴定的HVGs进行scale -------------------------------------------------
pbmc <- ScaleData(pbmc) 

#  * * * 移除不想要的差异来源 -------------------------------------------------------
# 怎样移除不想要的差异来源？比如说不想线粒体污染导致的差异影响整体差异分布，可以按下列方式进行scale
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")


# 线性降维 --------------------------------------------------------------------
# 降维的真实目的是尽可能减少具有相关性的变量数目，例如原来有700个样本，也就是700个维度/变量，但实际上根据样本中的基因表达量来归类，它们或多或少都会有一些关联，
# 这样有关联的一些样本就会聚成一小撮，表示它们的”特征距离“相近。最后直接分析这些”小撮“，就用最少的变量代表了实际最多的特征
# 使用scale后的数据进行 PCA降维。它默认选择之前鉴定的差异基因（2000个）作为input，但可以使用features进行设置；默认分析的主成分数量是50个
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc)) # 这里就使用默认值，结果保存在reductions 这个接口中。

# 检查下PCA结果
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)  #或者用 print(pbmc@reductions$pca, dims = 1:3, nfeatures = 5)

#通过 VizDimReduction, DimPlot和 DimHeatmap 对PCA结果进行可视化；
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")

DimPlot(pbmc, reduction = "pca")   #提取两个主成分（默认前两个，可以修改dims选项）绘制散点图。 这个图中每个点就是一个样本，它根据PCA的结果坐标进行画图，这个坐标保存在了cell.embeddings中：
# head(pbmc[['pca']]@cell.embeddings)[1:2,1:2]
# #还支持定义分组：根据pbmc@meta.data中的ident分组
# table(pbmc@meta.data$orig.ident)  #当然这里只有一个分组

DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)  #每个细胞和基因都根据PCA结果得分进行了排序，默认画前30个基因（nfeatures设置），1个主成分(dims设置)，细胞数没有默认值（cells设置）；其中balanced表示正负得分的基因各占一半
#也可以 DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)


# 确定数据集的维数 ----------------------------------------------------------------
# 既然每个主成分都表示具有相关性的一小撮变量（称之为”metafeature“，元特征集），那么问题来了：降维后怎么选择合适数量的主成分来代表整个数据集？
# 可以通过JackStraw或者碎石图来判断；# JackStrawPlot更耗时，碎石图更快。

#  * * * JackStraw法 -------------------------------------------------------
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot(pbmc, dims = 1:15)

#  * * * 碎石图法 -------------------------------------------------------------
ElbowPlot(pbmc)  #它是根据每个主成分对总体变异水平的贡献百分比排序得到的图，我们主要关注“肘部”的PC，它是一个转折点，一般认为拐点前的 PC(主成分)可以比较好地代表总体变化。
# 开始不确定时要多选一些主成分，选太少后续分析可能会出错。



# 细胞聚类 --------------------------------------------------------------------
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)  #在细胞数量为3000左右时设为0.4-1.2 会有比较好的结果
# Look at cluster IDs of the first 5 cells
head(Idents(pbmc), 5)  ## 结果聚成几类，用Idents查看
table(Idents(sce)) #也可以查看分成了几类。



#-----------------------------------------------------------------------
# 线性降维PCA的一个特点就是：它将高维中不相似的数据点放在较低维度时，会尽可能多地保留原始数据的差异，因此导致这些不相似的数据相距甚远，大多的数据重叠放置
# tSNE兼顾了高维降维+低维展示，降维后的那些不相似的数据点依然可以放得靠近一些，保证了点既在局部分散，又在全局聚集
# 后来开发的UMAP相比tSNE又能展示更多的维度，但是tSNE的图更好看一些(参考文献：https://sci-hub.tw/https://doi.org/10.1038/nbt.4314)

# 另一种降维方法—非线性降维（UMAP/tSNE） ------------------------------------------------
pbmc <- RunUMAP(pbmc, dims = 1:10)  #UMAP是基于Python的，如果要用它，必须先保证有一个python的安装环境
DimPlot(pbmc, reduction = "umap")  # 还可以设置label = TRUE让数字显示在每个cluster上

saveRDS(pbmc, file = "C:/Users/HASEE/Desktop/Seurat/pbmc_tutorial.rds") #对计算过程很费时的结果保存一下

# 找差异表达基因(cluster biomarkers) ---------------------------------------------
# 找差异表达基因的目的是根据表达量不同这个特征而分出不同细胞群的基因们（就是找具有生物学意义的HVGs，即biomarkers），marker可以理解成标记，就是一看到有这些基因，就说明是这个细胞群体。
# 利用FindAllMarkers()函数对单独的一个细胞群与其他几群细胞进行比较，找到正、负表达biomarker（这里的正负有点上调、下调基因的意思；正marker表示在我这个cluster中表达量高，而在其他的cluster中低）

cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)  ## 找到cluster1中的marker基因。min.pct的意思是：一个基因在任何两群细胞中的占比最低不能低于多少，当然这个可以设为0，但会大大加大计算量
head(cluster1.markers, n = 5)

cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)  #要比较cluster5和cluster0、cluster3的marker基因：
head(cluster5.markers, n = 5)


pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)# 也可以利用FindAllMarkers（对所有cluster都比较一下，并只挑出来正表达marker）
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)   ## 这一步过滤好好理解！(进行了分类、排序、挑前2个)
cluster1.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))



# 多种可视化方法 -----------------------------------------------------------------
# VlnPlot、FeaturePlot、DoHeatmap
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)  #用小提琴图对某些基因在全部cluster中表达量进行绘制

FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A")) #最常用的可视化方法，将基因表达量投射到降维聚类结果中

top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)  #对给定细胞和基因绘制热图，例如对每个cluster绘制前20个marker
DoHeatmap(pbmc, features = top10$gene) + NoLegend()


# 赋予每个cluster细胞类型 ---------------------------------------------------------
# 重点在于背景知识的理解，能够将marker与细胞类型对应起来，例如这里从各个cluster中找到的marker对应到了不同的细胞（这一过程是比较考验课题理解程度的）
new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "B", "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()

saveRDS(pbmc, file = "C:/Users/HASEE/Desktop/Seurat/pbmc3k_final.rds")