基礎(chǔ)流程（cellranger）(細(xì)胞管理員)

170336847_1_20190906045326454.jpg

cellranger 數(shù)據(jù)拆分

cellranger mkfastq可用于將單細(xì)胞測(cè)序獲得的 BCL 文件拆分為可以識(shí)別的 fastq 測(cè)序數(shù)據(jù)

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

-–run ：是下機(jī)數(shù)據(jù) BCL 所在的路徑沮尿；
-–samplesheet ：樣品信息列表--共三列（lane id ,sample name ,index name)
注意要控制好核心數(shù)和內(nèi)存數(shù)
運(yùn)行產(chǎn)出結(jié)果存在于 out 目錄中

cellranger 數(shù)據(jù)統(tǒng)計(jì)

cellranger count是 cellranger 最主要也是最重要的功能：完成細(xì)胞和基因的定量肢扯，也就是產(chǎn)生了我們用來(lái)做各種分析的基因表達(dá)矩陣聘裁。

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 ：產(chǎn)生的結(jié)果都在這個(gè)文件中，可以取幾號(hào)樣品（如 sample345）不翩；
fastqs ：由 cellranger mkfastq 產(chǎn)生的 fastqs 文件夾所在的路徑；fastqs ：由 cellranger mkfastq 產(chǎn)生的 fastqs 文件夾所在的路徑麻裳；
indices：sample index：SI-3A-A1慌盯；
transcriptome：參考轉(zhuǎn)錄組文件路徑；
cells：預(yù)期回復(fù)的細(xì)胞數(shù)掂器；

下游分析

cellranger count 計(jì)算的結(jié)果只能作為錯(cuò)略觀測(cè)的結(jié)果亚皂，如果需要進(jìn)一步分析聚類細(xì)胞，還需要進(jìn)行下游分析国瓮，這里使用官方推薦 R 包（Seurat 3.0）

軟件安裝

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

生成 Seruat 對(duì)象

library(dplyr)
library(Seurat)
#下載PBMC數(shù)據(jù)集
pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
#使用原始數(shù)據(jù)(非規(guī)范化數(shù)據(jù))初始化Seurat對(duì)象
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc

這里讀取的是單細(xì)胞 count 結(jié)果中的矩陣目錄灭必；
在對(duì)象生成的過(guò)程中，做了初步的過(guò)濾乃摹；
留下所有在>=3 個(gè)細(xì)胞中表達(dá)的基因 min.cells = 3禁漓；
留下所有檢測(cè)到>=200 個(gè)基因的細(xì)胞 min.genes = 200。
(為了除去一些質(zhì)量差的細(xì)胞)(為了除去一些質(zhì)量差的細(xì)胞)

標(biāo)準(zhǔn)預(yù)處理流程

使用原始數(shù)據(jù)(非規(guī)范化數(shù)據(jù))初始化Seurat對(duì)象孵睬。操作符可以向?qū)ο笤獢?shù)據(jù)添加列播歼。這是一個(gè)偉大的地方儲(chǔ)存QC統(tǒng)計(jì)。

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

這一步 mit-開(kāi)頭的為線粒體基因掰读，這里將其進(jìn)行標(biāo)記并統(tǒng)計(jì)其分布頻率

