单细胞分析1（一文读懂单细胞测序分析流程）

基础流程（cellranger）(细胞管理员)

170336847_1_20190906045326454.jpg

cellranger 数据拆分

cellranger mkfastq可用于将单细胞测序获得的 BCL 文件拆分为可以识别的 fastq 测序数据

cellranger makefastq   --run=[ ]   --samplesheet=[sample.csv] --jobmode=local --localcores=20 --localmem=80

-–run ：是下机数据 BCL 所在的路径；
-–samplesheet ：样品信息列表--共三列（lane id ,sample name ,index name)
注意要控制好核心数和内存数
运行产出结果存在于 out 目录中

cellranger 数据统计

cellranger count是 cellranger 最主要也是最重要的功能：完成细胞和基因的定量，也就是产生了我们用来做各种分析的基因表达矩阵。

cellranger count -–id=sample345 -–transcriptome=/opt/refdata-cellranger-GRCh38-1.2.0/GRCh38 -–fastqs=/home/jdoe/runs/HAWT7ADXX/outs/fastq_path -–indices=SI-3A-A1 –-cells=1000

id ：产生的结果都在这个文件中，可以取几号样品（如 sample345）；
fastqs ：由 cellranger mkfastq 产生的 fastqs 文件夹所在的路径；fastqs ：由 cellranger mkfastq 产生的 fastqs 文件夹所在的路径；
indices：sample index：SI-3A-A1；
transcriptome：参考转录组文件路径；
cells：预期回复的细胞数；

下游分析

cellranger count 计算的结果只能作为错略观测的结果，如果需要进一步分析聚类细胞，还需要进行下游分析，这里使用官方推荐 R 包（Seurat 3.0）

软件安装

install.packages('devtools')
devtools::install_github(repo = 'satijalab/seurat', ref = 'release/3.0')
library(Seurat)

生成 Seruat 对象

library(dplyr)
library(Seurat)
#下载PBMC数据集
pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
#使用原始数据(非规范化数据)初始化Seurat对象
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc

这里读取的是单细胞 count 结果中的矩阵目录；
在对象生成的过程中，做了初步的过滤；
留下所有在>=3 个细胞中表达的基因 min.cells = 3；
留下所有检测到>=200 个基因的细胞 min.genes = 200。
(为了除去一些质量差的细胞)(为了除去一些质量差的细胞)

标准预处理流程

使用原始数据(非规范化数据)初始化Seurat对象。操作符可以向对象元数据添加列。这是一个伟大的地方储存QC统计。

pbmc[["percent.mt"]] <- PercentageFeatureSet(object = pbmc, pattern = "^MT-")

这一步 mit-开头的为线粒体基因，这里将其进行标记并统计其分布频率

