生信-单细胞数据处理

写在前面

第一次尝试用markdown来记录自己的学习经历

写作来源

https://mp.weixin.qq.com/s/8keKxz_Q_DHt9ASNKHFjtg#如何利用芯片数据,表型数据,探针数据生成单细胞对象
https://www.jianshu.com/p/5a785a7afb7d?utm_campaign=hugo#单细胞的基本分析流程

https://zhuanlan.zhihu.com/p/28844468#为什么要进行单细胞测序,以及一些简单的单细胞测序前的处理

https://www.jianshu.com/p/5a785a7afb7d?utm_campaign=hugo#单细胞后聚类的问题

https://mp.weixin.qq.com/s/fE2yPbVvD0R5I8qnNh4F1w#基本就是b站课程的总结了

https://www.jianshu.com/p/a67fb39a213a#tsne的缺点

https://mp.weixin.qq.com/s/7xPbdwvZHbJEntY_Y3q5pA 为什么说RPKM是错误的

吐槽一下jimmy的代码有时候真的是看的头疼

https://mp.weixin.qq.com/s/w0kNPq3_W33kvMHx9Nm5Pg#比较新的10x单细胞数据(Seruat的一些处理)

https://mp.weixin.qq.com/s/kp1OpZqsMb0BmNDkH2OHQQ#生信技能树的祖传的10x单细胞单样本分析

https://mp.weixin.qq.com/s/wiykclqla1Kzt7GlGSlOFQ #很多材料汇总

https://www.jianshu.com/p/284d61ba5d7c#刘小泽的代码,部分可以调用

https://mp.weixin.qq.com/s/UOSRK_pK2XFV2F41rDq0Eg# 过滤信息部分解读

https://www.jianshu.com/p/b46b6b6d344f#超多图讲解

重磅综述:三万字长文读懂单细胞RNA测序分析的最佳实践教程 (原理、代码和评述) - 知乎 https://zhuanlan.zhihu.com/p/108918012
自动化注释-singR的原理讲解:https://mp.weixin.qq.com/s?src=11&timestamp=1623321455&ver=3122&signature=ht0gOMOCKi9Turnm5FQNwEWxYiu1zZDHw1Wg*AAfS-MgJ3Q9qm6qtKUM6TsH7cHc6w4fEniO6NzjZYvKujgp7tpFJGCHPnxaBBzlqWDApDt0ymU04PeaI6LyRU8X4d0s&new=1
#生信技能树的CNS图表复现-很重要
https://mp.weixin.qq.com/s/1O1zuwLyM6_W0hZm5I26UA

开始

为什么要用单细胞测序

细胞之间存在差异,比如肿瘤的中心,周围,淋巴转移灶有不同的DNA与RNA表现与表观遗传

传统基于bulk的测序,是基于多细胞平均水平

原始数据产出形式

更多的细胞(以10X为主打)

更多的基因(以smart-seq2为主打)

单细胞基本流程

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-lEpOyMR2-1623147257930)(492DD9B55BF24E46A4741A2C0A0EE3A5)]

一些背景知识点

为什么不使用RPKM等类似的统计量

首先RPKM是基于基因长度与总测序文库大小得到的规范化的数值,而这是不合理的。没有很好的利用到转录本的概念。

其实总平均转录本丰度应该只跟基因多少有关。
所以引入了一个平均转录本长度和总转录本长度来代替文库大小的感觉

CV(标准差/均值)

很好的去除了量纲和尺度的影响

单细胞的样本配置

10X:少的话两个即可,但每个样本大概有3-8K多个细胞,每个细胞15-50K的read,1到10万的reads文库对应着200到1000个基因

Smart-seq2:每个样本细胞数量通常是96的倍数,500个左右就很厉害了,然后每个细胞的reads就可以很多,百万级别都没有问题。所以检测到基因数量也很多

质控部分

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4758103/
1.检测每个细胞中检测到的unique genes数量(这种情况,一般低质量的细胞或空的液滴"droplet"中基因数量也会较少;如果一个液滴中有两个细胞"doublets"或者存在多个细胞"multiplets",这样会导致检测到的基因数量出奇的高)

2.比对到线粒体基因组的reads数量,因为一般低质量或死亡的细胞中会广泛存在线粒体基因组污染

标化问题

问题1:为什么在NormalizeData()之后,还需要进行ScaleData操作?

