#对neocortex(新皮质)样本的处理
#加载包
library(Seurat)
library(dplyr)
library(Matrix)
library(RANN)
library(igraph)
library(matrixStats)
library(qlcMatrix)
source("doFastPCA.R")
#加载文件load(file="/kriegsteinlab/data1/aparna/homefiles/cleanedobjects_bysample_June2019/neocortex.RData")
#NormalizeDate()消除不同细胞测序深度的影响
library(Seurat)
NormalizeData(object, assay=NULL, normalization.method = "LogNormalize", scale.factor = 10000)
函数默认将每个细胞的文库大小设置成为10000个reads大小
normalization.method = "LogNormalize"参数:
单细胞数据当中有很多基因的reads数很多,甚至上千,但是有很多基因却是个位数甚至0,那这种数据离散程度也是很大的,但是我们会发现当我们对1000取以10为底的对数时,就变成了3,对10取以10为底的对数时,就变成了1,这样就实现了降低数据离散程度的目的。但是这也存在一定的问题,如果一个基因的reads数为0,那岂不是不能取对数?我们对所有的值都加上1不就能够解决这个问题?
Neocortex <- NormalizeData(object = Neocortex, normalization.method = "LogNormalize")
#FindVariableFeatures()挑选高变基因,使用conuts不是data
#nfratures=2000:2000个高编辑与储存在var.features
Neocortex <- FindVariableGenes(object = Neocortex, mean.function = ExpMean,dispersion.function = LogVMR, x.low.cutoff = 0.0125, x.high.cutoff = 3, y.cutoff = 0.5)
#读取ribogene.txt文件
ribo <- read.table("/kriegsteinlab/data1/aparna/homefiles/ribogenes.txt", sep=" ", row.names=1)
#rownames函数用来拆看或更改“行” 函数.
#intersect()用来查找两个对象的交集
ribo.genes <- intersect(rownames(Neocortex@data),rownames(ribo))
#colSum()用于计算矩阵或数组列的综合
percent.ribo<-colSums(expm1(Neocortex@data[ribo.genes,]))/colSums(expm1(Neocortex@data))
#给seurat对象添加metadata用AddMetaData函数,添加ribo基因百分比,第三个参数是col.name,列名
Neocortex <- AddMetaData(Neocortex, percent.ribo, "percent.ribo")
#单细胞基因表达counts矩阵数据经过NormalizeData()归一化处理后,还需要进行scale标准化。ScaleData()函数将归一化的基因表达转换为Z分数(值以 0 为中心,方差为 1)。 它存储在 seurat_obj[['RNA']]@scale.data,用于下游的PCA降维。 默认是仅在高可变基因上运行标准化。
- 使每个基因在所有细胞的表达量均值为 0
- 使每个基因在所有细胞中的表达量方差为 1
#vars.to.regress="Individual"表示要回归的变量是Individual
Neocortex <- ScaleData(object = Neocortex, genes.use = Neocortex@var.genes, display.progress = FALSE, do.center=TRUE, vars.to.regress="Individual")
#
matrix_scaled <- Neocortex@scale.data
matrix_pca <- doFastPCA(t(matrix_scaled),50,center.use = F, scale.use = F,iterative = F)
ev <- matrix_pca$sdev^2
max(which(matrix_pca$sdev^2>(sqrt(length(row.names(Neocortex@data))/length(colnames(Neocortex@data))) + 1)^2))
FindAllMarkers()比较一个cluster和其他cluster差异基因的表达