# 将QC指标可视化为小提琴图
VlnPlot(object = pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

对 pbmc 对象做小提琴图，分别为基因数，细胞数和线粒体占比.


170336847_2_20190906045326641.jpg

接下来，根据图片中基因数和线粒体数，分别设置过滤参数，这里基因数 200-2500，线粒体百分比为小于 5%

pbmc <- subset(x = pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

数据标准化

pbmc <- NormalizeData(object = pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- NormalizeData(object = pbmc)

鉴定高度变化基因

pbmc <- FindVariableFeatures(object = pbmc, selection.method = "vst", nfeatures = 2000)
#找出10个变化最大的基因
top10 <- head(x = VariableFeatures(object = pbmc), 10)

绘制 variable features 的带有和不带有标签的图形

plot1 <- VariableFeaturePlot(object = pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
CombinePlots(plots = list(plot1, plot2))
170336847_3_20190906045327501.jpg

数据归一化

all.genes <- rownames(x = pbmc)
pbmc <- ScaleData(object = pbmc, features = all.genes)

线性降维

pbmc <- RunPCA(object = pbmc, features = VariableFeatures(object = pbmc))

这里有多种方法展示 pca 结果，本文采用最简单的方法

DimPlot(object = pbmc, reduction = "pca")
untitled.png

鉴定数据集的可用维度

pbmc <- JackStraw(object = pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(object = pbmc, dims = 1:20)
JackStrawPlot(object = pbmc, dims = 1:15)
untitled.png

虚线以上的为可用维度，你也可以调整 dims 参数，画出所有 pca 查看

细胞聚类

pbmc <- FindNeighbors(object = pbmc, dims = 1:10)
pbmc <- FindClusters(object = pbmc, resolution = 0.5)

这里的 dims 为上一步计算所用的维度数，而 resolution 参数控制聚类的数目，这里用默认就好

执行非线性降维

这里注意，这一步聚类有两种聚类方法(umap/tSNE)，两种方法都可以使用，但不要混用，这样，后面的结算结果会将先前的聚类覆盖掉，只能保留一个.
本文采用基于图论的聚类方法

pbmc <- RunUMAP(object = pbmc, dims = 1:10)
DimPlot(object = pbmc, reduction = "umap")
untitled.png

完成聚类后，一定要记住保存数据，不然重新计算可要头疼了

saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")

如下为TSNE聚类

TSNE聚类分析

pcSelect=20
pbmc <- FindNeighbors(object = pbmc, dims = 1:pcSelect)                #计算邻接距离
pbmc <- FindClusters(object = pbmc, resolution = 0.5)                  #对细胞分组,优化标准模块化
pbmc <- RunTSNE(object = pbmc, dims = 1:pcSelect)                      #TSNE聚类
pdf(file="06.TSNE.pdf",width=6.5,height=6)
TSNEPlot(object = pbmc, do.label = TRUE, pt.size = 2, label = TRUE)    #TSNE可视化
dev.off()
write.table(pbmc$seurat_clusters,file="06.tsneCluster.txt",quote=F,sep="\t",col.names=F)

寻找每个聚类中显著表达的基因

cluster1.markers <- FindMarkers(object = pbmc, ident.1 = 1, min.pct = 0.25)
head(x = cluster1.markers, n = 5)

这样是寻找单个聚类中的显著基因

cluster5.markers <- FindMarkers(object = pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(x = cluster5.markers, n = 5)

这样寻找所有聚类中显著基因，计算速度很慢，需要等待.find markers for every cluster compared to all remaining cells, report only the positive ones

pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)

有多种方法统计基因的显著性

FeaturePlot(object = pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ",
    "PPBP", "CD8A"))
untitled.png
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(object = pbmc, features = top10$gene) + NoLegend()
untitled.png

针对前十个基因做热图。

寻找基因 marker 并对细胞类型进行注释

全自动细胞类型注释

SingleR:一个全自动细胞注释的 R 包，用法很简单
软件安装

devtools::install_github('dviraran/SingleR')

这可能需要很长时间，尽管主要是因为Seurat的安装。

创建 SingleR 对象

官方有多种方法创建该对象，参考SingleR - create object
我们这里由于已经具有了 Seurat 对象，所以可以采用直接转化的方法

library(SingleR)

singler = CreateSinglerObject(counts, annot = NULL, project.name, min.genes = 0,
  technology = "10X", species = "Human", citation = "",
  ref.list = list(), normalize.gene.length = F, variable.genes = "de",
  fine.tune = T, do.signatures = T, clusters = NULL, do.main.types = T,
  reduce.file.size = T, numCores = SingleR.numCores)

singler$seurat = seurat.object # (optional)
singler$meta.data$orig.ident = seurat.object@meta.data$orig.ident # the original identities, if not supplied in 'annot'

if using Seurat v3.0 and over use:

singler$meta.data$xy = seurat.object@reductions$tsne@cell.embeddings # the tSNE coordinates
singler$meta.data$clusters = seurat.object@active.ident # the Seurat clusters (if 'clusters' not provided)

对于S4对象，需要手动寻找数据

counts<-seurat.object@assays$RNA@counts
clusters<-seurat.object@meta.data$seurat_clusters

fine.tune 如果设置为 T，会消耗大量时间，这一步是对数据小差异的进一步细化，可以不计算
do.signatures 这个也会消耗大量时间，做单细胞基因集丰度分析，可以先设置为 F

对象载入完成就可以保存好去官方网站进行可视化分析了

singler.new = convertSingleR2Browser(singler)
saveRDS(singler.new,file=paste0(singler.new@project.name,'.rds')
untitled.png

伪时间分析

伪时间分析建议采用 monocle3.0 软件

软件安装

source("http://bioconductor.org/biocLite.R")
biocLite("monocle")
devtools::install_github("cole-trapnell-lab/DDRTree", ref="simple-ppt-like")
devtools::install_github("cole-trapnell-lab/L1-graph")

这一步在Seurat3.0的安装过程中已经安装过的就不必安装了

install.packages("reticulate")
library(reticulate)
py_install('umap-learn', pip = T, pip_ignore_installed = T) # Ensure the latest version of UMAP is installed
py_install("louvain")
devtools::install_github("cole-trapnell-lab/monocle-release", ref="monocle3_alpha")

伪时间分析

library(Seurat)
library(monocle)
Seurat.obj<-readRDS("**.rds")

如果使用的是seurat2.4版本，可以使用monocle的importCDS命令直接导入，如果是3.0版本，需要进行如下手动导入数据

这里采用的是官方教程中所需要的三个文件，细胞矩阵，细胞注释表和基因注释表表

data <- as(as.matrix(Seurat.obj@assays$RNA@data), 'sparseMatrix')
pd<-new("AnnotatedDataFrame", data = Seurat.obj@meta.data)
fd<-new("AnnotatedDataFrame", data = data.frame(gene_short_name = row.names(data), row.names = row.names(data)))
cds <- newCellDataSet(data, phenoData = pd, featureData = fd)

给其中一列数据重命名

names(pData(cds))[names(pData(cds))=="RNA_snn_res.0.5"]="Cluster"

添加细胞聚类数据

pData(cds)$cell_type2 <- plyr::revalue(as.character(pData(cds)$Cluster),c("0" = 'Fibroblasts',"1" = 'Fibroblasts',"2" = 'Fibroblasts',"3" = 'Fibroblasts',"4" = 'Fibroblasts',"5" = 'NK',"6" = 'Fibroblasts',"7" = 'Macrophage',"8" = 'NK',"9" = 'Macrophage',"10" = 'EC',"11" = 'Fibroblasts',"12" = 'EC'))
cell_type_color <- c("Fibroblasts" = "#E088B8","NK" = "#46C7EF","Macrophage" = "#EFAD1E","EC" = "#8CB3DF")

伪时间分析流程

cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
cds <- preprocessCDS(cds, num_dim = 20)
cds <- reduceDimension(cds, reduction_method = 'UMAP')
cds <- partitionCells(cds)
cds <- learnGraph(cds,  RGE_method = 'SimplePPT')
plot_cell_trajectory(cds,color_by = "cell_type2") + scale_color_manual(values = cell_type_color)

选择特定细胞进行伪时间分析

get_correct_root_state <- function(cds, cell_phenotype, root_type){
  cell_ids <- which(pData(cds)[, cell_phenotype] == root_type)
  closest_vertex <-cds@auxOrderingData[[cds@rge_method]]$pr_graph_cell_proj_closest_vertex
  closest_vertex <- as.matrix(closest_vertex[colnames(cds), ])
  root_pr_nodes <-V(cds@minSpanningTree)$name[as.numeric(names(which.max(table(closest_vertex[cell_ids,]))))]
}
MPP_node_ids = get_correct_root_state(cds,cell_phenotype ='cell_type2', 'Fibroblasts')
cds <- orderCells(cds, root_pr_nodes = MPP_node_ids)
plot_cell_trajectory(cds)
untitled.png

伪时间分析


untitled.png

特定细胞分析

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 213,254评论 6赞 492
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 90,875评论 3赞 387
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 158,682评论 0赞 348
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 56,896评论 1赞 285
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 66,015评论 6赞 385
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 50,152评论 1赞 291
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 39,208评论 3赞 412
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 37,962评论 0赞 268
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 44,388评论 1赞 304
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,700评论 2赞 327
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,867评论 1赞 341
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 34,551评论 4赞 335
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 40,186评论 3赞 317
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,901评论 0赞 21
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,142评论 1赞 267
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 46,689评论 2赞 362
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 43,757评论 2赞 351

推荐阅读更多精彩内容

  • 10X单细胞测序分析软件:Cell ranger
    单细胞测序有着漫长的过去，却只有短暂的历史---谁说的！ 说她漫长是因为到如今也有十几年的历史了，说她段短暂是因为...
    周运来就是我阅读 56,518评论 45赞 123
  • 10×单细胞测序分析练习（一）
    目前单细胞测序获取细胞的方法主要有两种：（1）Smart-seq：测序细胞量少(因为一次通过一个细胞)，但是测序的...
    生信start_site阅读 22,243评论 6赞 55
  • 正念2
    尼采所说：婚姻生活犹如长期的对话。当你要迈进婚姻生活时，一定要先这样反问自己：你是否能和这位女子在白头偕老时，仍能...
    蜻蜓严阅读 182评论 0赞 1
  • 文本分类-分类器构建与测试
    二、构建分类器## 2.1用logistic回归训练出分类器利用文本预处理得到的训练集的特征向量集、标签集trai...
    Shira0905阅读 878评论 0赞 0
  • 一年
    这一年 那一年 一年又一年 细数今年，回忆那年 展望明年，畅想后年 祈祷祝福每一年 一年 装着365天 走过春夏秋...
    一棹碧涛阅读 633评论 7赞 17