单细胞 R语言单细胞数据处理_mob6454cc78d412的技术博客_

重磅综述：三万字长文读懂单细胞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 单细胞里的对象说明：https://www.jianshu.com/p/1db7c28249d4

开始

为什么要用单细胞测序

细胞之间存在差异，比如肿瘤的中心，周围，淋巴转移灶有不同的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降解
有文章建议线粒体取15%以下较为合理，但有文章取到了30%（Single-cell RNA sequencing reveals distinct tumor microenvironmental patterns in lung adenocarcinoma）
核糖体我看都比较低，一般取0-5

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
### 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/'),
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 <- 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,使降维加快
使用sctransform标化时，最好先对单个样本进行标化，而后再合并，参考链接https://github.com/satijalab/seurat/issues/8331，https://satijalab.org/seurat/archive/v4.3/sctransform_v2_vignette

pro='merge'
# 过滤标准
raw_sce1 <- subset(raw_sce, subset = nFeature_RNA > 200 & nCount_RNA > 1000 & percent.mt < 20)
raw_sce1
sce=raw_sce1
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聚类后，会形成一个个簇。每个簇有很多细胞，对于细胞，会计算细胞与参考组的某个类型的细胞类型相关性，然后取分位数。作为与该簇的相关性，这部分用的是斯皮尔曼相关，也就是秩相关。那么对于秩相关而言不是取不取对数值都无所谓的么？怎么有文章说应该不取对数值。对于簇而言，有时对于细胞注释会有不同的单一结果，此时去掉占比少的。然后再进行进一步的细分计算剔除再计算，这么一想好像还是挺合理的
老鼠与人不一样；不建议进行自动化注释

细胞手动注释

先将细胞注释成几个大类

其他文献常见单细胞聚类初图

单细胞 R语言单细胞数据处理_数据

单细胞 R语言单细胞数据处理_单细胞 R语言_02

单细胞组织样本与血液样本的区别

以下检索答案来自于chatGPT:
单细胞血液数据和组织数据在样本来源和研究对象等方面有一些区别：

样本来源：单细胞血液数据来自于外周血样本，而单细胞组织数据来自于组织样本，例如肝脏、肾脏、肺部等。

细胞类型：单细胞血液数据主要包含来自血液中的各种免疫细胞（如T细胞、B细胞、单核细胞等），以及一些其他细胞类型（如红细胞、血小板等）。而单细胞组织数据则包含特定组织中的各种细胞类型，如肝脏中的肝细胞、肾脏中的肾小管细胞等。

细胞状态：由于环境和功能的不同，单细胞血液数据中的细胞可能处于不同的状态，如激活状态、静止状态等。而单细胞组织数据中的细胞可能受到组织微环境的影响，处于特定的分化或功能状态。

研究目的：单细胞血液数据通常用于研究免疫系统的功能、免疫细胞的分化和活化过程等。而单细胞组织数据则可以用于研究特定组织中不同细胞类型的分布、相互作用以及在疾病状态下的变化等。

综上所述，单细胞血液数据和组织数据在来源、细胞类型、状态和研究目的等方面有所不同，研究人员在选择数据来源和分析方法时需要考虑这些差异。

为什么NK细胞与T细胞难以区分

在生物信息学中，单细胞T细胞和NK细胞之间难以区分的主要原因是它们在基因组水平上的相似性较高。T细胞和NK细胞都是免疫系统中的重要成分，具有类似的功能和表型特征。在单细胞转录组数据中，由于这两种细胞类型之间的转录组差异不是很大，很难通过基因表达模式来清晰地区分它们。

为了更准确地区分单细胞T细胞和NK细胞，研究人员通常会结合细胞表面标记物（如CD3、CD56等）和功能基因（如IFNG、GZMB等）的表达模式来进行区分。此外，使用细胞亚群分析的方法，如聚类分析和降维技术（如t-SNE、PCA等），也可以帮助识别和区分这两种细胞类型。

JVM内存使用率和容器内存使用率

java 语言, 开发者不能直接控制程序运行内存, 对象的创建都是由类加载器一步步解析, 执行与生成与内存区域中的; 并且jvm有自己的垃圾回收器对内存区域管理, 回收; 但是我们已经可以通过一些工具来在程序运行时查看对应的jvm内存使用情况, 帮助更好的分析与优化我们的代码;注: 查看系统里java进程信息// 查看当前机器上所有运行的java进程名称与pid(进程编号) jps -l

文章目录

写在前面

写作来源

开始