单细胞第一节课

数据来源

##系统报错改为英文
Sys.setenv(LANGUAGE = "en")
##禁止转化为因子
options(stringsAsFactors = FALSE)
##清空环境
rm(list=ls())
###4.0 安装包
install.packages('Seurat') #我采用的本地安装
###加载所需要的包
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)
##读取10x的数据
scRNA.counts=Read10X("D:/genetic_r/study/geo/study7/GSE152048_BC10.matrix/BC10/")
class(scRNA.counts)
#[1] "dgCMatrix"
#attr(,"package")
#[1] "Matrix"
###创建Seurat对象
scRNA = CreateSeuratObject(scRNA.counts ,min.cells = 3,project="os", min.features = 300)
#mincell 每个基因至少要在3个细胞中表达
# min.features 每个细胞至少有多少基因表达 一个细胞中基因最小表达量
#查看样本的细胞数量
table(scRNA@meta.data$orig.ident) #os 17481     
##计算质控指标
#过滤线粒体基因 
#计算细胞中线粒体基因比例
?PercentageFeatureSet
y=scRNA[[]]
scRNA[["percent.mt"]] <- PercentageFeatureSet(scRNA, pattern = "^MT-")
x <-PercentageFeatureSet(scRNA, pattern = "^MT-")
#计算红细胞比例
HB.genes <- c("HBA1","HBA2","HBB","HBD","HBE1","HBG1","HBG2","HBM","HBQ1","HBZ")
HB_m <- match(HB.genes, rownames(scRNA@assays$RNA)) 
HB.genes <- rownames(scRNA@assays$RNA)[HB_m] 
HB.genes <- HB.genes[!is.na(HB.genes)] 
scRNA[["percent.HB"]]<-PercentageFeatureSet(scRNA, features=HB.genes) 
#head(scRNA@meta.data)
col.num <- length(levels(scRNA@active.ident))
####Feature、count、线粒体基因、红细胞基因占比可视化。
violin <- VlnPlot(scRNA,
                  features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.HB"), 
                  cols =rainbow(col.num), 
                  pt.size = 0.01, #不需要显示点，可以设置pt.size = 0
                  ncol = 4) + 
  theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank()) 
###把图片画到画板上面
violin

plot_zoom

换一个小数据练习

##系统报错改为英文
Sys.setenv(LANGUAGE = "en")
##禁止转化为因子
options(stringsAsFactors = FALSE)
##清空环境
rm(list=ls())
###4.0 安装包
install.packages('Seurat')
getwd()
###加载所需要的包
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)
####清除环境，就是把环境清空了
rm(list=ls())

##读取10x的数据
scRNA.counts=Read10X("D:/genetic_r/study/geo/lesson7/GSE152048_BC2.matrix/BC2/")

class(scRNA.counts)
"dgCMatrix"
attr(,"package")
[1] "Matrix"  稀疏矩阵

m1 <- scRNA.counts[1:20,1:20]
m2 <- as.matrix(m1)
###创建Seurat对象
?CreateSeuratObject
scRNA = CreateSeuratObject(scRNA.counts ,min.cells = 3,project="os", min.features = 300)
#mincell 每个基因至少要在3个细胞中表达
# min.features 每个细胞至少有多少基因表达 一个细胞中基因最小表达量
view(scRNA)
View(scRNA)



### s3 data.frame [,] list[[]]  character[] matirx 
#两个提取原则 1 排列对象 2 名称
###s4 层级结构

##第一种s4对象的提取方法  点击白框
phe=scRNA@meta.data
count=scRNA@assays[["RNA"]]@counts
z=scRNA@assays[["RNA"]]@counts@Dimnames[[2]]

#####第二提取s4对象的方法 @ $交替使用
count=scRNA@assays$RNA@counts
scRNA@meta.data


#查看样本的细胞数量
table(scRNA@meta.data$orig.ident) #os 5937      
##计算质控指标
#过滤线粒体基因 
#计算细胞中线粒体基因比例
?PercentageFeatureSet
y=scRNA[[]]
scRNA[["percent.mt"]] <- PercentageFeatureSet(scRNA, pattern = "^MT-")
x <-PercentageFeatureSet(scRNA, pattern = "^MT-")
#计算红细胞比例
HB.genes <- c("HBA1","HBA2","HBB","HBD","HBE1","HBG1","HBG2","HBM","HBQ1","HBZ")
HB_m <- match(HB.genes, rownames(scRNA@assays$RNA)) 
HB.genes <- rownames(scRNA@assays$RNA)[HB_m] 
HB.genes <- HB.genes[!is.na(HB.genes)] 
scRNA[["percent.HB"]]<-PercentageFeatureSet(scRNA, features=HB.genes) 
#head(scRNA@meta.data)
head(scRNA@meta.data)
col.num <- length(levels(scRNA@active.ident))
####Feature、count、线粒体基因、红细胞基因占比可视化。
violin <- VlnPlot(scRNA,
                  features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.HB"), 
                  cols =rainbow(col.num), 
                  pt.size = 0.01, #不需要显示点，可以设置pt.size = 0
                  ncol = 4) + 
  theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank()) 

violin
###这几个指标之间的相关性。 把图画到画板上，然后手动保存
plot1=FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2=FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot3=FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "percent.HB")
pearplot <- CombinePlots(plots = list(plot1, plot2, plot3), nrow=1, legend="none") 
plot1
####看画板
plot2
####看画板
plot3
####看画板
pearplot
####看画板
#自己选择性保存图片，但是 pearplot的图片必须要看，因为这个是做质控用的
##我们可以看到，nFeature_RNA的范围在0到8000之内，percent.mt代表线粒体含量
###我们默认线粒体含量至少要小于20%，这是根据生物学知识得出的默认阈值。红细胞的数目要至少小于5%
###至于nFeature_RNA和nCount_RNA的阈值怎么确定，这个要结合 pearplot的图来判断。我们质控的目标就是删除离异值。而且注意阈值尽可能取的宽松一下，防止后面分析想要的细胞得不到。
###接下来从pearplot的图片来做质控---剔除离异值
##nFeature_RNA选择大于300 小于2000的 nFeature_RNA选择小于100000，percent.mt小于10，percent.HB小于5
scRNA1 <- subset(scRNA, subset = nFeature_RNA > 300& nFeature_RNA < 2000 & percent.mt < 10 & percent.HB < 3 & nCount_RNA < 10000)
scRNA
scRNA1
violin1 <- VlnPlot(scRNA1,
                  features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.HB"), 
                  cols =rainbow(col.num), 
                  pt.size = 0.01, #不需要显示点，可以设置pt.size = 0
                  ncol = 4) + 
  theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank()) 

violin1
####过滤完之后 我们就要对数据进行均一化，使用NormalizeData这个函数。
###注意均一化是用NormalizeData，标准化是用ScaleData
?NormalizeData
scRNA11 <- NormalizeData(scRNA1, normalization.method = "LogNormalize", scale.factor = 10000)
###好了，这一节数据加载、质控的内容就算是做完了。
###在我们关闭rstudio之前 先把环境中运行好的数据保存一下
###数据将保存在之前设定好的路径中。还有保存的scRNA1，不是scRNA，因为scRNA1才是过滤好的数据。
save(scRNA11,file='scRNA11.Rdata')