1.下载和整理数据
示例https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE231920
需要3个文件：
GSM7306054_sample1_barcodes.tsv.gz
GSM7306054_sample1_features.tsv.gz
GSM7306054_sample1_matrix.mtx.gz
下下来是个压缩包
1.1 解包文件
untar("GSE231920_RAW.tar",exdir = "input")
1.2 单细胞文件组织的要求
这三个文件是10X的标准文件，需要放在同一个文件夹里，并且名字是固定的，不能有前缀。
“barcodes.tsv.gz”：存储的是barcodes，相当于细胞的编号，是表达矩阵的列名。
“features.tsv.gz”：存储的是基因名称，是表达矩阵的行名，也可以是”genes.tsv.gz”。
“matrix.mtx.gz”:存储的是每个位置的数值，是表达矩阵的内容，仅存储了非零的数值。
1.3 修改文件名称
library(stringr)
nn = str_remove(dir("input/"),"GSM7306054_sample1_")
nn
## [1] "barcodes.tsv.gz" "features.tsv.gz" "matrix.mtx.gz"
file.rename(paste0("input/",dir("input/")),
paste0("input/",nn))
## [1] TRUE TRUE TRUE
dir("input/")
## [1] "barcodes.tsv.gz" "features.tsv.gz" "matrix.mtx.gz"
2.读取并创建Seurat对象
2.1读取文件
rm(list = ls())
#清空右上角的所有变量，方便反复调试代码
library(Seurat)
library(patchwork)
library(tidyverse)
ct = Read10X("input/")
#读取标准文件，接受的参数是文件夹名称。文件夹里的三个文件合到一起才是完整的单细胞表达矩阵。
dim(ct)
## [1] 33538 10218
两个数字分别是行数(基因数)和列数(细胞数)。
2.2 稀疏矩阵
class(ct)
## [1] "dgCMatrix"
## attr(,"package")
## [1] "Matrix"
ct[c("CD3D","TCL1A","MS4A1"), 1:30]
## 3 x 30 sparse Matrix of class "dgCMatrix"
##
## CD3D . . 5 . . 1 . . 2 . 3 . . . 4 . . 3 . . . . . . . 5 3 . . .
## TCL1A . . . . . . . . . . . . 3 . . . . . . . . 1 1 . . . . . . .
## MS4A1 4 . . . . . . . . 4 . . 1 . 1 . 2 . . . . 2 4 7 5 . . . 1 .
稀疏矩阵是存储0值比较多的数据用的，用“.”表示0，可以节省空间。单细胞矩阵里面有大量的0值。
2.3 创建Seurat对象
seu.obj <- CreateSeuratObject(counts = ct,
min.cells = 3, #一个基因至少要在3个细胞里有表达，才被保留
min.features = 200) #一个细胞里至少要200个基因有表达，才被保留
dim(seu.obj)
## [1] 20648 10105
2.4 细胞抽样
!!!这一步是为了为了节省计算资源，我们减一下细胞的数量，实战时不能减!!!!
set.seed(1234) #set.seed(1234)是设计随机种子，作用是让后面的随机抽样变得可重复，只要多次运行时，setseed的括号里的数字没变，抽样的结果就不会变。
seu.obj = subset(seu.obj,downsample = 3000)
2.5 对象
广义的“对象”是Rstudio右上角所有的数据，向量数据框矩阵列表都是对象。
狭义的“对象”是由R包的作者定义的，以固定模式组织的数据，里面的结构都是规定好的。
可以用class查看
class(seu.obj)
## [1] "Seurat"
## attr(,"package")
## [1] "SeuratObject"
#意思是：这是一个Seurat对象，出自于SeuratObject这个包。
提取它的子集，有@和$两个符号,具体该用哪一个可以试试，一般第一层次都是@，后面的就不一定了
2.6 探索Seurat对象的meta信息
seu.obj@meta.data %>% head()
## orig.ident nCount_RNA nFeature_RNA
## AAACCCACAGTCGGTC-1 SeuratProject 4243 1256
## AAACGAAAGAATCGCG-1 SeuratProject 7307 2577
## AAAGAACAGCTTACGT-1 SeuratProject 8154 1881
## AAAGAACAGGTTCATC-1 SeuratProject 8223 2182
## AAAGAACAGTCTGTAC-1 SeuratProject 3884 1377
## AAAGAACTCCACCTCA-1 SeuratProject 3997 1307
或
seu.obj[[]] %>% head()
## orig.ident nCount_RNA nFeature_RNA
## AAACCCACAGTCGGTC-1 SeuratProject 4243 1256
## AAACGAAAGAATCGCG-1 SeuratProject 7307 2577
## AAAGAACAGCTTACGT-1 SeuratProject 8154 1881
## AAAGAACAGGTTCATC-1 SeuratProject 8223 2182
## AAAGAACAGTCTGTAC-1 SeuratProject 3884 1377
## AAAGAACTCCACCTCA-1 SeuratProject 3997 1307
行名是barcodes,也就是表达矩阵里面的列名，相当于细胞的编号。
orig.ident是细胞的原始分类，如果是单样本数据，在前面CreateSeuratObject时，如果指定project参数，就会显示在这一列，整列都是同一个单词，因为我们前面没指定，所以显示默认值SeuratProject。
nCount_RNA是每个细胞中所有基因的表达量之和，也就是表达矩阵的列和。
nFeature_RNA是每个细胞中表达量不为零的基因的数量。
exp = seu.obj@assays$RNA$counts #这是表达矩阵
exp[1:4,1:4]
## 4 x 4 sparse Matrix of class "dgCMatrix"
## AAACCCACAGTCGGTC-1 AAACGAAAGAATCGCG-1 AAAGAACAGCTTACGT-1
## AL627309.1 . . .
## AL627309.3 . . .
## AL669831.5 . . .
## FAM87B . . .
## AAAGAACAGGTTCATC-1
## AL627309.1 .
## AL627309.3 .
## AL669831.5 .
## FAM87B .
sum(exp[,1])
## [1] 4243
table(exp[,1]!=0) #统计TRUE和FALSE的个数
##
## FALSE TRUE
## 19392 1256
3.质控
3.1 质控的指标及原因
这里是对细胞进行的质控，指标是：（1）线粒体基因含量不能过高；（2）nFeature_RNA 不能过高或过低
nFeature_RNA是每个细胞中检测到的基因数量。nCount_RNA是细胞内检测到的分子总数。nFeature_RNA过低，表示该细胞可能已死/将死或是空液滴。太高的nCount_RNA和/或nFeature_RNA表明“细胞”实际上可以是两个或多个细胞。结合线粒体基因count数除去异常值，即可除去大多数双峰/死细胞/空液滴，因此它们过滤是常见的预处理步骤。 参考自：https://www.biostars.org/p/407036/
3.2计算线粒体基因比例
人类的线粒体基因都是以MT-开头
seu.obj[["percent.mt"]] <- PercentageFeatureSet(seu.obj, pattern = "^MT-")
head(seu.obj@meta.data)
## orig.ident nCount_RNA nFeature_RNA percent.mt
## AAACCCACAGTCGGTC-1 SeuratProject 4243 1256 6.292717
## AAACGAAAGAATCGCG-1 SeuratProject 7307 2577 2.572875
## AAAGAACAGCTTACGT-1 SeuratProject 8154 1881 4.083885
## AAAGAACAGGTTCATC-1 SeuratProject 8223 2182 4.377964
## AAAGAACAGTCTGTAC-1 SeuratProject 3884 1377 4.763131
## AAAGAACTCCACCTCA-1 SeuratProject 3997 1307 3.402552
VlnPlot(seu.obj,
features = c("nFeature_RNA",
"nCount_RNA",
"percent.mt"),
ncol = 3,pt.size = 0.5)
小提琴的宽窄代表对应纵坐标的细胞数量多少。如果以后遇到点的数量太多，把小提琴都给盖上的情况，那就pt.size = 0，大小为0就是不画点。
根据上面的小提琴图里可以看到所有细胞的三个指标的分布情况，我们过滤就是把三个指标比大部队数值大的点给过滤掉。例如这个数据过滤标准可以是：
seu.obj = subset(seu.obj,
nFeature_RNA < 6000 &
nCount_RNA < 30000 &
percent.mt < 18)
标准不是唯一的，大点儿小点儿问题不大。如果没有尖尖的线，不过滤也可以。有的公共数据拿到时就已经是过滤好的。尖尖部分对应的纵坐标数值太大，属于离群值，就不要了。 一般过滤掉尖尖部分就可以了，不要过滤的太狠，标准见下：
过滤后
4.降维聚类分群
4.1 理解降维这件事
对于表达矩阵来说，一个基因就是一个维度，有几万个维度，太多了，需要减少，称之为降维。
PCA是主成分分析，完成线性的降维，会把几万个维度转换为几十个新的维度，而UMAP和tSNE是非线性的降维，线性降维不够彻底，整点高级的，进一步降维。
UMAP 和tSNE二选一就行，没有必须用哪个的说法，只能说UMAP用的多一些，图一般会紧凑一些。
4.2 file.exists——存在即跳过
file.exists("filename")
#存在的文件名会返回TRUE，而不存在与工作目录下的文件名则会返回FALSE
因此下面的代码是，结合if语句实现了分情况讨论：
如果工作目录下不存在”obj.Rdata”这个文件，则运行大括号里的代码(最后一句是save，所以运行完也就保存了”obj.Rdata”)； 如果工作目录下存在”obj.Rdata”这个文件(说明大括号里的代码已经运行过了)，则不运行大括号里的代码，也就跳过了这个耗时比较长的限速步骤
总之这个技巧可以避免多次重复运行限速步骤，非常实用。但需要注意⚠，一旦输入数据有改动，比如说前面的过滤标准改了，那么就要删除工作目录下的obj.Rdata文件，才能重新运行这段代码。
f = "obj.Rdata"
if(!file.exists(f)){
seu.obj = seu.obj %>%
NormalizeData() %>%
FindVariableFeatures() %>%
ScaleData(features = rownames(.)) %>%
RunPCA(features = VariableFeatures(.)) %>%
FindNeighbors(dims = 1:15) %>% #选择前15个维度，15是个比较折中的值，选择的数量越多计算量越大，而且信息越冗余；选的太少又不具有代表性。一般是10-20
FindClusters(resolution = 0.5) %>% #分群的分辨率是0.5，这个参数的选择范围是0.1-1.5之间，数值越大分出来的群越多，数值越小分出来的点越少
RunUMAP(dims = 1:15) %>%
RunTSNE(dims = 1:15)
save(seu.obj,file = f)
}
load(f)
ElbowPlot(seu.obj)
p1 <- DimPlot(seu.obj, reduction = "umap",label = T)+NoLegend();p1
标记#的两个地方可能会需要改动
一个是dims = 1:15,代表选择前15个维度，15是个比较折中的值，选择的数量越多计算量越大，而且信息越冗余；选的太少又不具有代表性。一般是10-20，根据ElbowPlot来选择拐点的横坐标(拐点就是在这个点之前纵坐标下降比较快，在这个点之后纵坐标下降比较慢)。不是很重要，只要拐点在15之前，选15就一点问题都没有。
一个是resolution = 0.5，代表分群的分辨率是0.5，这个参数的选择范围是0.1-1.5之间，数值越大分出来的群越多，数值越小分出来的点越少。0.5也是一个比较折中的值，如果后面注释不困难，就不用改动。
如果决定要改动那么就不能只改代码，要把obj.Rdata从工作目录下删掉，再重新运行修改后的代码。
