单细胞分析-seurat包

R Markdown

安装环境

导入数据

setwd("D:/scRNA-seq")

pbmc.data <- Read10X(data.dir = 'hg19')#导入10x矩阵

#对数据进行初步赛选

pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
#数据导入，筛选，将一个基因必须在三个细胞内表达，一个细胞必须要有200个基因片段。
pbmc  #自己看看有多少个基因和细胞是有效的被耐人到pbmc中。

##计算线粒体基因表达情况–质控–计算线粒体比例，不好的时候线粒体可能表达偏高的，’^MT-’–线粒体基因名字特征
###过滤一些低质量的细胞同时查看赛选后的样本分布情况展示每个群中细胞的基因数和UMI数根据分群结果，展示每个群中细胞数和基因数的分布,基因数过双细胞现象，线粒体基因过高数据也不行，去掉过高和过低的值
###质控的指标：UMI、基因数、线粒体MUI的比例 低质量的细胞其基因数很低、存在线粒体污染；双细胞的基因数很高； 比较简单的方法：直接去掉离群值（过高和过低）

pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA"), ncol = 3,pt.size =0)
##pt.size =0时点消失

根据数据特征过滤

#计算四分位点，得到四分之三和四分之一的值，进行过滤差异值为四分之三处减去四分之一处，然后去掉两个分支外1.5IQR以外的值
Q1 <- quantile(pbmc$nFeature_RNA)[2] 
Q3 <- quantile(pbmc$nFeature_RNA)[4]
IQR <- Q3 - Q1
upper <- Q3+1.5*IQR
lower <- Q1-1.5*IQR
pbmc <- subset(pbmc, subset = nFeature_RNA > lower & nFeature_RNA < upper & percent.mt < 5)
print(quantile(pbmc$nFeature_RNA))
pbmc

#数据的标准化
##标准化之后的数据存放在pbmc[[“RNA”]]@data
##不确定参数设置时，可以选择默认参数：pbmc <- NormalizeData(object = pbmc)

pbmc <- NormalizeData(object = pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc[["RNA"]]@data[1:3,1:3]

#识别高度可变的特征（特征选择）—矩阵中存在这个高度可变的特征能够代表这个一种细胞

pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
#默认情况会选择前2000个基因，nfeatures可以设置选择的基因数
top10 <- head(VariableFeatures(pbmc), 10)
#前十的选出来
top10

##展示选择了哪些基因
###黑色为没有选上的基因，红色位选上的基因，命名为高表达基因top10

# plot variable features with and without labels

plot1 <- VariableFeaturePlot(pbmc)

plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE) 
CombinePlots(plots = list(plot1, plot2))

#对数据进行缩放

all.genes <- rownames(pbmc)#将所有基因赋值到变量all.gennes中
pbmc <- ScaleData(pbmc, features = all.genes) #数据缩放后的结果保存在pbmc[["RNA"]]@scale.data
sacle_dat=as.data.frame(pbmc[["RNA"]]@scale.data)
sacle_dat[1:3]

#数据量大的时候使用高度特征缩放，特征少就都缩放没事

线性降维—将相同的细胞类型放到一起去

pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc),verbose = FALSE) # 通过设置 features来选择用于降维的基因，这里选择的时高度特征基因
DimPlot(object = pbmc, reduction = "pca")
#PCA分组
ElbowPlot(pbmc)
#选择一些拐点处的pca进行分析，变速变小或不变说明后面的pca影响越小

细胞聚类

pbmc <- FindNeighbors(pbmc, dims = 1:10)#根据影响选择1：10个维度分组
pbmc <- FindClusters(pbmc, resolution = 0.7)  #resolution设置分辨率，即分群的数目，需不断调试；同样的数据resolution越大，分的群越多

#TSUN/UMAP的可视化

pbmc <- RunTSNE(pbmc, dims = 1:10)
# If you haven't installed UMAP, you can do so via reticulate::py_install(packages = 'umap-learn')
pbmc <- RunUMAP(object = pbmc, dims = 1:10)
#reduction 设置进行可视化的方法 TSNE/UMAP
DimPlot(pbmc, reduction = "tsne") 
DimPlot(object = pbmc, label = T,label.size = 5,reduction = "umap") 

寻找差异表达的特征–找差异性基因

# 寻找某个cluster1相对于cluster2和3的marker gene; ident.2默认是除了ident1以外的所有cluster
markers1 <- FindMarkers(object = pbmc, ident.1 = 1, ident.2 = c(2,3), min.pct = 0.25)
# 计算每个cluster相对于其他cluster 的marker gene 
pbmc.markers <- FindAllMarkers(object = pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25) #logfc.threshold设置的差异倍数
head(pbmc.markers)
write.table(pbmc.markers,"Marker_list.txt",row.names = F,quote = F,sep = "\t")

可视化

VlnPlot(object = pbmc, features =c("MS4A1", "CD79A"),pt.size = 0)

###上图发现差异很小可能由于分组不一样所以导致结果不一样，两个中表达量差不多

##CD79A的TSNE图谱发现大部分

FeaturePlot(object = pbmc,features = c("CD79A"),cols = c("#C0C0C0", "#E41C12"),reduction = "tsne",label = T)
print('通过图谱显示可以发现CD79A基因基本出现在3好区域')

#更换颜色

FeaturePlot(object = pbmc,features = c("FCGR3A"),cols =c('#1fca04','#cf5a34','#adda5f'),reduction = "tsne",label = T)
print('通过图谱显示可以发现CD79A基因基本出现在3好区域')

#当鉴定完细胞类型以后，需要将细胞类型的名称重新映射到分群的图上。尝试替换cluster的名称，可自行定义。

new.cluster.ids=top10[1:9]
names(new.cluster.ids)=levels(pbmc)
pbmc<-RenameIdents(pbmc,new.cluster.ids)
DimPlot(object = pbmc, label = T,label.size = 5,reduction = "umap",pt.size = 0.5)