|
學(xué)習(xí)單細(xì)胞轉(zhuǎn)錄組肯定先來一遍Seurat官網(wǎng)的標(biāo)準(zhǔn)流程。 數(shù)據(jù)來源于Peripheral Blood Mononuclear Cells (PBMC)，共2700個(gè)單細(xì)胞， Illumina NextSeq 500平臺(tái)。下載鏈接在這：https://cf./samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz 解壓到目錄下，有以下三個(gè)文件，也是10X的標(biāo)準(zhǔn)文件 barcodes.tsv ,genes.tsv, matrix.mtx 一 加載R包 數(shù)據(jù)集 Seurat可以直接用Read10X函數(shù)讀取cellranger的結(jié)果數(shù)據(jù)，使用pbmc數(shù)據(jù)初始化Seurat對(duì)象 library(dplyr)library(Seurat) library(patchwork) # Load the PBMC dataset pbmc.data <- Read10X(data.dir = "./data/hg19/") #解壓縮后的路徑 # Initialize the Seurat object with the raw (non-normalized data). pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200) pbmc 查看count matrix矩陣 #查看數(shù)據(jù)dim(pbmc.data) # 32738 2700 # Lets examine a few genes in the first thirty cells pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30] 3 x 30 sparse Matrix of class "dgCMatrix" [[ suppressing 30 column names 'AAACATACAACCAC-1’, 'AAACATTGAGCTAC-1’, 'AAACATTGATCAGC-1’ ... ]] CD3D 4 . 10 . . 1 2 3 1 . . 2 7 1 . . 1 3 . 2 3 . . . . . 3 4 1 5 TCL1A . . . . . . . . 1 . . . . . . . . . . . . 1 . . . . . . . . MS4A1 . 6 . . . . . . 1 1 1 . . . . . . . . . 36 1 2 . . 2 . . . . 其中的.表示0，即no molecules detected。seurat使用一個(gè)稀疏矩陣來保存count matrix，節(jié)省存儲(chǔ)空間。 二 數(shù)據(jù)預(yù)處理 2.1 QC 一般使用以下三個(gè)標(biāo)準(zhǔn)，也可以參考commonly used
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-") # Show QC metrics for the first 5 cells head(pbmc@meta.data, 5)## orig.ident nCount_RNA nFeature_RNA percent.mt ## AAACATACAACCAC-1 pbmc3k 2419 779 3.0177759 ## AAACATTGAGCTAC-1 pbmc3k 4903 1352 3.7935958 ## AAACATTGATCAGC-1 pbmc3k 3147 1129 0.8897363 ## AAACCGTGCTTCCG-1 pbmc3k 2639 960 1.7430845 ## AAACCGTGTATGCG-1 pbmc3k 980 521 1.2244898 注意PercentageFeatureSet函數(shù)可以計(jì)算每個(gè)CELL中的指定基因子集的計(jì)數(shù)百分比。 小提琴圖：查看基因數(shù)目, UMI數(shù)目, 線粒體基因占比 # Visualize QC metrics as a violin plotVlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3) 基因數(shù)目, 線粒體基因占比與UMI數(shù)目相關(guān)性 # FeatureScatter is typically used to visualize feature-feature relationships, but can be used# for anything calculated by the object, i.e. columns in object metadata, PC scores etc. plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt") plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") plot1 + plot2 過濾
2.2 數(shù)據(jù)標(biāo)準(zhǔn)化默認(rèn)標(biāo)準(zhǔn)化方法為L(zhǎng)ogNormalize，標(biāo)化后的數(shù)據(jù)存在pbmc[["RNA"]]@data中。 pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)#pbmc <- NormalizeData(pbmc) 如果我記不住這些[[ 和 @ 怎么辦？ 使用最“簡(jiǎn)單，通用”的方式，str(pbmc) 查看一下數(shù)據(jù)結(jié)構(gòu)，然后根據(jù)結(jié)構(gòu)對(duì)應(yīng)使用@ 或者 $ 。 str(pbmc)#截取部分展示用 Formal class 'Seurat' [package "SeuratObject"] with 13 slots ..@ assays :List of 1 .. ..$ RNA:Formal class 'Assay' [package "SeuratObject"] with 8 slots .. .. .. ..@ counts :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots .. .. .. .. .. ..@ i : int [1:2238732] 29 73 80 148 163 184 186 227 229 230 ... .. .. .. .. .. ..@ p : int [1:2639] 0 779 2131 3260 4220 4741 5522 6304 7094 7626 ... .. .. .. .. .. ..@ Dim : int [1:2] 13714 2638 .. .. .. .. .. ..@ Dimnames:List of 2 .. .. .. .. .. .. ..$ : chr [1:13714] "AL627309.1" "AP006222.2" "RP11-206L10.2" "RP11-206L10.9" ... .. .. .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ... .. .. .. .. .. ..@ x : num [1:2238732] 1 1 2 1 1 1 1 41 1 1 ... .. .. .. .. .. ..@ factors : list() .. .. .. ..@ data :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots .. .. .. .. .. ..@ i : int [1:2238732] 29 73 80 148 163 184 186 227 229 230 ... .. .. .. .. .. ..@ p : int [1:2639] 0 779 2131 3260 4220 4741 5522 6304 7094 7626 ... .. .. .. .. .. ..@ Dim : int [1:2] 13714 2638 .. .. .. .. .. ..@ Dimnames:List of 2 .. .. .. .. .. .. ..$ : chr [1:13714] "AL627309.1" "AP006222.2" "RP11-206L10.2" "RP11-206L10.9" ... .. .. .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ... .. .. .. .. .. ..@ x : num [1:2238732] 1.64 1.64 2.23 1.64 1.64 ... .. .. .. .. .. ..@ factors : list() .. .. .. ..@ scale.data : num[0 , 0 ] .. .. .. ..@ key : chr "rna_" .. .. .. ..@ assay.orig : NULL .. .. .. ..@ var.features : chr(0) .. .. .. ..@ meta.features:'data.frame':13714 obs. of 0 variables .. .. .. ..@ misc : list() ..@ meta.data :'data.frame':2638 obs. of 5 variables: .. ..$ orig.ident : Factor w/ 1 level "pbmc3k": 1 1 1 1 1 1 1 1 1 1 ... .. ..$ nCount_RNA : num [1:2638] 2419 4903 3147 2639 980 ... .. ..$ nFeature_RNA: int [1:2638] 779 1352 1129 960 521 781 782 790 532 550 ... .. ..$ percent.mt : num [1:2638] 3.02 3.79 0.89 1.74 1.22 ... .. ..$ percent.HB : num [1:2638] 0 0 0 0 0 0 0 0 0 0 ... 根據(jù)層級(jí)結(jié)構(gòu)找到 data 即可， pbmc@assays$RNA@data （標(biāo)準(zhǔn)化后的數(shù)據(jù)） 另：pbmc@assays$RNA@counts為原始count數(shù)據(jù) 簡(jiǎn)單看一下標(biāo)準(zhǔn)化前后的數(shù)據(jù)par(mfrow = c(1,2))hist(colSums(pbmc$RNA@counts),breaks = 50) hist(colSums(pbmc$RNA@data),breaks = 50) 2.3 高變基因（特征選擇）選擇數(shù)據(jù)集中高變異的特征子集（在某些細(xì)胞中高表達(dá)，在其他細(xì)胞中低表達(dá)）。通過Seurat內(nèi)置的FindVariableFeatures()函數(shù)，計(jì)算每一個(gè)基因的均值和方差，默認(rèn)選擇高變的2000個(gè)基因用于下游分析。 標(biāo)示Top10的基因 以及 標(biāo)示目標(biāo)基因 pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)# Identify the 10 most highly variable genes top10 <- head(VariableFeatures(pbmc), 10) # plot variable features with and without labels plot1 <- VariableFeaturePlot(pbmc) #標(biāo)示TOP10的基因 plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE) #標(biāo)示感興趣的基因 plot3 <- LabelPoints(plot = plot1, points = c(top10,"HLA-DPB1","CCL4"), repel = TRUE) plot2 + plot3 2.4 歸一化使用ScaleData()函數(shù)線性轉(zhuǎn)換（"歸一化"）數(shù)據(jù)，使得每一個(gè)基因在所有cell中的表達(dá)均值為0，方差為1. 結(jié)果在pbmc[["RNA"]]@scale.data #同樣可以通過str的方式根據(jù)數(shù)據(jù)結(jié)構(gòu)獲取。 all.genes <- rownames(pbmc)pbmc <- ScaleData(pbmc, features = all.genes) 有需要的話也可以移出不必要的來源數(shù)據(jù)？ pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")三 PCA 降維 print(pbmc[["pca"]], dims = 1:5, nfeatures = 5) PC_ 1 Positive: CST3, TYROBP, LST1, AIF1, FTL Negative: MALAT1, LTB, IL32, IL7R, CD2 PC_ 2 Positive: CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1 Negative: NKG7, PRF1, CST7, GZMA, GZMB PC_ 3 Positive: HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPA1 Negative: PPBP, PF4, SDPR, SPARC, GNG11 PC_ 4 Positive: HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1 Negative: VIM, IL7R, S100A6, S100A8, IL32 PC_ 5 Positive: GZMB, FGFBP2, S100A8, NKG7, GNLY Negative: LTB, IL7R, CKB, MS4A7, RP11-290F20.3 3.2 可視化展示 1）dims指定展示幾個(gè)PCs ， nfeatures指定每個(gè)PC展示多少基因 VizDimLoadings(pbmc, dims = 1:2,nfeatures = 20, reduction = "pca") 2） PCA點(diǎn)圖 DimPlot(pbmc, reduction = "pca")3） PCA熱圖 DimHeatmap(pbmc,dims = 1:4, #幾個(gè)PCs cells = 600, #多少cell ncol = 2, #圖形展示幾列 balanced = TRUE) 四 確定維度 如何確定使用多少PCA呢？以下2種方法科研作為參考 4.1 JackStrawPlot使用JackStrawPlot()函數(shù)對(duì)PCA結(jié)果進(jìn)行可視化，"重要"的PC 將會(huì)顯示在虛線上方且P值較低，本示例中PC10-PC12后顯著性下降較明顯。較耗時(shí)。 # NOTE: This process can take a long time for big datasets, comment out for expediency. More# approximate techniques such as those implemented in ElbowPlot() can be used to reduce # computation time pbmc <- JackStraw(pbmc, num.replicate = 100) pbmc <- ScoreJackStraw(pbmc, dims = 1:20) JackStrawPlot(pbmc, dims = 1:15) 4.2 Elbow圖展示每個(gè)主成分對(duì)數(shù)據(jù)方差的解釋情況（百分比表示），并進(jìn)行排序。發(fā)現(xiàn)第9個(gè)主成分是一個(gè)拐點(diǎn)，后續(xù)的主成分(PC)變化不大。 ElbowPlot(pbmc)A：以上結(jié)果“供參考”； B：Seurat官網(wǎng)鼓勵(lì)用戶使用不同數(shù)量的PC（10、15，甚至50）重復(fù)下游分析，雖然結(jié)果通常沒有顯著差異。 C：建議設(shè)置此參數(shù)時(shí)偏高一些，較少維度進(jìn)行下游分析可能會(huì)對(duì)結(jié)果產(chǎn)生一些負(fù)面影響。 4.3 顯著相關(guān)基因這個(gè)也可以作為后面分析選擇基因的一個(gè)參考。 #Returns a set of genes, based on the JackStraw analysis, that have statistically significant associations with a set of PCs.?PCASigGenes head(PCASigGenes(pbmc,pcs.use=2,pval.cut = 0.7)) [1] "PPBP" "LYZ" "S100A9" "IGLL5" "GNLY" "FTL" 五 聚類 基于PCA結(jié)果進(jìn)行聚類； 大約3K細(xì)胞的單細(xì)胞數(shù)據(jù)集，將resolution參數(shù)設(shè)置在0.4-1.2之間，數(shù)據(jù)集增加resolution 值對(duì)應(yīng)增加。這個(gè)參數(shù)的設(shè)置目前沒有找到標(biāo)準(zhǔn)，有清楚的小伙伴歡迎后臺(tái)告知，謝謝。 pbmc <- FindNeighbors(pbmc, dims = 1:10)pbmc <- FindClusters(pbmc, resolution = 0.5) # Look at cluster IDs of the first 5 cells head(Idents(pbmc), 5) head(pbmc@active.ident,5) #同上 AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1 AAACCGTGTATGCG-1 0 3 2 1 6 Levels: 0 1 2 3 4 5 6 7 8 5.1 查看指定cluster的cellhead(subset(as.data.frame(pbmc@active.ident),pbmc@active.ident=="2"))pbmc@active.identAAACATTGATCAGC-1 2AAACGCACTGGTAC-1 2 AAAGCCTGTATGCG-1 2 AAATCAACTCGCAA-1 2 AAATGTTGCCACAA-1 2 AACACGTGCAGAGG-1 2 5.2 提取某一cluster細(xì)胞subpbmc<-subset(x = pbmc,idents="2")subpbmc head(subpbmc@active.ident,5) AAACATTGATCAGC-1 AAACGCACTGGTAC-1 AAAGCCTGTATGCG-1 AAATCAACTCGCAA-1 AAATGTTGCCACAA-1 2 2 2 2 2 Levels: 2 六 非線性降維聚類 建議使用與聚類分析相同的pc維度進(jìn)行非線性降維度分析，如tSNE和UMAP，并可視化展示。 6.1 UMAP# If you haven't installed UMAP, you can do so via reticulate::py_install(packages =# 'umap-learn') pbmc <- RunUMAP(pbmc, dims = 1:10) # note that you can set `label = TRUE` or use the LabelClusters function to help label # individual clusters DimPlot(pbmc, reduction = "umap") 6.2 tSNEpbmc <- RunTSNE(pbmc, dims = 1:10)head(pbmc@reductions$tsne@cell.embeddings) DimPlot(pbmc, reduction = "tsne")
6.3 比較# note that you can set `label = TRUE` or use the LabelClusters function to help label# individual clusters plot1<-DimPlot(pbmc, reduction = "umap",label = TRUE) plot2<-DimPlot(pbmc, reduction = "tsne",label = TRUE) plot1 + plot2
可以看出，兩者的降維可視化的結(jié)構(gòu)是一致的，UMAP方法相對(duì)更加緊湊。 七 差異表達(dá)基因 Seurat可以通過FindMarkers函數(shù) 和 FindAllMarkers函數(shù)尋找不同cluster的差異表達(dá)基因。 min.pct參數(shù):設(shè)定在兩個(gè)細(xì)胞群中任何一個(gè)被檢測(cè)到的百分比，通過此設(shè)定不檢測(cè)很少表達(dá)基因來縮短程序運(yùn)行時(shí)間。默認(rèn)0.1 thresh.test參數(shù)：設(shè)定在兩個(gè)細(xì)胞群中基因差異表達(dá)量。可以設(shè)置為0 ，程序運(yùn)行時(shí)間會(huì)更長(zhǎng)。 max.cells.per.ident參數(shù)：每個(gè)類群細(xì)胞抽樣設(shè)置；也可以縮短程序運(yùn)行時(shí)間。 7.1 特定cluster是one-others的差異分析方法，由ident.1來制定cluster，本例就是cluster2與其余的cluster來做比較。 # find all markers of cluster 2cluster2.markers <- FindMarkers(pbmc, ident.1 = 2, min.pct = 0.25) head(cluster2.markers, n = 5)## p_val avg_log2FC pct.1 pct.2 p_val_adj ## IL32 2.593535e-91 1.2154360 0.949 0.466 3.556774e-87 ## LTB 7.994465e-87 1.2828597 0.981 0.644 1.096361e-82 ## CD3D 3.922451e-70 0.9359210 0.922 0.433 5.379250e-66 ## IL7R 1.130870e-66 1.1776027 0.748 0.327 1.550876e-62 ## LDHB 4.082189e-65 0.8837324 0.953 0.614 5.598314e-61 7.2 指定cluster本例為cluster5 和 cluster0_cluster3的差異 # find all markers distinguishing cluster 5 from clusters 0 and 3cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25) head(cluster5.markers, n = 5) p_val avg_log2FC pct.1 pct.2 p_val_adjFCGR3A 8.331882e-208 4.261784 0.975 0.040 1.142634e-203 CFD 1.932644e-198 3.423863 0.938 0.036 2.650429e-194 IFITM3 2.710023e-198 3.876058 0.975 0.049 3.716525e-194 CD68 1.069778e-193 3.013656 0.926 0.035 1.467094e-189 RP11-290F20.3 4.218926e-190 2.722303 0.840 0.016 5.785835e-186 7.3 FindAllMarkers每個(gè)cluster分別與其他所有cluster進(jìn)行比較，展示各個(gè)cluster的前2個(gè)基因 # find markers for every cluster compared to all remaining cells, report only the positive onespbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25) pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC) # A tibble: 18 x 7# Groups: cluster [9]p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr> 1 1.88e-117 1.08 0.913 0.588 2.57e-113 0 LDHB 2 5.01e- 85 1.34 0.437 0.108 6.88e- 81 0 CCR7 3 0 5.57 0.996 0.215 0 1 S100A9 4 0 5.48 0.975 0.121 0 1 S100A8 5 1.93e- 80 1.27 0.981 0.65 2.65e- 76 2 LTB 6 2.91e- 58 1.27 0.667 0.251 3.98e- 54 2 CD2 7 0 4.31 0.939 0.042 0 3 CD79A 8 1.06e-269 3.59 0.623 0.022 1.45e-265 3 TCL1A 9 3.60e-221 3.21 0.984 0.226 4.93e-217 4 CCL5 10 4.27e-176 3.01 0.573 0.05 5.85e-172 4 GZMK 11 3.51e-184 3.31 0.975 0.134 4.82e-180 5 FCGR3A 12 2.03e-125 3.09 1 0.315 2.78e-121 5 LST1 13 3.17e-267 4.83 0.961 0.068 4.35e-263 6 GZMB 14 1.04e-189 5.28 0.961 0.132 1.43e-185 6 GNLY 15 1.48e-220 3.87 0.812 0.011 2.03e-216 7 FCER1A 16 1.67e- 21 2.87 1 0.513 2.28e- 17 7 HLA-DPB1 17 7.73e-200 7.24 1 0.01 1.06e-195 8 PF4 18 3.68e-110 8.58 1 0.024 5.05e-106 8 PPBP ?FindAllMarkers查看更多參數(shù)。 cluster0.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25,test.use = "roc", # 檢驗(yàn)的方式 only.pos = TRUE) #只輸出pos的基因 其中test.use 可選參數(shù)有：wilcox（默認(rèn)），bimod，roc，t，negbinom，poisson，LR，MAST，DESeq2。 7.4 可視化可以通過seurat的內(nèi)置函數(shù)查看重點(diǎn)基因的基本情況： 1）VlnPlot: 基因表達(dá)概率分布 VlnPlot(pbmc, features = c("MS4A1", "CD79A"))# you can plot raw counts as well VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
"CD8A"))
3）DoHeatmap()查看重點(diǎn)基因細(xì)胞和cluster的表達(dá)熱圖 top10 <- pbmc.markers %>%group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC) DoHeatmap(pbmc, features = top10$gene) + NoLegend()
此外，還建議使用RidgePlot()、CellScatter()和DotPlot()查看數(shù)據(jù)集的情況。 這也可以作為后續(xù)手動(dòng)注釋的一些參考。 參考資料：https:///seurat/articles/pbmc3k_tutorial.html#standard-pre-processing-workflow-1 ◆ ◆ ◆ ◆ ◆ 精心整理（含圖PLUS版）|R語言生信分析，可視化（R統(tǒng)計(jì)，ggplot2繪圖，生信圖形可視化匯總）
|
|
|