前面NormalizeData进行的归一化主要考虑了测序深度/文库大小的影响(reads *10000 divide by total reads, and then log),但实际上归一化的结果还是存在基因表达量分布区间不同的问题;而ScaleData做的就是在归一化的基础上再添zero-centres (mean/sd =》 z-score),将各个表达量放在了同一个范围中,真正实现了”表达量的同一起跑线“

tsne与pca

从理论上讲,tsne更适用于单细胞数据,而pca更适用于传统rna-seq数据。为什么呢?一方面PCA能捕捉比较大的实验方差,而其实单细胞更主要的并不是大的差异,而是微小的差异。所以tsne的设计要点,类内尽可能近,类间尽可能远,就成了比较好的选择

为什么要去除核糖体和线粒体

如果线粒体RNA过高,也同样预示着细胞有破损。因为当细胞破损时,细胞质RNA会跑出来,但是线粒体RNA由于有线粒体膜的包裹,不会溢出。因此,当细胞膜有破损时,线粒体RNA所占比例会很高。注意:当细胞出现apoptosis, necrosis的时候,也会有这种现象。

核糖体RNA占比较高时,可能是因为细胞内出现了较多的RNA降解。在全长单细胞转录组中,3’偏好性可用于检测细胞内是否存在大量RNA降解

UMI,indrops,SEQC,zUMIs

单分子标签(unique molecular identifier,UMI)是人工设计的、碱基序列随机组合的一小段DNA序列,用于从浩瀚的NGS数据中分辨哪些reads来源于同一祖先分子。主要是为了提高低灵敏度

原理与流程

1.首先提取组织,裂解组织形成单个细胞

2.利用液滴或者板计术对单个细胞进行捕获

3.barcode用于区分样本,UMI用于区分序列来源扩增还是独立分子

4.对原始矩阵,对细胞进行过滤:进行空载分子、线粒体、多分子

5.对mrna进行过滤:可以以比例,也可以以个数,如果以个数的话不能大于最小细胞亚群数量,但最小细胞亚群数量不应该在后面才生成的么

6.也行还得考虑在建库前,有碎裂RNA混入样本中

7.标准化与归一化、log转换

为什么需要标准化,由于操作流程的可变性,同一个细胞不同的两种测序也可能获得不同的差异,由技术差异带来的。所以需要标准化,CPM是常用的标准化方法。除了线性标准化外,还有非线性标准化(可能更适用于板,因为不同板的批次效应更明显)。

标准化与归一化是不一样的,CPM的是标准化,z统计量是归一化。

8.数据校正和整合(整合不同于批次效应)

数据标准化是去除实验过程中随机性的影响 (count sampling)。但是,标准化后的数据可能仍然包含有与研究不相干的因素带来的影响。数据校正的目的就是进一步去除技术因素和非关注的生物学混杂因素,例如批次、dropout或细胞周期不同带来的影响 (如何火眼金睛鉴定那些单细胞转录组中的混杂因素)。需要注意的是这些混淆因素并不总是需要校正

9.缺失值填充

仅在进行轨迹推断和校正的生物学混杂因素不影响感兴趣的生物学过程时才校正这些因素的影响。

基于板的数据集预处理时需要校正count数的影响,建议采用非线性标准化方法或downsampling方法进行标准化。

10.特征选择、降维和可视化

应用方向

1.大规模细胞图谱构建

2.细胞亚群细化&稀有细胞类型鉴定

3.肿瘤方向研究

4.干细胞发育及分化研究

5.免疫方向研究

6.神经科学研究

7.发育生物学

8.生物标志物/疾病分型

R包需求等

环境的设置

# 总结一下,可以先用if判断再进行设置
if(length(getOption("CRAN"))==0) options(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")

if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)

if(length(getOption("BioC_mirror"))==0) options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")

R包的安装

# 快速安装cran包
cran_pkgs <- c("ggfortify","FactoMineR","factoextra")
for (pkg in cran_pkgs){
  if (! require(pkg,character.only=T) ) {
    install.packages(pkg,ask = F,update = F)
    require(pkg,character.only=T)
   }
}
# Bioconductor包
library(BiocManager)
bioc_pkgs <- c("scran","TxDb.Mmusculus.UCSC.mm10.knownGene","org.Mm.eg.db","genefu","org.Hs.eg.db","TxDb.Hsapiens.UCSC.hg38.knownGene")
for (pkg in bioc_pkgs){
  if (! require(pkg,character.only=T) ) {
    BiocManager::install(pkg,ask = F,update = F)
    require(pkg,character.only=T)
   }
 }

