目标：
参考：
SingleR单细胞类型自动注释与celldex参考数据包
说白了就是根据细胞的表达基因确定细胞类型。总的来说，两种方法：手动注释和软件注释。我在实际处理的过程中，主要还是手动注释，软件的注释结果只作为一个参考，来辅助证实手动注释的结果是准确无误的。
相信看过几篇单细胞文献的小伙伴基本都会知道几种常见细胞大类的marker，我们在注释的时候用这些marker就可以基本确定细胞大类了。另外，以SingleR为例，对于细胞大类的注释还算准确，当然也没有到很准确的程度，我试过更换SingleR里面不同的参考集，最后得到的大类注释结果一致性不到80%。对于细胞小类的注释，软件就更加无法胜任了，我几乎没有看过高分的文献会用SingleR的小类注释结果。当然不排除随着单细胞研究的普及和深入，以后能有更准确的软件出现。
人工注释
人工注释需要借助文献检索marker或者结合常用的注释数据库，例如两个常用的数据库：
CellMarker（http://bio-bigdata.hrbmu.edu.cn/CellMarker/）；
panglaoDB（A Single Cell Sequencing Resource ForGene Expression Data，https://panglaodb.se/）。
人工注释比较适合有经验的科研工作者，但随着单细胞的研究越来越多，可提供给我们的细胞类型的marker信息也越来越丰富，人工注释比较耗费精力，优点在于准确性相对较好。
手动注释的话，用小提琴图、热图展示
SingleR
使用内置参考进行注释（最简便的）
使用SingleR的最简单方法是使用内置参考对细胞进行注释。通过专用的检索功能提供了7个参考数据集（主要来自大量RNA-seq或微阵列数据）。
singleR自带的7个参考数据集，需要联网才能下载，其中5个是人类数据，2个是小鼠的数据：
BlueprintEncodeData Blueprint (Martens and Stunnenberg 2013) and Encode (The ENCODE Project Consortium 2012) （人）
DatabaseImmuneCellExpressionData The Database for Immune Cell Expression(/eQTLs/Epigenomics)(Schmiedel et al. 2018)（人）
HumanPrimaryCellAtlasData the Human Primary Cell Atlas (Mabbott et al. 2013)（人）
MonacoImmuneData, Monaco Immune Cell Data - GSE107011 (Monaco et al. 2019)（人）
NovershternHematopoieticData Novershtern Hematopoietic Cell Data - GSE24759（人）
ImmGenData the murine ImmGen (Heng et al. 2008) （鼠）
MouseRNAseqData a collection of mouse data sets downloaded from GEO (Benayoun et al. 2019).鼠）
SingleR包的工具包括：
SingleR()函数，用于注释
plotScoreHeatmap(）使用热图比较label score
pruneScores() 对低质量的注释进行“修剪”（即注释为NA）
plotScoreDistribution()观察修剪是否合适，默认为偏离中值3个MAD（评价数据离散程度的统计学指标）
一、SingleR自动注释流程
1、准备输入数据
同一格式要求：
（1）矩阵(matrix)/稀疏矩阵(dgCMatrix)/数据框(data.frame)均可，或者是SummarizedExperiment对象（默认指定为logcounts slot）；
（2）必须是经标准化并log转换的表达矩阵，对应Seurat对象的data slot。未知细胞类型、待注释的单细胞表达矩阵
SingleR包的安装
使用devtools包进行安装，如果安装失败，可下载进行本地安装，需提前安装好报错的各种依赖包。
# 使用devtools包进行安装
devtools::install_github('dviraran/SingleR')
# this might take long, though mostly because of the installation of Seurat.
# 安装celldex
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("celldex")
1、选择/下载参考数据集
下载一个参考数据集，并保存在一个RDS文件中。
可视化：
plotScoreHeatmap(pred)
plotDeltaDistribution(pred,size=0.5)
热图的每一列是细胞，
从头开始的细胞类型注释
参考：
单细胞分析实录(5): Seurat标准流程
单细胞分析实录(6): 去除批次效应/整合数据
单细胞分析实录(7): 差异表达分析/细胞类型注释
SingleR如何使用自定义的参考集
看了很多教程，都是从聚类开始的，TOP大佬流程完整，这次跟着TOP大佬从头做。
数据是TOP给的看TOP的原帖百度云下载，暂时不知道是哪篇文章的。下面就是先跑一遍Seurat标准流程。
######## 1. 导入数据，创建Seurat对象
library(Seurat)
library(tidyverse)
testdf=read.table("test_20210105.txt",header = T,row.names = 1)
test.seu=CreateSeuratObject(counts = testdf)
test.seu
#An object of class Seurat
#33538 features across 6746 samples within 1 assay
#Active assay: RNA (33538 features)
#①1 assay表示有一套数据，假如用Seurat里面的函数去过批次效应，
#这里会有2个assay，另外一个是去批次（整合）之后的；
#②cell hashing的tag表达矩阵生成Seurat对象，这时的assay为"HTO"，不叫"RNA"；
#其他情况类似
#测试数据有33538个基因，6746个细胞。除此之外，还要关注一下另外两层信息：
#test.seu@meta.data这个数据框用来存储元数据，每一个细胞都有多个属性；
#test.seu[["RNA"]]@counts这个稀疏矩阵用来存储原始UMI表达矩阵。
head(test.seu@meta.data)
orig.ident nCount_RNA nFeature_RNA
A_AAACCCAAGGGTCACA A 3714 1151
A_AAACCCAAGTATAACG A 1855 816
A_AAACCCAGTCTCTCAC A 1530 823
A_AAACCCAGTGAGTCAG A 11145 1087
A_AAACCCAGTGGCACTC A 2289 834
A_AAACGAAAGCCAGAGT A 3714 990
#我这里CB的前面人为加上了样本来源，用下划线连接，orig.ident是自动识别得到的
test.seu[["RNA"]]@counts[1:4,1:4]
######## 2. 简单过滤
#接下来，我们根据每个细胞内部线粒体基因表达占比、检测到的基因数、检测的UMI总数这三个方面来对细胞进行简单的过滤。
#先计算细胞内线粒体基因表达占比，类似的核糖体基因(大多为RP开头)也能这样计算，
#还要注意不要将线粒体基因的MT-写成了MT，不然就把别的基因也算进去了：
test.seu[["percent.mt"]] <- PercentageFeatureSet(test.seu, pattern = "^MT-")
#正则表达式，表示以MT-开头；test.seu[["percent.mt"]]这种写法会在meta.data矩阵加上一列
#这里我已经根据预先设定好的阈值过滤了，代码如下
test.seu <- subset(test.seu, subset = nCount_RNA > 1000 &
nFeature_RNA < 5000 &
percent.mt < 30 &
nFeature_RNA > 600)
#过滤之后的数值分布如下，用到VlnPlot()函数，该绘图函数里面的feature参数可以是meta.data矩阵的某一列，
#也可以是某一个基因，很多文章都用这种图展示marker gene
VlnPlot(test.seu,features = c("nCount_RNA", "nFeature_RNA", "percent.mt"),pt.size = 0)
#nFeature_RNA/nCount_RNA不能太小（空液滴），不能太大（doublet、测序技术限制）
######## 3. LogNormalize消除文库大小的影响
test.seu <- NormalizeData(test.seu, normalization.method = "LogNormalize", scale.factor = 10000)
#标准化之后的矩阵存储在test.seu[["RNA"]]@data
test.seu[["RNA"]]@data[1:4,1:4]
######## 4. 找Variable基因
#因为单细胞表达矩阵很稀疏（很多0），选high variable基因的目的可以找到包含信息最多的基因（很多基因的表达差不多都是0），
#同时极大提升软件运行速度
test.seu <- FindVariableFeatures(test.seu, selection.method = "vst", nfeatures = 2000)
#这些基因存储在VariableFeatures(test.seu)，有时候可能需要人为指定high variable基因，可以这样：
VariableFeatures(test.seu)="specific genes"
######## 5. scale归一化表达矩阵
#（基于前面得到的data矩阵）这一步之后，所有基因的表达值的分布就差不多了，不然表达值不在一个数量级，
#对后续降维聚类影响挺大。新的矩阵存储在test.seu[["RNA"]]@scale.data里面。
test.seu <- ScaleData(test.seu, features = rownames(test.seu))
#默认只对上一步选出来的基因scale，这里调整为所有基因，
#是为了方便以后画热图（画热图一般会用scale之后的z-score）
######## 6. 降维聚类
#基于前面得到的high variable基因的scale矩阵
test.seu <- RunPCA(test.seu, npcs = 50, verbose = FALSE)
test.seu <- FindNeighbors(test.seu, dims = 1:30)
test.seu <- FindClusters(test.seu, resolution = 0.5)
test.seu <- RunUMAP(test.seu, dims = 1:30)
test.seu <- RunTSNE(test.seu, dims = 1:30)
#Run开头的函数降维，Find开头的函数聚类，一般就这几步，相对固定。
#PCA将原来2000维的数据降到50维，dims参数表示使用多少个主成分（一般20左右就可以了，多几个少几个对结果影响不大），
#resolution参数表达聚类的分辨率，这个值大于0，一般都是在0-1范围里面调整，越大得到的cluster越多，
#这个值可以反复调整，并不会改变降维的结果（也就是tsne、umap图的二维坐标）。
#这一步之后的数据是这样的
test.seu
#An object of class Seurat
#33538 features across 6746 samples within 1 assay
#Active assay: RNA (33538 features)
#3 dimensional reductions calculated: pca, umap, tsne
# 几种降维方式都会呈现出来
#聚类之后test.seu@meta.data多了两列，RNA_snn_res.0.5记录了你用的分辨率，
#最终的聚类结果保存在seurat_clusters中
> head(test.seu@meta.data)
orig.ident nCount_RNA nFeature_RNA percent.mt RNA_snn_res.0.5 seurat_clusters
A_AAACCCAAGGGTCACA A 3714 1151 9.585353 6 6
A_AAACCCAAGTATAACG A 1855 816 12.776280 0 0
A_AAACCCAGTCTCTCAC A 1530 823 14.248366 13 13
A_AAACCCAGTGAGTCAG A 11145 1087 2.853297 3 3
A_AAACCCAGTGGCACTC A 2289 834 15.640017 2 2
A_AAACGAAAGCCAGAGT A 3714 990 5.654281 0 0
######## 7. tsne/umap展示结果
library(cowplot)
test.seu$patient=str_replace(test.seu$orig.ident,"_.*$","")
p1 <- DimPlot(test.seu, reduction = "tsne", group.by = "patient", pt.size=0.5)
p2 <- DimPlot(test.seu, reduction = "tsne", group.by = "ident", pt.size=0.5, label = TRUE,repel = TRUE) #后面两个参数用来添加文本标签
p3 <- DimPlot(test.seu, reduction = "umap", group.by = "patient", pt.size=0.5)
p4 <- DimPlot(test.seu, reduction = "umap", group.by = "ident", pt.size=0.5, label = TRUE,repel = TRUE)
fig_tsne <- plot_grid(p1, p2, labels = c('patient','ident'),align = "v",ncol = 2)
ggsave(filename = "tsne.pdf", plot = fig_tsne, device = 'pdf', width = 30, height = 15, units = 'cm')
fig_umap <- plot_grid(p3, p4, labels = c('patient','ident'),align = "v",ncol = 2)
ggsave(filename = "umap.pdf", plot = fig_umap, device = 'pdf', width = 30, height = 15, units = 'cm')
#ident表示每个细胞的标签，聚类之后就是聚类的结果，在一些特定场景可以更换。
#在umap图中，cluster之间的距离更明显
#从上面的图可以看出不同样本其实是有批次效应的，下一讲我会介绍两种去批次效应的方法。
#这里补充一下什么是批次效应以及什么时候需要去批次，什么时候单独整合
#批次效应直观的体现在UMAP和tSNE图中，如果细胞群大致重合但又不是完全重叠，这就是批次效应的体现。
#去批次就是要把他们变得完全重合在一起。同类型组织才需要去批次，因为细胞组成或者趋势是相同的。
#但是如果是不同患者的肿瘤样本的话，肿瘤细胞本身异质性就很大，这样如果去批次的话，反而影响了细胞本身的异质性，
#这种情况下简单合并就好。
######## 8.去除批次效应/整合数据
#接下来我会介绍Seurat v3的标准整合流程、Seurat结合Harmony 的整合流程
8.去除批次效应/整合数据
这部分说是用上一次数据，但是我看到不太一样，准备还是用经典数据GSE139324里的两组数据GSM4138110和GSM4138111作为测试数据。参考：
1. Seurat v3的标准整合流程
对于不同样本先分别运行Seurat的标准流程到找Variable基因这一步
rm(list=ls())
library(dplyr)
library(Seurat)
library(patchwork)
#下面的10X标准三文件不是.gz格式，还需要解压成tsv。
#对于不同样本先分别运行Seurat的标准流程到找Variable基因这一步
### testA ----
setwd("D:/去批次T")
testA.seu <- Read10X(data.dir ="D:/去批次T/GSM4138110")
testA.seu <- CreateSeuratObject(counts = testA.seu, project = "testA.seu",min.cells = 3, min.features = 200)
testA.seu <- NormalizeData(testA.seu, normalization.method = "LogNormalize", scale.factor = 10000)
testA.seu <- FindVariableFeatures(testA.seu, selection.method = "vst", nfeatures = 2000)
### testB ----
setwd("D:/去批次T")
testB.seu <- Read10X(data.dir ="D:/去批次T/GSM4138111")
testB.seu <- CreateSeuratObject(counts = testB.seu, project = "testB.seu",min.cells = 3, min.features = 200)
testB.seu <- NormalizeData(testB.seu, normalization.method = "LogNormalize", scale.factor = 10000)
testB.seu <- FindVariableFeatures(testB.seu, selection.method = "vst", nfeatures = 2000)
#然后就是主要的整合步骤，object.list参数是由多个Seurat对象构成的列表，如下：
### Integration ----
testAB.anchors <- FindIntegrationAnchors(object.list = list(testA.seu,testB.seu), dims = 1:20)
testAB.integrated <- IntegrateData(anchorset = testAB.anchors, dims = 1:20)
#需要注意的是：上面的整合步骤相对于harmony整合方法，对于较大的数据集（几万个细胞），非常消耗内存和时间；
#当存在某一个Seurat对象细胞数很少（印象中200以下这样子），会报错，这时建议用第二种整合方法
#这一步之后就多了一个整合后的assay（原先有一个RNA的assay），整合前后的数据分别存储在这两个assay中
> testAB.integrated
An object of class Seurat
18411 features across 3019 samples within 2 assays
Active assay: integrated (2000 features, 2000 variable features)
1 other assay present: RNA
> dim(testAB.integrated[["RNA"]]@counts)
[1] 16411 3019
> dim(testAB.integrated[["RNA"]]@data)
[1] 16411 3019
> dim(testAB.integrated[["integrated"]]@counts) #因为是从RNA这个assay的data矩阵开始整合的，所以这个矩阵为空
[1] 0 0
> dim(testAB.integrated[["integrated"]]@data)
[1] 2000 3019