# 將QC指標(biāo)可視化為小提琴圖
VlnPlot(object = pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

對(duì) pbmc 對(duì)象做小提琴圖秘狞，分別為基因數(shù)，細(xì)胞數(shù)和線粒體占比.

170336847_2_20190906045326641.jpg

接下來(lái)蹈集，根據(jù)圖片中基因數(shù)和線粒體數(shù)烁试，分別設(shè)置過(guò)濾參數(shù)，這里基因數(shù) 200-2500拢肆，線粒體百分比為小于 5%

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

數(shù)據(jù)標(biāo)準(zhǔn)化

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

鑒定高度變化基因

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

繪制 variable features 的帶有和不帶有標(biāo)簽的圖形

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

170336847_3_20190906045327501.jpg

數(shù)據(jù)歸一化

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

線性降維

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

這里有多種方法展示 pca 結(jié)果减响，本文采用最簡(jiǎn)單的方法

DimPlot(object = pbmc, reduction = "pca")

untitled.png

鑒定數(shù)據(jù)集的可用維度

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

untitled.png

虛線以上的為可用維度靖诗，你也可以調(diào)整 dims 參數(shù)，畫出所有 pca 查看

細(xì)胞聚類

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

這里的 dims 為上一步計(jì)算所用的維度數(shù)支示，而 resolution 參數(shù)控制聚類的數(shù)目刊橘，這里用默認(rèn)就好

執(zhí)行非線性降維

這里注意，這一步聚類有兩種聚類方法(umap/tSNE)颂鸿，兩種方法都可以使用伤为，但不要混用，這樣据途，后面的結(jié)算結(jié)果會(huì)將先前的聚類覆蓋掉绞愚，只能保留一個(gè).
本文采用基于圖論的聚類方法

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

untitled.png

完成聚類后，一定要記住保存數(shù)據(jù)颖医，不然重新計(jì)算可要頭疼了

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

如下為TSNE聚類

TSNE聚類分析

pcSelect=20
pbmc <- FindNeighbors(object = pbmc, dims = 1:pcSelect)                #計(jì)算鄰接距離
pbmc <- FindClusters(object = pbmc, resolution = 0.5)                  #對(duì)細(xì)胞分組,優(yōu)化標(biāo)準(zhǔn)模塊化
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)

尋找每個(gè)聚類中顯著表達(dá)的基因

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

這樣是尋找單個(gè)聚類中的顯著基因

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

這樣尋找所有聚類中顯著基因位衩，計(jì)算速度很慢，需要等待.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)

有多種方法統(tǒng)計(jì)基因的顯著性

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

針對(duì)前十個(gè)基因做熱圖熔萧。

尋找基因 marker 并對(duì)細(xì)胞類型進(jìn)行注釋

全自動(dòng)細(xì)胞類型注釋

SingleR:一個(gè)全自動(dòng)細(xì)胞注釋的 R 包糖驴，用法很簡(jiǎn)單
軟件安裝

devtools::install_github('dviraran/SingleR')

這可能需要很長(zhǎng)時(shí)間，盡管主要是因?yàn)镾eurat的安裝佛致。

創(chuàng)建 SingleR 對(duì)象

官方有多種方法創(chuàng)建該對(duì)象贮缕，參考SingleR - create object
我們這里由于已經(jīng)具有了 Seurat 對(duì)象，所以可以采用直接轉(zhuǎn)化的方法

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)

對(duì)于S4對(duì)象俺榆，需要手動(dòng)尋找數(shù)據(jù)

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

fine.tune 如果設(shè)置為 T感昼，會(huì)消耗大量時(shí)間，這一步是對(duì)數(shù)據(jù)小差異的進(jìn)一步細(xì)化罐脊，可以不計(jì)算
do.signatures 這個(gè)也會(huì)消耗大量時(shí)間定嗓，做單細(xì)胞基因集豐度分析，可以先設(shè)置為 F

對(duì)象載入完成就可以保存好去官方網(wǎng)站進(jìn)行可視化分析了

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

untitled.png

偽時(shí)間分析

偽時(shí)間分析建議采用 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的安裝過(guò)程中已經(jīng)安裝過(guò)的就不必安裝了

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")

偽時(shí)間分析

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

如果使用的是seurat2.4版本萍桌，可以使用monocle的importCDS命令直接導(dǎo)入宵溅，如果是3.0版本，需要進(jìn)行如下手動(dòng)導(dǎo)入數(shù)據(jù)

這里采用的是官方教程中所需要的三個(gè)文件上炎，細(xì)胞矩陣恃逻，細(xì)胞注釋表和基因注釋表表

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)

給其中一列數(shù)據(jù)重命名

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

添加細(xì)胞聚類數(shù)據(jù)

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")

偽時(shí)間分析流程

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)

選擇特定細(xì)胞進(jìn)行偽時(shí)間分析

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

偽時(shí)間分析

untitled.png

特定細(xì)胞分析