文章目录

  • 写在前面
    • 写作来源
    • 开始
      • 为什么要用单细胞测序
      • 原始数据产出形式
      • 单细胞基本流程
      • 一些背景知识点
        • 为什么不使用RPKM等类似的统计量
        • CV（标准差/均值）
        • 单细胞的样本配置
        • 质控部分
        • 标化问题
        • tsne与pca
        • 为什么要去除核糖体和线粒体
        • UMI，indrops`, `SEQC`, `zUMIs
        • 原理与流程
        • 应用方向
      • R包需求等
        • 环境的设置
        • R包的安装
      • 代码实操
        • 关于单个10x的祖传代码
      • 实战部分（前面稍微看下就可以了）
        • 下载数据
        • 将多个数据拆分并处理成seruat可以处理的文件格式
        • 批量读入写入10x数据并合并
        • 可视化特征
        • 筛选特征,过滤样本，并做标化
        • PCA与TSNE降维与聚类
        • 降完维且聚类后的图
        • 找细胞群的标志基因
        • 细胞亚群自动化注释
          • 老鼠
          • 人

写在前面

第一次尝试用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_scehp_sce1 <- subset(hp_sce, subset = nFeature_RNA > 200 & nCount_RNA > 1000 & percent.mt < 20)
hp_sce1sce=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('./')sampleslibrary(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.bigtable(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.bigraw_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_sce1sce=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_clustershpca.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')

生信-单细胞数据处理相关推荐

  1. 肿瘤/非肿瘤/单基因/单细胞/非编码：史上最全生信分析攻略！！！

    解读生信之美,探讨每篇文献背后的逻辑 非肿瘤专栏:条条大路通罗马 1.4+非肿瘤生信分析+铁死亡/焦亡/自噬/代谢/免疫的万能钥匙 短评:适合一些热门机制如铁死亡/焦亡/自噬等在非肿瘤疾病中的研究 2 ...

  2. 肿瘤NGS生信知识来源-博客、公众号、网站

    记录一些人,防止失联 1.肿瘤方向生信工程师博客 2.肿瘤NGS相关公众号 3.肿瘤NGS相关咨询网站 4.编程基础信息来源 5.其他值得推荐的 6.主流肿瘤NGS基因测序公司 7.其他信息来源    ...

  3. 有了易生信，导师再也不用担心我的单细胞转录组整合分析啦

    2019年10月9日,单细胞转录组再等Nature.题为Decoding human fetal liver haematopoiesis的研究,对受孕后4周至17周的人胚胎肝脏.卵黄囊.肾脏和皮肤组 ...

  4. 易生信群体和单细胞转录组专题第6期于5月10日在北京开课了

    群体转录组是我们最常接触到的一种高通量测序数据类型,其实验方法成熟,花费较低,分析思路简洁清晰,是入门生信,解决最常见问题的首选. 单细胞分析是近几年的明星技术,多次被Nature.Science评为 ...

  5. 生信入门（六）——单细胞分析（Seurat）

    生信入门(六)--单细胞分析(Seurat) 文章目录 生信入门(六)--单细胞分析(Seurat) 一.数据导入 1.数据来源 2.数据导入 二.标准预处理 1.QC和选择细胞进行进一步分析 2.规 ...

  6. 有没有人带？这些都是学习生信的一大助力！

    经常能看到某人5篇SCI! 某实验室10篇SCI! 科研学习过程中,一年多篇文章的人是如何保持科研创新能力的? 及时关注学界动态,高效获取优质资源非常重要. 但是大量的科研资讯,热点也层出不穷,那么作 ...

  7. perl语言入门第七版 电子版_百迈客带您走近生信分析【入门篇】

    年末促销倒计时:59天 百迈客推出年末活动促销,发文有礼.推广有礼以及多种产品钜惠来袭,百迈客为您倾情打造科研福利,您还在等什么?快快行动起来领取您的超级奖励吧!(详情请见"决战2020!品 ...

  8. ​北京大学吴华君组诚聘医学/生信助理研究员和博士后

    一.研究方向及合作导师 吴华君,北京大学医学部精准医疗多组学研究中心PI.博士生导师,北京大学肿瘤医院双聘研究员.2012年在中科院遗传与发育所获得博士学位,随后在美国哈佛大学从事博士后研究,在Can ...

  9. 微生态、生信和植物领域最新资讯合集，不看你就亏大啦！！！

    宏基因组/微生物组是当今世界科研最热门的研究领域之一,中科院科研人员创立"宏基因组"公众号,入选科研圈评选"2019年度学术媒体优质公众号",联合海内外同行共同 ...

最新文章

  1. php使用redis命令,PHP 使用 Redis
  2. Codeforces Round #703 (Div. 2) E. Paired Payment 最短路 + 思维
  3. 单例设计模式–内省和最佳实践
  4. 正方形矩阵求对角线之和
  5. ip地址转换pta题目_PTA「实验2-3-5 输出华氏-摄氏温度转换表」
  6. codeforces 758 A
  7. python中的缩进是长度统一吗_python缩进长度是否统一
  8. JZOJ5787轨道（容斥+DP）
  9. MagicAjax的用法, 每10秒刷新, 更改等待loading效果
  10. C++与Python混合编程
  11. robot framework 添加selenium2library显示红色
  12. dvr行业的linux
  13. 如何使用浏览器的网页全文翻译工具
  14. Windows 最全CMD命令，带死机修复系统命令
  15. Word粘贴时出现“文件未找到：MathPage.WLL”的解决方案
  16. 最大规模开源说话人识别语料集——VoxCeleb
  17. php 读取纯真书库,PHP读取纯真IP数据库的函数
  18. C++初学必练基础题【第二期】
  19. 浙江理工大学英语平台Unipus自动答题
  20. 解决‘from sklearn.inspection import DecisionBoundaryDisplay‘行代码爆红问题

热门文章

  1. 背景与字体的搭配经验
  2. 字符串库函数（1）Strlen，strcpy，strcat，strcmp
  3. 美图秀秀开发插件生成的图片都有哪些格式？
  4. Attention 二 创新篇
  5. 用canvas画一个炫酷的粒子动画倒计时
  6. 前端工程师ps学习笔记
  7. 【安全牛学习笔记】TearDrop
  8. 两全险的主险是什么意思？
  9. 程序员通过google 赶快来赚美金~ hot
  10. 苹果电脑python编程里面怎么切到中文_mac下的 idle为何不能输入中文?该如何解决?...