以PBMC数据集进行探索
1 创建对象
library(Seurat)
pbmc.data <- Read10X(data.dir ="./")
class(pbmc.data)
dim(pbmc.data)
2 创建对象
pbmc <- CreateSeuratObject(counts = pbmc.data,
project = "pbmc3k", min.cells = 3, min.features = 200)
Seurat.PNG
Seurat通过CreateSeuratObject函数创建对象后,将我们导入的UMI count原始稀疏矩阵储存在pbmc@assays[["RNA"]]@counts,此外Seurat自动计算每个细胞总的UMI count,即每一列数字之和,储存在pbmc@meta.data[["nCount_RNA"]];计算每个细胞总的基因数,每一列非0的行数,储存在pbmc@meta.data[["nFeature_RNA"]]
我们也可以自己计算验证:
2.1 计算nCount_RNA
2.1.1 手动计算每一个细胞的UMI count数
matrix<-as.matrix(pbmc@assays[["RNA"]]@counts)
a<-apply(matrix,2,sum) #2对列操作,1对行操作,sum求和
head(a)
#AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1
# 2419 4903 3147
#AAACCGTGCTTCCG-1 AAACCGTGTATGCG-1 AAACGCACTGGTAC-1
# 2639 980 2163
2.1.2 与Seurat自动计算的结果作对比
head(pbmc@meta.data[["nCount_RNA"]])
#[1] 2419 4903 3147 2639 980 2163
head(pbmc[["nCount_RNA"]])
# nCount_RNA
#AAACATACAACCAC-1 2419
#AAACATTGAGCTAC-1 4903
#AAACATTGATCAGC-1 3147
#AAACCGTGCTTCCG-1 2639
#AAACCGTGTATGCG-1 980
#AAACGCACTGGTAC-1 2163
pbmc@meta.data[["nCount_RNA"]]结果为向量,pbmc[["nCount_RNA"]]结果为矩阵
2.2 计算nFeature_RNA
手动计算前六个细胞的基因Feature_RNA数
#提取并将稀疏矩阵转换为普通矩阵
b=as.matrix(pbmc@assays[["RNA"]]@counts)
e=NULL
for (i in 1:6) {c=nrow(subset(b,b[,i] != 0));e=append(e,c)}
print(e)
#[1] 779 1352 1129 960 521 781
与Seurat自动计算的结果作对比
head(pbmc@meta.data[["nFeature_RNA"]])
#[1] 779 1352 1129 960 521 781
2.2 计算线粒体基因UMI count百分比
手动计算
matrix<-as.matrix(pbmc@assays[["RNA"]]@counts)
a<-apply(matrix,2,sum)
b<-matrix[grep("MT-", row.names(matrix),ignore.case = T),]
l=apply(b,2,sum)
e=NULL
for (i in 1:6) {f=l[i]*100/a[i];e=append(e,f);}
print(e)
#AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1
# 3.0177759 3.7935958 0.8897363
#AAACCGTGCTTCCG-1 AAACCGTGTATGCG-1 AAACGCACTGGTAC-1
# 1.7430845 1.2244898 1.6643551
与Seurat自动计算的结果作对比
pbmc[["percent.mt"]] <- PercentageFeatureSet(object = pbmc, pattern = "^MT-")
head(pbmc@meta.data[["percent.mt"]])
#[1] 3.0177759 3.7935958 0.8897363 1.7430845 1.2244898
#[6] 1.6643551
捕获.PNG
2.3 标准化
手动计算第一列(细胞)标准化后的值
g<-as.matrix(pbmc@assays[["RNA"]]@counts)[,1]
a<-log(1+10000g/sum(g))
min(a);max(a);mean(a)
#[1] 0
#[1] 5.753142
#[1] 0.1170937
Seurat自动计算
pbmc <- NormalizeData(pbmc,
normalization.method = "LogNormalize",
scale.factor = 10000)
h<-pbmc@assays[["RNA"]]@data[,1]
min(h);max(h);mean(h)
#[1] 0
#[1] 5.753142
#[1] 0.1170937
标准化类似于RNA-seq的FPKM,计算公式为:In( 1 + ( UMIA ÷ UMITotal ) × 10000 ),R语言中log表示自然对数In(以e为底的对数),此处的 UMITotal即为每个细胞对应的nCount_RNA。
pbmc@assays[["RNA"]]@data在未标准化之前,储存的是原始稀疏矩阵,标准化后,储存的为标准化的矩阵。
2.3 归一化
目前还不知道手动计算公式。
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
捕获.PNG
归一化后的矩阵储存在pbmc@assays[["RNA"]]@scale.data
归一化与标准化区别:
- 标准化是对每一列(细胞)操作,归一行对行(基因)操作
- 归一化后,每一行有一半数值变为负值,标准化后都是≥0的值
网友评论