-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathscRNA_basic_from_huage
54 lines (52 loc) · 3.17 KB
/
scRNA_basic_from_huage
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
#project见 huage_prac
Sys.setenv(LANGUAGE='en')
options(strinngAsFactors = FALSE)
rm(list = ls())
setwd()
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)
library(ggplot2)
#读取某路径下的所有文件:利用Read10X函数的特性(目录下所有的,并且根据vector的名字直接命名barcode)
dir = c("young_1/","old_1")
names(dir) = c("young","old")
dir
scRNA.all = Read10X(data.dir = dir)
scRNA = CreateSeuratObject(scRNA.all,min.cells = 3, min.features = 40)
#table查看分类,交替使用@$符号进行对V4对象的提取,差异分析中也用到deq?
table([email protected]$orig.ident)
#质控包括去除低质量细胞和双细胞,前者包括总测序深度、线粒体基因比例等
?PercentageFeatureSet
#正则表达式,grep函数的使用,线粒体基因命名的规范 MT-
scRNA[["precent.mt"]] = PercentageFeatureSet(scRNA, pattern = "^MT-")
#红细胞相关基因 固定的? 基因集 约为10个? 一般细胞数量少,可以不做
#Feature 每个细胞中检测到的基因数,不可过少,总共2w+个基因
#Count 每个细胞总的UMI数,可以理解为测序深度。类似所有基因的reads?
#根据小提琴图,确定阈值选择,一般线粒体基因<20%,nFeature_RNA > 300 & nFeature_RNA < 7000,
scRNA_QC = subset(scRNA, subset = nFeature_RNA > 300 & nFeature_RNA < 7000 & precent.mt <20 & nCount_RNA >1000)
#细胞数量得到下降
?NormalizeData
#表达量x10000/total UMI +1 后 取ln、.目的是将测序深度归为1
scRNA_QC = NormalizeData(scRNA_QC, normalization.method = "LogNormalize", scale.factor = 10000 )
#归一化后的数据储存在assays-RNA-layers中counts下面的data
#降维: 每个gene是一个维度。考虑的是变化大的gene,去繁从简
#2维umap、3维umap。见于2024CNS 。分别代表所有维度的降维后的特征1、2、3
#降维聚类的基本pipeline:找高变基因、ScaleData、RunPCA、FindNeighbors、FindClusters
#高变gene储存在assays-RNA-layers-meta.data的var.features中
?FindVariableFeatures
scRNA_QC = FindVariableFeatures(scRNA_QC, selection.method = "vst", nfeatures = 3000 )
#数据放缩 scale,目的是标准化数据,变成均值为0,方差为1,以将所有基因的权重变为相同重要的(避免管家基因的影响等)
#默认是对高变基因进行放缩,储存在data下的scaledata中。
scRNA_QC = ScaleData(scRNA_QC , verbose = TRUE)
#PCA主成分分析降维,默认是50,储存在reduction中
scRNA_QC = RunPCA(scRNA_QC , features = VariableFeatures(scRNA_QC) , verbose = TRUE, npcs = 30)
#展示每个pca维度对整体的贡献度。越往前的pc包含的特征信息越多。越往后的PC越可能有更多的技术噪音,前面更可能反映生物信息
#三维PCA图显示前三个pc的占比,及各样本的空间分布情况。
ElbowPlot(scRNA_QC , ndims = 30, reduction = "pca" )
#选择合适的PC数进行后续的降维
scRNA_QC = RunUMAP(scRNA_QC , reduction = "pca", dims = 1:20)
#降维后进行聚类。聚类:高维的无向有权图,邻接矩阵
scRNA_QC = FindNeighbors(scRNA_QC , reduction = "pca", dims = 1:20)
scRNA_QC = FindClusters(scRNA_QC , resolution = 0.5)
DimPlot(scRNA_QC ,reduction = "umap", label = T )