代码实操

关于单个10x的祖传代码

### ---------------
###
### Create: Jianming Zeng
### Date: 2020-08-27 15:00:26
### Email: jmzeng1314@163.com
### Blog: http://www.bio-info-trainee.com/
### Forum:  http://www.biotrainee.com/thread-1376-1-1.html
### CAFS/SUSTC/Eli Lilly/University of Macau
### Update Log: 2020-08-27  First version
###
### ---------------

rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
pro='S1'
# 搞清楚你的10x单细胞项目的cellranger输出文件夹哦
hp_sce <- CreateSeuratObject(Read10X('scRNAseq_10_s1/filtered_feature_bc_matrix/'),
                             pro)

hp_sce
rownames(hp_sce)[grepl('^mt-',rownames(hp_sce))]
rownames(hp_sce)[grepl('^Rp[sl]',rownames(hp_sce))]

hp_sce[["percent.mt"]] <- PercentageFeatureSet(hp_sce, pattern = "^mt-")
fivenum(hp_sce[["percent.mt"]][,1])
rb.genes <- rownames(hp_sce)[grep("^Rp[sl]",rownames(hp_sce))]
C<-GetAssayData(object = hp_sce, slot = "counts")
percent.ribo <- Matrix::colSums(C[rb.genes,])/Matrix::colSums(C)*100
hp_sce <- AddMetaData(hp_sce, percent.ribo, col.name = "percent.ribo")

plot1 <- FeatureScatter(hp_sce, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(hp_sce, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))

