单细胞day2
原创
注意此笔记内容总结于生信星球公众号内容
要注意一点，每个样本细胞数量不会超过一万，因此必须选用16g以上内存的电脑，庆幸换了电脑，yeah
1下载和整理数据
直接在浏览器搜索GEO号
前三个文件属于一个样本，因为编号都是一致的
之后页面拉到最后直接选择Download，这时会下载一个压缩包，得需要解压之后才可以用
1.1解压缩
> #查看Seurat的版本
> packageVersion("Seurat")
[1] '5.1.0'
> #解压文件
> untar("GSE231920_RAW.tar",exdir = "input")
> dir("input")
[1] "GSM7306054_sample1_barcodes.tsv.gz" "GSM7306054_sample1_features.tsv.gz"
[3] "GSM7306054_sample1_matrix.mtx.gz"
> dir("input") #罗列一下input文件夹里面的内容
[1] "GSM7306054_sample1_barcodes.tsv.gz" "GSM7306054_sample1_features.tsv.gz"
[3] "GSM7306054_sample1_matrix.mtx.gz" 1.2单细胞文件组织的要求
对文件是有要求的，必须固定而且不能有前缀
“barcodes.tsv.gz”：存储的是barcodes，相当于细胞的编号，是表达矩阵的列名。 “features.tsv.gz”：存储的是基因名称，是表达矩阵的行名。 “matrix.mtx.gz”:存储的是每个位置的数值，是表达矩阵的内容，仅存储了非零的数值。
这样就会产生一个疑问，那如果是多样本咋办，多个样本就需要每个样本都有一个文件夹，文件夹里面的命名也必须和这个保持一致
1.3 改名
> library(stringr)
> x = str_remove(dir("input/"),"GSM7306054_sample1_")
> x
[1] "barcodes.tsv.gz" "features.tsv.gz" "matrix.mtx.gz"
> xx = str_remove(dir("input/"),"GSM7306054_sample1_")
> xx
[1] "barcodes.tsv.gz" "features.tsv.gz" "matrix.mtx.gz"
> #这个R代码片段的功能是将“input”目录中的所有文件重命名。具体地，它将“input”目录下的每个文件重命名为新名称xx中指定的名称。
> file.rename(paste0("input/",dir("input/")),
+ paste0("input/",xx))
[1] TRUE TRUE TRUE
>
> dir("input/") #检查一下改名是否成功
[1] "barcodes.tsv.gz" "features.tsv.gz" "matrix.mtx.gz" 2读取并且创建Seurat对象
2.1读取文件
读取文件用的是Read10X这个函数，接受的参数是文件夹名称，文件夹里面三个数据合并在一起才是完整的单细胞表达矩阵，每个都只是存储了一部分
> rm(list = ls())
> #清空右上角的所有变量，方便反复调试代码
> library(Seurat)
> library(patchwork)
> library(tidyverse)
> ct = Read10X("input/") #重新指定到别的变量，方便修改且不影响原来的
> dim(ct)
[1] 33538 10218
> dim(ct) #返回的数字分别是行数和列数
[1] 33538 102182.2 R语言补充知识，稀疏矩阵
> class(ct) #看一下数据类型，发现返回了好多
[1] "dgCMatrix"
attr(,"package")
[1] "Matrix"
> ct[c("CD3D","TCL1A","MS4A1"), 1:30] # 提取"CD3D","TCL1A","MS4A1"这几行，的1-30列
3 x 30 sparse Matrix of class "dgCMatrix"
[[ suppressing 30 column names 'AAACCCACACAAATAG-1', 'AAACCCACAGGCTCTG-1', 'AAACCCACAGGTCCCA-1' ... ]]
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)
> dim(seu.obj) #检查行列数
[1] 20648 3000行数没变，抽取了3000列
2.5 R语言补充知识，对象
广义的“对象”是Rstudio右上角所有的数据，向量数据框矩阵列表都是对象。
狭义的“对象”是由R包的作者定义的，以固定模式组织的数据，里面的结构都是规定好的。
> class(seu.obj) #Seurat对象，出自于SeuratObject这个包
[1] "Seurat"
attr(,"package")
[1] "SeuratObject"
> class(ct)
[1] "dgCMatrix"
attr(,"package")
[1] "Matrix"根据上面的那个，下面这个很好理解
要提取它的子集，有@和$两个符号,具体可以尝试一下看看
2.6 探索Seurat对象的meta信息
> seu.obj@meta.data %>% head() #提取meta信息
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创建Seurat对象时，自动生成的，共有3列，分别是orig.ident,nCount_RNA,nFeature_RNA。
行名（barcodes）：行名是细胞的唯一标识符（细胞条形码），对应于表达矩阵中的列名。
例如，AAACCCACACAAATAG-1是一个细胞的条形码。
orig.ident表示细胞的原始分类，通常用于表示样本来源。如果在创建Seurat对象时指定了project参数，这一列会被赋值为project参数的值。如果没有指定project参数，默认值是SeuratProject。
例如，如果project = "Sample1"，那么orig.ident列的值会全部是Sample1。
nCount_RNA表示每个细胞中所有基因的表达量之和，即表达矩阵中每一列的总和。这是一个细胞中所有转录本的总数，用于衡量细胞的总转录活性。例如，如果一个细胞中有5个基因的表达量分别为1, 2, 3, 4, 5，那么该细胞的nCount_RNA值就是1+2+3+4+5 = 15。
nFeature_RNA表示每个细胞中表达量不为零的基因的数量。这是一个细胞中检测到的基因数，用于衡量细胞的基因表达复杂性。例如，如果一个细胞中有5个基因的表达量分别为1, 0, 3, 0, 5，那么该细胞的nFeature_RNA值就是3（因为有三个基因的表达量不为零）。
2.7 管道符号
%>%是向后传递，类似于linux里面的 |
3 质控
3.1 质控的指标
线粒体基因含量不能过高；线粒体基因的表达水平在单细胞RNA测序中是一个重要的质控指标。通常，线粒体基因的表达量占总表达量的百分比用于衡量细胞的质量。高比例的线粒体基因表达通常表明细胞质量较差。在线粒体基因表达高的情况下，可能是因为细胞已经进入了凋亡或坏死过程。
线粒体基因表达占总表达量的比例超过20%或30%的细胞会被认为是低质量细胞
nFeature_RNA 不能过高或过低:nFeature_RNA过低,说明该细胞中只有少量基因有表达。这可能表示该细胞质量较差。如果nFeature_RNA过高，可能是由于双细胞。
3.2 计算线粒体基因比例
人类的线粒体基因都是以MT-开头，下面的pattern = “^MT-”代表正则表达式匹配以MT-开头的基因。
> seu.obj[["percent.mt"]] <- PercentageFeatureSet(seu.obj, pattern = "^MT-")
> head(seu.obj@meta.data)
orig.ident nCount_RNA
AAACCCACAGTCGGTC-1 SeuratProject 4243
AAACGAAAGAATCGCG-1 SeuratProject 7307
AAAGAACAGCTTACGT-1 SeuratProject 8154
AAAGAACAGGTTCATC-1 SeuratProject 8223
AAAGAACAGTCTGTAC-1 SeuratProject 3884
AAAGAACTCCACCTCA-1 SeuratProject 3997
nFeature_RNA percent.mt
AAACCCACAGTCGGTC-1 1256 6.292717
AAACGAAAGAATCGCG-1 2577 2.572875
AAAGAACAGCTTACGT-1 1881 4.083885
AAAGAACAGGTTCATC-1 2182 4.377964
AAAGAACAGTCTGTAC-1 1377 4.763131
AAAGAACTCCACCTCA-1 1307 3.402552可以看到运行接收后percent.mt这个信息加到meta.data里面了
然后画一个小提琴图来制定一下筛选标准，筛掉离群点
> VlnPlot(seu.obj,
+ features = c("nFeature_RNA",
+ "nCount_RNA",
+ "percent.mt"),
+ ncol = 3,pt.size = 0.5) #设置为3，表示将三个特征的小提琴图排列在一行的三列中显示
Warning: Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.warning不用管
也可以不显示点
> VlnPlot(seu.obj,
+ features = c("nFeature_RNA",
+ "nCount_RNA",
+ "percent.mt"),
+ ncol = 3,pt.size = 0)
Warning: Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.
通过这个图确定了点聚集的位置，之后进行筛选
> seu.obj = subset(seu.obj,
+ nFeature_RNA < 6000 &
+ nCount_RNA < 30000 &
+ percent.mt < 18)
> VlnPlot(seu.obj,
+ features = c("nFeature_RNA",
+ "nCount_RNA",
+ "percent.mt"),
+ ncol = 3,pt.size = 0.2)
Warning: Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.
> 筛选完成之后在进行小提琴图的绘制，发现尖尖的线消失了
标准不是唯一的，大点儿小点儿问题不大。如果没有尖尖的线，不过滤也可以。有的公共数据拿到时就已经是过滤好的。一般不要过滤的太狠，以前有个神人，照抄代码percent.mt < 5，结果他的是心肌细胞，线粒体含量本来就高，这样一过滤一万细胞剩了几百，离大谱啦！注意一下
4降维聚类分析
4.1为什么要降维
降维是通过数学或统计方法将高维数据映射到低维空间的过程，目的是保留数据中的主要变化和模式，同时减少数据的复杂性和计算负担。这使得数据能够更容易地可视化、分析和理解，帮助发现隐藏在数据中的结构和关系。例如，如果有数以万计的基因，每个细胞的表达数据就是一个数以万计的维度。这种高维度的数据不仅难以可视化，也增加了计算复杂度和存储需求。降维可以将数据投影到更低维度的空间中，保留尽可能多的数据变异性和信息，同时减少数据的复杂度和计算负担。
4.2用什么降维
主成分分析（PCA）：通过线性变换将数据映射到低维空间，保留数据中方差最大的方向。
t分布邻域嵌入（t-SNE）：非线性降维方法，能够保留高维空间中数据点之间的局部相似性关系。
UMAP（Uniform Manifold Approximation and Projection）：一种高效的非线性降维技术，能够更好地保留数据的全局和局部结构。
一般先用PCA进行线性降维，线性降维不够彻底，UMAP和tSNE是非线性的降维，所以要进一步利用UMAP和tSNE降维，UMAP用的更多一点
4.3 R语言知识补充
file.exists是判断一个文件在工作目录下是否存在的函数。
> dir()
[1] "GSE231920_RAW.tar" "input" "install.R" "run1.R"
[5] "single cell.Rproj"
> file.exists("day5.Rproj")
[1] FALSE
> file.exists("GSE231920_RAW.tar")
[1] TRUE结合if语句实现了分情况讨论
if(!file.exists(f)){#判断是否存在，存在跳过大括号内容，不存在运行生成一个。
seu.obj = seu.obj %>%
NormalizeData() %>%
FindVariableFeatures() %>%
ScaleData(features = rownames(.)) %>%
RunPCA(features = VariableFeatures(.)) %>%
FindNeighbors(dims = 1:15) %>%
FindClusters(resolution = 0.5) %>%
RunUMAP(dims = 1:15) %>%
RunTSNE(dims = 1:15)
save(seu.obj,file = f)
}
load(f)
ElbowPlot(seu.obj)ElbowPlot(seu.obj) 绘制肘部图，其中横轴表示聚类数目（k值），纵轴表示对应的聚类性能指标（如惯性或误差平方和）。通过观察图形，可以找到一个适合的聚类数目，这个数目通常是在肘部点（拐点）附近的聚类数目。
一个是dims = 1:15,代表选择前15个维度，15是个比较折中的值，选择的数量越多计算量越大；选的太少又不能代表所有的数据。一般是10-20，根据ElbowPlot来选择拐点的横坐标(拐点就是在这个点之前纵坐标下降比较快，在这个点之后纵坐标下降比较慢)。不是很重要，只要拐点在15之前，选15就一点问题都没有。一个是resolution = 0.5，代表分群的分辨率是0.5，这个参数的选择范围是0.1-1.5之间，数值越大分出来的群越多，数值越小分出来的点越少。0.5也是一个比较折中的值，如果后面注释不困难，就不用改动
如果决定要改动那么就不能只改代码，要把obj.Rdata从工作目录下删掉，再重新运行修改后的代码。
p1 <- DimPlot(seu.obj, reduction = "umap",label = T)+NoLegend();p1
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
腾讯云开发者
