注意此笔记内容总结于生信星球公众号内容
要注意一点,每个样本细胞数量不会超过一万,因此必须选用16g以上内存的电脑,庆幸换了电脑,yeah
直接在浏览器搜索GEO号
前三个文件属于一个样本,因为编号都是一致的
之后页面拉到最后直接选择Download,这时会下载一个压缩包,得需要解压之后才可以用
> #查看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"
对文件是有要求的,必须固定而且不能有前缀
“barcodes.tsv.gz”:存储的是barcodes,相当于细胞的编号,是表达矩阵的列名。 “features.tsv.gz”:存储的是基因名称,是表达矩阵的行名。 “matrix.mtx.gz”:存储的是每个位置的数值,是表达矩阵的内容,仅存储了非零的数值。
这样就会产生一个疑问,那如果是多样本咋办,多个样本就需要每个样本都有一个文件夹,文件夹里面的命名也必须和这个保持一致
> 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"
读取文件用的是Read10X这个函数,接受的参数是文件夹名称,文件夹里面三个数据合并在一起才是完整的单细胞表达矩阵,每个都只是存储了一部分
> rm(list = ls())
> #清空右上角的所有变量,方便反复调试代码
> library(Seurat)
> library(patchwork)
> library(tidyverse)
> ct = Read10X("input/") #重新指定到别的变量,方便修改且不影响原来的
> dim(ct)
[1] 33538 10218
> dim(ct) #返回的数字分别是行数和列数
[1] 33538 10218
> 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值比较多。
> seu.obj <- CreateSeuratObject(counts = ct,
+ min.cells = 3, #一个基因至少要在3个细胞里有表达,才被保留
+ min.features = 200) #一个细胞里至少要200个基因有表达,才被保留
> dim(seu.obj) #检查行列数
[1] 20648 10105
行列数都变少了,和花花老师的可以对上
为了节省计算资源
> set.seed(1234) #set.seed(1234)是设计随机种子,作用是让后面的随机抽样变得可重复,只要多次运行时,setseed的括号里的数字没变,抽样的结果就不会变。
> seu.obj = subset(seu.obj,downsample = 3000)
> dim(seu.obj) #检查行列数
[1] 20648 3000
行数没变,抽取了3000列
广义的“对象”是Rstudio右上角所有的数据,向量数据框矩阵列表都是对象。
狭义的“对象”是由R包的作者定义的,以固定模式组织的数据,里面的结构都是规定好的。
> class(seu.obj) #Seurat对象,出自于SeuratObject这个包
[1] "Seurat"
attr(,"package")
[1] "SeuratObject"
> class(ct)
[1] "dgCMatrix"
attr(,"package")
[1] "Matrix"
根据上面的那个,下面这个很好理解
要提取它的子集,有@和$两个符号,具体可以尝试一下看看
> 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(因为有三个基因的表达量不为零)。
%>%是向后传递,类似于linux里面的 |
线粒体基因含量不能过高;线粒体基因的表达水平在单细胞RNA测序中是一个重要的质控指标。通常,线粒体基因的表达量占总表达量的百分比用于衡量细胞的质量。高比例的线粒体基因表达通常表明细胞质量较差。在线粒体基因表达高的情况下,可能是因为细胞已经进入了凋亡或坏死过程。
线粒体基因表达占总表达量的比例超过20%或30%的细胞会被认为是低质量细胞
nFeature_RNA 不能过高或过低:nFeature_RNA过低,说明该细胞中只有少量基因有表达。这可能表示该细胞质量较差。如果nFeature_RNA过高,可能是由于双细胞。
人类的线粒体基因都是以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,结果他的是心肌细胞,线粒体含量本来就高,这样一过滤一万细胞剩了几百,离大谱啦!注意一下
降维是通过数学或统计方法将高维数据映射到低维空间的过程,目的是保留数据中的主要变化和模式,同时减少数据的复杂性和计算负担。这使得数据能够更容易地可视化、分析和理解,帮助发现隐藏在数据中的结构和关系。例如,如果有数以万计的基因,每个细胞的表达数据就是一个数以万计的维度。这种高维度的数据不仅难以可视化,也增加了计算复杂度和存储需求。降维可以将数据投影到更低维度的空间中,保留尽可能多的数据变异性和信息,同时减少数据的复杂度和计算负担。
主成分分析(PCA):通过线性变换将数据映射到低维空间,保留数据中方差最大的方向。
t分布邻域嵌入(t-SNE):非线性降维方法,能够保留高维空间中数据点之间的局部相似性关系。
UMAP(Uniform Manifold Approximation and Projection):一种高效的非线性降维技术,能够更好地保留数据的全局和局部结构。
一般先用PCA进行线性降维,线性降维不够彻底,UMAP和tSNE是非线性的降维,所以要进一步利用UMAP和tSNE降维,UMAP用的更多一点
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 删除。