VlnPlot(hp_sce, features = c("percent.ribo", "percent.mt"), ncol = 2)
VlnPlot(hp_sce, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
VlnPlot(hp_sce, features = c("percent.ribo", "nCount_RNA"), ncol = 2)
hp_sce

hp_sce1 <- subset(hp_sce, subset = nFeature_RNA > 200 & nCount_RNA > 1000 & percent.mt < 20)
hp_sce1


sce=hp_sce1
sce
sce <- NormalizeData(sce, normalization.method =  "LogNormalize",
                     scale.factor = 10000)
GetAssay(sce,assay = "RNA")
sce <- FindVariableFeatures(sce,
                            selection.method = "vst", nfeatures = 2000)
# 步骤 ScaleData 的耗时取决于电脑系统配置(保守估计大于一分钟)
sce <- ScaleData(sce)
sce <- RunPCA(object = sce, pc.genes = VariableFeatures(sce))
DimHeatmap(sce, dims = 1:12, cells = 100, balanced = TRUE)
ElbowPlot(sce)
sce <- FindNeighbors(sce, dims = 1:15)
sce <- FindClusters(sce, resolution = 0.2)
table(sce@meta.data$RNA_snn_res.0.2)
sce <- FindClusters(sce, resolution = 0.8)
table(sce@meta.data$RNA_snn_res.0.8)


set.seed(123)
sce <- RunTSNE(object = sce, dims = 1:15, do.fast = TRUE)
DimPlot(sce,reduction = "tsne",label=T)
phe=data.frame(cell=rownames(sce@meta.data),
               res2=sce@meta.data$RNA_snn_res.0.2,
               res8=sce@meta.data$RNA_snn_res.0.8,
               cluster =sce@meta.data$seurat_clusters)
head(phe)
table(phe$cluster)
tsne_pos=Embeddings(sce,'tsne')
DimPlot(sce,reduction = "tsne",label=T,split.by ='orig.ident')

head(phe)
table(phe$cluster)
head(tsne_pos)
dat=cbind(tsne_pos,phe)
save(dat,file=paste0(pro,'_for_tSNE.pos.Rdata'))
load(file=paste0(pro,'_for_tSNE.pos.Rdata'))
library(ggplot2)
p=ggplot(dat,aes(x=tSNE_1,y=tSNE_2,color=cluster))+geom_point(size=0.95)
p=p+stat_ellipse(data=dat,aes(x=tSNE_1,y=tSNE_2,fill=cluster,color=cluster),
                 geom = "polygon",alpha=0.2,level=0.9,type="t",linetype = 2,show.legend = F)+coord_fixed()
print(p)
theme= theme(panel.grid =element_blank()) +   ## 删去网格
  theme(panel.border = element_blank(),panel.background = element_blank()) +   ## 删去外层边框
  theme(axis.line = element_line(size=1, colour = "black"))
p=p+theme+guides(colour = guide_legend(override.aes = list(size=5)))
print(p)
ggplot2::ggsave(filename = paste0(pro,'_tsne_res0.2.pdf'))

plot1 <- FeatureScatter(sce, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(sce, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))
ggplot2::ggsave(filename = paste0(pro,'_CombinePlots.pdf'))

VlnPlot(sce, features = c("percent.ribo", "percent.mt"), ncol = 2)
ggplot2::ggsave(filename = paste0(pro,'_mt-and-ribo.pdf'))
VlnPlot(sce, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
ggplot2::ggsave(filename = paste0(pro,'_counts-and-feature.pdf'))
VlnPlot(sce, features = c("percent.ribo", "nCount_RNA"), ncol = 2)


table(sce@meta.data$seurat_clusters)

for( i in unique(sce@meta.data$seurat_clusters) ){
  markers_df <- FindMarkers(object = sce, ident.1 = i, min.pct = 0.25)
  print(x = head(markers_df))
  markers_genes =  rownames(head(x = markers_df, n = 5))
  VlnPlot(object = sce, features =markers_genes,log =T )
  ggsave(filename=paste0(pro,'_VlnPlot_subcluster_',i,'_sce.markers_heatmap.pdf'))
  FeaturePlot(object = sce, features=markers_genes )
  ggsave(filename=paste0(pro,'_FeaturePlot_subcluster_',i,'_sce.markers_heatmap.pdf'))
}

sce.markers <- FindAllMarkers(object = sce, only.pos = TRUE, min.pct = 0.25,
                              thresh.use = 0.25)

DT::datatable(sce.markers)
write.csv(sce.markers,file=paste0(pro,'_sce.markers.csv'))
library(dplyr)
top10 <- sce.markers %>% group_by(cluster) %>% top_n(10, avg_logFC)
DoHeatmap(sce,top10$gene,size=3)
ggsave(filename=paste0(pro,'_sce.markers_heatmap.pdf'))
# FeaturePlot( sce,  top10$gene )
# ggsave(filename=paste0(pro,'_sce.markers_FeaturePlot.pdf'),height = 49)


library(SingleR)
sce_for_SingleR <- GetAssayData(sce, slot="data")
mouseImmu <- ImmGenData()
pred.mouseImmu <- SingleR(test = sce_for_SingleR, ref = mouseImmu, labels = mouseImmu$label.main)

mouseRNA <- MouseRNAseqData()
pred.mouseRNA <- SingleR(test = sce_for_SingleR, ref = mouseRNA, labels = mouseRNA$label.fine )

cellType=data.frame(seurat=sce@meta.data$seurat_clusters,
                    mouseImmu=pred.mouseImmu$labels,
                    mouseRNA=pred.mouseRNA$labels)
sort(table(cellType[,2]))
sort(table(cellType[,3]))
table(cellType[,2:3])
save(sce,cellType, file=paste0(pro,'_output.Rdata') )

rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
pro='S1'
library(SingleR)
load(file=paste0(pro,'_output.Rdata') )

cg=names(tail(sort(table(cellType[,2]))))
cellname=cellType[,2]
cellname[!cellname %in% cg ]='other'
table(sce@meta.data$seurat_clusters,cellname)
table(cellname)
# Epithelial cells, 0,2,5,7,8,12
# Endothelial cells, 10
# Fibroblasts , 1,3,6,9,14,15
# Macrophages , 4,13,16
# NKT , 11
# other ,4,11
# Stromal cells,
cl=sce@meta.data$seurat_clusters
ps=ifelse(cl %in% c(0,2,5,7,8,12),'epi',
       ifelse(cl %in% c(1,3,6,9,14,15),'fibro',
              ifelse(cl %in% c(4,13,16),'macro',
                     ifelse(cl %in% c(10),'endo',
                            ifelse(cl %in% c(11),'NKT','other'
                            )))))
table(ps)
cellname=paste(cl,ps)
sce@meta.data$cellname=cellname
DimPlot(sce,reduction = "tsne",label=T,group.by = 'cellname')
ggplot2::ggsave(filename = paste0(pro,'_tsne_res0.2_singleR_raw.pdf'))
dat$cluster=cellname
library(ggplot2)
p=ggplot(dat,aes(x=tSNE_1,y=tSNE_2,color=cluster))+geom_point(size=0.95)
p=p+stat_ellipse(data=dat,aes(x=tSNE_1,y=tSNE_2,fill=cluster,color=cluster),
                 geom = "polygon",alpha=0.2,level=0.9,type="t",linetype = 2,show.legend = F)+coord_fixed()
print(p)
theme= theme(panel.grid =element_blank()) +   ## 删去网格
  theme(panel.border = element_blank(),panel.background = element_blank()) +   ## 删去外层边框
  theme(axis.line = element_line(size=1, colour = "black"))
p=p+theme+guides(colour = guide_legend(override.aes = list(size=5)))
print(p)
ggplot2::ggsave(filename = paste0(pro,'_tsne_res0.2_singleR_pretty.pdf'))

实战部分(前面稍微看下就可以了)

下载数据

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE135927 这边在原始数据的custom自定义下载,也可以下载全部数据后,然后只筛选自己想要的部分,注意,下载完之后,可以选择手动解压或者代码解压,这边因为懒得查了,就直接手动解压,解压后把解压前的文件删除,每一个10x文件夹应该都对应着三个文件:一个是barcode文件(确定reads来源的细胞,有时可能为空,有时可能为多个细胞)、feature文件,matrix文件

将多个数据拆分并处理成seruat可以处理的文件格式

Seruat:read10x能读入的文件真的是只能读入tsv.gz,…等,且不能有任何的前缀名,所以这一步我们需要自己去处理,借用https://www.jianshu.com/p/284d61ba5d7c的代码

注意:这边的代码,我修改了部分,因为有时别人上传的文件命名规则不一样,基本是在folder的命名上,还有就是原博主不是用的最新的seruat

### 
# 列出当前目录下所有开头是GSM的文件
fs=list.files('./','^GSM')
# 然后获取四个样本信息
library(stringr)
samples=str_split(fs,'_',simplify = T)[,1]
# 设置一个循环,对每个样本信息做同样的事:
#(1)找到包含这个样本的文件(用grepl)
# (2)设置对应的目录名(str_split+paste)然后创建目录(用dir.create)
# (3)将文件放到对应目录(采用的是file.rename)并重命名文件

lapply(unique(samples),function(x){
  y=fs[grepl(x,fs)]
  folder=paste(str_split(y[1],'_',simplify = T)[,1],
               collapse = '')
  dir.create(folder,recursive = T)
  file.rename(y[1],file.path(folder,"barcodes.tsv.gz"))
  file.rename(y[2],file.path(folder,"features.tsv.gz"))
  file.rename(y[3],file.path(folder,"matrix.mtx.gz"))
})

批量读入写入10x数据并合并

注意:这边合并的规则,是按你需要合并的样本数指定的y,即根据额外的样本

rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
setwd("E:\\单细胞\\测试数据集\\")

# 分别读取每个10x样本的结果文件夹
if(F){ 
  samples=list.files('./')
  samples
  library(Seurat)
  sceList = lapply(samples,function(pro){ 
    folder=file.path('./',pro) 
    CreateSeuratObject(counts = Read10X(folder), 
                       project = pro )
  })
  sce.big <- merge(sceList[[1]], 
                   y = c(sceList[[2]]), 
                   add.cell.ids = samples, 
                   project = "ls_12")
  sce.big
  table(sce.big$orig.ident)
  save(sce.big,file = 'sce.big.merge.ls_12.Rdata')
  
}

可视化特征

注意:修改了核糖体的源代码,发现结果是一样的,所以源代码被我注释掉了,这步可视化的目的应该是跟后面的筛选特征相关,关于图怎么看,可以参考上面的链接

rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
## 合并成为一个R对象文件

load(file = 'sce.big.merge.ls_12.Rdata')
raw_sce=sce.big

raw_sce
rownames(raw_sce)[grepl('^mt-',rownames(raw_sce),ignore.case = T)]
rownames(raw_sce)[grepl('^Rp[sl]',rownames(raw_sce),ignore.case = T)]
# 计算线粒体,以及核糖体的比例,这边可能由于核糖体比较特殊,不能直接利用percentage函数?使用了一个迂回战术?
raw_sce[["percent.mt"]] <- PercentageFeatureSet(raw_sce, pattern = "^MT-")
fivenum(raw_sce[["percent.mt"]][,1])
raw_sce[["percent.ribo"]] <- PercentageFeatureSet(raw_sce, pattern = "^RP[SL]")
# rb.genes <- rownames(raw_sce)[grep("^RP[SL]",rownames(raw_sce))]
# C<-GetAssayData(object = raw_sce, slot = "counts")
# percent.ribo <- Matrix::colSums(C[rb.genes,])/Matrix::colSums(C)*100
# raw_sce <- AddMetaData(raw_sce, percent.ribo, col.name = "percent.ribo")

# 可视化特征,并可以在后面过滤条件做准备
plot1 <- FeatureScatter(raw_sce, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(raw_sce, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))

VlnPlot(raw_sce, features = c("percent.ribo", "percent.mt"), ncol = 2)
VlnPlot(raw_sce, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
VlnPlot(raw_sce, features = c("percent.ribo", "nCount_RNA"), ncol = 2)
raw_sce

筛选特征,过滤样本,并做标化

注意:过滤过程灵活变化,具体可能需要多看文献,NormalizeData是全局标化,假入有一个TPM的count矩阵,那么就可以不需要使用Seurat::NormalizeData()操作了,因为TPM已经根据测序深度进行了标准化,只需要进行log降一下维度即可。如果后续进行ScaleData操作,程序会检测是否使用了Seurat提供的标准化方法,如果我们使用自己的的标准化数据,那么就可能出现一个warning提醒,不过到时候不想被提醒,可以设置check.for.norm =F,findvariablefeatures 鉴定HVGs,使降维加快

pro='merge'
# 过滤标准
raw_sce1 <- subset(raw_sce, subset = nFeature_RNA > 200 & nCount_RNA > 1000 & percent.mt < 20)
raw_sce1

sce=raw_sce1
sce
sce <- NormalizeData(sce, normalization.method =  "LogNormalize", 
                     scale.factor = 10000)
GetAssay(sce,assay = "RNA")
sce <- FindVariableFeatures(sce, 
                            selection.method = "vst", nfeatures = 2000) #只是选取,并没有改变数据,或者在元数据的基础上多了一列
VariableFeaturePlot(sce)
# 步骤 ScaleData 的耗时取决于电脑系统配置(保守估计大于一分钟)
sce <- ScaleData(sce) 
# 标化后的材料放在sce@assays$RNA@scale.data,也就是说做PCA的时候用的是先找可变然后再做标化。那么问题来了先标化后可变或先可变再标化有无区别
# 关于这部分我大概试了一下先找可变后再标化和先标化再找可变然后观察标化后的材料发现,先标化再找可变此时标化的数据是原始数据行数,也就是说并没有挑选出可变后的材料。但是我想探讨的是先应该标化还是先应该找可变,由于这部分方法是vst即找方差,因为标化过后的话,方差都一样,所以其实确实应该先做找差异的部分

PCA与TSNE降维与聚类

降维中如何选择合适的主成分就是用elbowplot确定,DimHeatmap探索的是异质性来源,所谓的异质性来源我个人认为是关于主成分中哪些基因在不同的细胞间差异比较大。我认为这findneighbors与findcluster则是根据降维后的数据聚类了,dimplot这个用的是降维的输出,tsne降维的信息存放在seurat内,而pca则是RNA_snn_res,tSNE的参数设置:perplexity < (细胞数-1)/3,建议perplexity = 细胞数 / 50;

探索异质性来源—DimHeatmap

每个细胞和基因都根据PCA结果得分进行了排序,默认画前30个基因(nfeatures设置),1个主成分(dims设置),细胞数没有默认值(cells设置)

# PCA,PCA的对象是上面储存于rna中的
sce <- RunPCA(object = sce, pc.genes = VariableFeatures(sce)) 
DimHeatmap(sce, dims = 1:12, cells = 100, balanced = TRUE)
ElbowPlot(sce)
sce <- FindNeighbors(sce, dims = 1:15)
sce <- FindClusters(sce, resolution = 0.2)
table(sce@meta.data$RNA_snn_res.0.2)
# tsne
set.seed(123)
sce <- RunTSNE(object = sce, dims = 1:15, do.fast = TRUE)
DimPlot(sce,reduction = "tsne",label=T)
table(sce@meta.data$seurat_clusters,sce@meta.data$orig.ident)

phe=data.frame(cell=rownames(sce@meta.data),
               cluster =sce@meta.data$seurat_clusters)
head(phe)
table(phe$cluster)
tsne_pos=Embeddings(sce,'tsne')
DimPlot(sce,reduction = "tsne",label=T,split.by ='orig.ident')

head(phe)
table(phe$cluster)
head(tsne_pos) 
dat=cbind(tsne_pos,phe)
save(dat,file=paste0(pro,'_for_tSNE.pos.Rdata')) 

降完维且聚类后的图

load(file=paste0(pro,'_for_tSNE.pos.Rdata'))  
library(ggplot2)
p=ggplot(dat,aes(x=tSNE_1,y=tSNE_2,color=cluster))+geom_point(size=0.95)
p=p+stat_ellipse(data=dat,aes(x=tSNE_1,y=tSNE_2,fill=cluster,color=cluster),
                 geom = "polygon",alpha=0.2,level=0.9,type="t",linetype = 2,show.legend = F)+coord_fixed()
print(p) 
theme= theme(panel.grid =element_blank()) +   ## 删去网格
  theme(panel.border = element_blank(),panel.background = element_blank()) +   ## 删去外层边框
  theme(axis.line = element_line(size=1, colour = "black")) 
p=p+theme+guides(colour = guide_legend(override.aes = list(size=5)))
print(p)
ggplot2::ggsave(filename = paste0(pro,'_tsne_res0.2.pdf'))

plot1 <- FeatureScatter(sce, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(sce, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))
ggplot2::ggsave(filename = paste0(pro,'_CombinePlots.pdf'))

VlnPlot(sce, features = c("percent.ribo", "percent.mt"), ncol = 2)
ggplot2::ggsave(filename = paste0(pro,'_mt-and-ribo.pdf'))
VlnPlot(sce, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
ggplot2::ggsave(filename = paste0(pro,'_counts-and-feature.pdf'))
VlnPlot(sce, features = c("percent.ribo", "nCount_RNA"), ncol = 2)


table(sce@meta.data$seurat_clusters)

找细胞群的标志基因

for( i in unique(sce@meta.data$seurat_clusters) ){
  markers_df <- FindMarkers(object = sce, ident.1 = i, min.pct = 0.25)
  print(x = head(markers_df))
  markers_genes =  rownames(head(x = markers_df, n = 5))
  VlnPlot(object = sce, features =markers_genes,log =T )
  ggsave(filename=paste0(pro,'_VlnPlot_subcluster_',i,'_sce.markers_heatmap.pdf'))
  FeaturePlot(object = sce, features=markers_genes )
  ggsave(filename=paste0(pro,'_FeaturePlot_subcluster_',i,'_sce.markers_heatmap.pdf'))
}

sce.markers <- FindAllMarkers(object = sce, only.pos = TRUE, min.pct = 0.25, 
                              thresh.use = 0.25)

DT::datatable(sce.markers)
write.csv(sce.markers,file=paste0(pro,'_sce.markers.csv'))
library(dplyr) 
top10 <- sce.markers %>% group_by(cluster) %>% top_n(10, avg_log2FC)
DoHeatmap(sce,top10$gene,size=3)
ggsave(filename=paste0(pro,'_sce.markers_heatmap.pdf'))

save(sce,sce.markers,file = 'sce.output.merge.ls_12.Rdata')

细胞亚群自动化注释

注释时如果你希望去了解原理的话,可以翻一翻之前singR的原理。我个人的理解是,singR收录了部分不同类型的单细胞数据,并进行存储。它收录的单细胞数据是基于其他文献而进行的收集。
当我们的数据经tsne聚类后,会形成一个个簇。每个簇有很多细胞,对于细胞,会计算细胞与参考组的某个类型的细胞类型相关性,然后取分位数。作为与该簇的相关性,这部分用的是斯皮尔曼相关,也就是秩相关。那么对于秩相关而言不是取不取对数值都无所谓的么?怎么有文章说应该不取对数值。对于簇而言,有时对于细胞注释会有不同的单一结果,此时去掉占比少的。然后再进行进一步的细分计算剔除再计算,这么一想好像还是挺合理的
老鼠与人不一样

老鼠

老鼠的ps部分需要进行修改,其实有一点点问题,并不是所有类都刚好属于某一个类别,而且当主成分个数小于head的6个时,并不是很完美的分类方法

library(SingleR)
sce_for_SingleR <- GetAssayData(sce, slot="data")
mouseImmu <- ImmGenData()
pred.mouseImmu <- SingleR(test = sce_for_SingleR, ref = mouseImmu, labels = mouseImmu$label.main)

mouseRNA <- MouseRNAseqData()
pred.mouseRNA <- SingleR(test = sce_for_SingleR, ref = mouseRNA, labels = mouseRNA$label.fine )

cellType=data.frame(seurat=sce@meta.data$seurat_clusters,
                    mouseImmu=pred.mouseImmu$labels,
                    mouseRNA=pred.mouseRNA$labels)
sort(table(cellType[,2]))
sort(table(cellType[,3]))
table(cellType[,2:3])
save(sce,cellType, file=paste0(pro,'_output.Rdata') )
# 这边拿到了sce即原始表达数据,cellType(即根据TSNE、小鼠免疫,小鼠RNA的聚类标签)

rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
pro='S1'
library(SingleR)
load(file=paste0(pro,'_output.Rdata') )
# 把没有被分到免疫前6个的细胞设为other
cg=names(tail(sort(table(cellType[,2]))))
cellname=cellType[,2]
cellname[!cellname %in% cg ]='other'
table(sce@meta.data$seurat_clusters,cellname)
table(cellname)
# Epithelial cells, 0,2,5,7,8,12
# Endothelial cells, 10
# Fibroblasts , 1,3,6,9,14,15
# Macrophages , 4,13,16
# NKT , 11
# other ,4,11
# Stromal cells,
cl=sce@meta.data$seurat_clusters
ps=ifelse(cl %in% c(0,2,5,7,8,12),'epi',
          ifelse(cl %in% c(1,3,6,9,14,15),'fibro',
                 ifelse(cl %in% c(4,13,16),'macro',
                        ifelse(cl %in% c(10),'endo',
                               ifelse(cl %in% c(11),'NKT','other'
                               )))))
table(ps)
cellname=paste(cl,ps)
sce@meta.data$cellname=cellname
DimPlot(sce,reduction = "tsne",label=T,group.by = 'cellname')
ggplot2::ggsave(filename = paste0(pro,'_tsne_res0.2_singleR_raw.pdf'))
dat$cluster=cellname
library(ggplot2)
p=ggplot(dat,aes(x=tSNE_1,y=tSNE_2,color=cluster))+geom_point(size=0.95)
p=p+stat_ellipse(data=dat,aes(x=tSNE_1,y=tSNE_2,fill=cluster,color=cluster),
                 geom = "polygon",alpha=0.2,level=0.9,type="t",linetype = 2,show.legend = F)+coord_fixed()
print(p)
theme= theme(panel.grid =element_blank()) +   ## 删去网格
  theme(panel.border = element_blank(),panel.background = element_blank()) +   ## 删去外层边框
  theme(axis.line = element_line(size=1, colour = "black"))
p=p+theme+guides(colour = guide_legend(override.aes = list(size=5)))
print(p)
ggplot2::ggsave(filename = paste0(pro,'_tsne_res0.2_singleR_pretty.pdf'))
rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)

library(SingleR)
library("celldex")
hpca.se <- HumanPrimaryCellAtlasData()
hpca.se
library(SingleR)
sce_for_SingleR <- GetAssayData(sce, slot="data")
clusters=sce@meta.data$seurat_clusters

hpca.se <- HumanPrimaryCellAtlasData()
pred.hesc <- SingleR(test = sce_for_SingleR, ref = hpca.se, labels = hpca.se$label.fine,
                     method = "cluster", clusters = clusters, 
                     assay.type.test = "logcounts", assay.type.ref = "logcounts")

cellType=data.frame(ClusterID=levels(sce@meta.data$seurat_clusters),
                      mouseRNA=pred.hesc$labels )
head(cellType)
sce@meta.data$singleR=cellType[match(clusters,cellType$ClusterID),'mouseRNA']

pro='first_anno'
DimPlot(sce,reduction = "tsne",label=T, group.by = 'singleR')
ggplot2::ggsave(filename = paste0(pro,'_tSNE_singleR.pdf'))
DimPlot(sce,reduction = "tsne",label=T,split.by ='orig.ident',group.by = 'singleR')
ggplot2::ggsave(filename = paste0(pro,'_tSNE_singleR_by_orig.ident.pdf'))

save(sce,file = 'last_sce.Rdata') 
  • 3
    点赞
  • 23
    收藏
  • 打赏
    打赏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
©️2022 CSDN 皮肤主题:数字20 设计师:CSDN官方博客 返回首页
评论

打赏作者

愿航

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值