代码学习neocortex/R_general/neocortex.R

#对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差异基因的表达