讓你的單細胞數據動起來！|iCellR(一)
- 2019 年 11 月 22 日
- 筆記
今天在翻閱single cell 的github時候，我看見了這個R包，允許我們處理各種來自單細胞測序技術的數據，如scRNA-seq，scVDJ-seq和CITE-Seq。
想看整套的學習流程還可以戳這裡：
https://vimeo.com/337822487
iCellR優點
Single（i）Cell R軟體包（iCellR）在分析pipeline的每個步驟提供前所未有的靈活性，包括標準化，聚類，降維，插補，可視化等。可以通過無監督和有監督聚類進行分析， 此外，該工具包提供2D和3D互動式可視化（一個震撼的交互型3D可視化R包 - 可直接轉ggplot2圖為3D），差異表達分析，基於細胞、基因和聚類的filter，數據合併，dropouts的標準化，數據插補方法，校正批次差異，找到標記基因的工具聚類和條件，預測細胞類型和偽時序分析（NBT|45種單細胞軌跡推斷方法比較，110個實際數據集和229個合成數據集）。感覺能做的真的好多啊！
安裝iCellR
# 從鏡像中進行安裝 install.packages("iCellR") # 也可以從github中進行安裝 #library(devtools) #install_github("rezakj/iCellR") # or #git clone https://github.com/rezakj/iCellR.git #R #install.packages('iCellR/', repos = NULL, type="source")下載演示數據
ion(samples = c("sample1","sample2","sample3"),# 設置需要下載的文件夾位置 setwd("/your/download/directory") # 通過URL進行數據下載 sample.file.url = "https://s3-us-west-2.amazonaws.com/10x.files/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz" # 下載文件 download.file(url = sample.file.url, destfile = "pbmc3k_filtered_gene_bc_matrices.tar.gz", method = "auto") # 解壓縮. untar("pbmc3k_filtered_gene_bc_matrices.tar.gz")更多的數據可以通過這裡獲得：https://genome.med.nyu.edu/results/external/iCellR/ 如何使用iCellR分析scRNA-seq數據 1. 載入iCellR包和下載的PBMC樣本數據。library("iCellR") my.data <- load10x("filtered_gene_bc_matrices/hg19/") ＃該目錄有barcodes.tsv，genes.tsv和matrix.mtx文件，將10x的數據讀入數據框 # 如果是SMART-seq2中有名的小鼠全器官數據，讀入方式如下 # my.data = read.delim("FACS/Kidney-counts.csv", sep=",", header=TRUE)2. 匯總數據dim(my.data) # [1] 32738 2700 # 將示例數據進行拆分為3個samples sample1 <- my.data[1:900] sample2 <- my.data[901:1800] sample3 <- my.data[1801:2700] # 合併示例文件my.data <- data.aggregat condition.names = c("WT","KO","Ctrl"))# 查看數據 head(my.data)[1:2] # WT_AAACATACAACCAC-1 WT_AAACATTGAGCTAC-1 #A1BG 0 0 #A1BG.AS1 0 0 #A1CF 0 0 #A2M 0 0 #A2M.AS1 0 0
3. 製作iCellR類的對象。
my.obj <- make.obj(my.data) my.obj ################################### ,--. ,-----. ,--.,--.,------. `--'' .--./ ,---. | || || .--. ' ,--.| | | .-. :| || || '--'.' | |' '--' --. | || || | `--' `-----' `----'`--'`--'`--' '--' ################################### An object of class iCellR version: 0.99.0 Raw/original data dimentions (rows,columns): 32738,2700 Data conditions in raw data: Ctrl,KO,WT (900,900,900) Row names: A1BG,A1BG.AS1,A1CF ... Columns names: WT_AAACATACAACCAC.1,WT_AAACATTGAGCTAC.1,WT_AAACATTGATCAGC.1 ... ################################### QC stats performed: FALSE , PCA performed: FALSE , CCA performed: FALSE Clustering performed: FALSE , Number of clusters: 0 tSNE performed: FALSE , UMAP performed: FALSE , DiffMap performed: FALSE Main data dimentions (rows,columns): 0 0 Normalization factors: ... Imputed data dimentions (rows,columns): 0 0 ############## scVDJ-Seq ########### VDJ data dimentions (rows,columns): 0 0 ############## CITE-Seq ############ ADT raw data dimentions (rows,columns): 0 0 ADT main data dimentions (rows,columns): 0 0 ADT columns names: ... ADT row names: ... ######## iCellR object made ########
4. 進行質控（QC）
my.obj <- qc.stats(my.obj) # 計算每個細胞的UMI數和基因數，還有細胞周期基因和線粒體基因（'mito.percent'）的百分比
5. 細胞周期預測
和上一個推送的細胞周期預測方法還是較為一致的，都是通過評分進行周期的預測，
my.obj <- cc(my.obj, s.genes = s.phase, g2m.genes = g2m.phase) head(my.obj@stats) # CellIds nGenes UMIs mito.percent #WT_AAACATACAACCAC.1 WT_AAACATACAACCAC.1 781 2421 0.030152829 #WT_AAACATTGAGCTAC.1 WT_AAACATTGAGCTAC.1 1352 4903 0.037935958 #WT_AAACATTGATCAGC.1 WT_AAACATTGATCAGC.1 1131 3149 0.008891712 #WT_AAACCGTGCTTCCG.1 WT_AAACCGTGCTTCCG.1 960 2639 0.017430845 #WT_AAACCGTGTATGCG.1 WT_AAACCGTGTATGCG.1 522 981 0.012232416 #WT_AAACGCACTGGTAC.1 WT_AAACGCACTGGTAC.1 782 2164 0.016635860 # S.phase.probability g2m.phase.probability S.Score #WT_AAACATACAACCAC.1 0.0012391574 0.0004130525 0.030569081 #WT_AAACATTGAGCTAC.1 0.0002039568 0.0004079135 -0.077860621 #WT_AAACATTGATCAGC.1 0.0003175611 0.0019053668 -0.028560560 #WT_AAACCGTGCTTCCG.1 0.0007578628 0.0011367942 0.001917225 #WT_AAACCGTGTATGCG.1 0.0000000000 0.0020387360 -0.020085210 #WT_AAACGCACTGGTAC.1 0.0000000000 0.0000000000 -0.038953135 # G2M.Score Phase #WT_AAACATACAACCAC.1 -0.0652390011 S #WT_AAACATTGAGCTAC.1 -0.1277015099 G1 #WT_AAACATTGATCAGC.1 -0.0036505733 G1 #WT_AAACCGTGCTTCCG.1 -0.0499511543 S #WT_AAACCGTGTATGCG.1 0.0009426363 G2M #WT_AAACGCACTGGTAC.1 -0.0680240629 G1 # 繪製細胞周期圖 pie(table(my.obj@stats$Phase))
6. Plot QC
默認情況下的繪圖功能都將創建互動式html文件，可以設置interactive = FALSE來關閉。
# 繪製 三個分組中UMIs, genes和 percent mito的beanplot分布圖，想單獨畫地話可以用`?stats.plot`help一下如何修改參數 stats.plot(my.obj, plot.type = "all.in.one", out.name = "UMI-plot", interactive = FALSE, cell.color = "slategray3", cell.size = 1, cell.transparency = 0.5, box.color = "red", box.line.col = "green")
# 繪製散點圖查看線粒體基因和基因分布，注意參數`plot.type` stats.plot(my.obj, plot.type = "point.mito.umi", out.name = "mito-umi-plot")###該圖為線粒體基因和UMI的相關表達分布（一般原則下會認為線粒體基因比例較高可能為破碎衰老細胞，可以設置cutoff值去掉，但如果研究為高代謝或本身耗能較多的細胞時需要考慮） stats.plot(my.obj, plot.type = "point.gene.umi", out.name = "gene-umi-plot")###此圖為UMI和GENE的表達分布（原則上遵循線性分布）
7. 過濾細胞
my.obj <- cell.filter(my.obj, min.mito = 0, max.mito = 0.05, min.genes = 200, max.genes = 2400, min.umis = 0, max.umis = Inf) #[1] "cells with min mito ratio of 0 and max mito ratio of 0.05 were filtered." 保留線粒體基因比例在0到0.05的細胞 #[1] "cells with min genes of 200 and max genes of 2400 were filtered." 保留基因數在200到2400之間的細胞 #[1] "No UMI number filter" 不通過UMI數和細胞ID來篩選 #[1] "No cell filter by provided gene/genes" #[1] "No cell id filter" #[1] "filters_set.txt file has beed generated and includes the filters set for this experiment."篩選之後會生成一個文件filter_set，含有實驗的篩選設置 # more examples 也可以通過基因和細胞ID來篩選 # my.obj <- cell.filter(my.obj, filter.by.gene = c("RPL13","RPL10")) # filter our cell having no counts for these genes # my.obj <- cell.filter(my.obj, filter.by.cell.id = c("WT_AAACATACAACCAC.1")) # filter our cell cell by their cell ids.# 查看剩下的細胞數(線粒體基因比例在0到0.05,基因數在200到2400之間) dim([email protected]) #[1] 32738 2637 # 過濾掉了63個細胞
8. 進行標準化
可以根據自己的研究對數據進行標準化。還可以使用iCellR以外的工具對數據進行標準化，並將數據導入iCellR。對於大多數單細胞研究，建議使用「ranked.glsf」標準化。這種標準化化比較適合具有大量零的矩陣。且在合併數據（操作參考步驟2：匯總數據）的情況下，也能對批次進行校正。
my.obj <- norm.data(my.obj, norm.method = "ranked.glsf", top.rank = 500) # best for scRNA-Seq # **更多的標準化方式** #my.obj <- norm.data(my.obj, norm.method = "ranked.deseq", top.rank = 500) #my.obj <- norm.data(my.obj, norm.method = "deseq") # best for bulk RNA-Seq #my.obj <- norm.data(my.obj, norm.method = "global.glsf") # best for bulk RNA-Seq #my.obj <- norm.data(my.obj, norm.method = "rpm", rpm.factor = 100000) # best for bulk RNA-Seq #my.obj <- norm.data(my.obj, norm.method = "spike.in", spike.in.factors = NULL) #my.obj <- norm.data(my.obj, norm.method = "no.norm") # if the data is already normalized
9. 進行二次QC
my.obj <- qc.stats(my.obj,which.data = "main.data") ### 該次QC的主要目的就是計算UMI數量，每個細胞的基因以及每個細胞和細胞周期基因的線粒體基因的百分比。從圖中可以看出其細胞周期相關基因表達較低，線粒體基因表達比例也較低 stats.plot(my.obj, plot.type = "all.in.one", out.name = "UMI-plot", interactive = F, cell.color = "slategray3", cell.size = 1, cell.transparency = 0.5, box.color = "red", box.line.col = "green", back.col = "white")
10. Scale data
my.obj <- data.scale(my.obj) # 對數據進行中心化和標準化處理
11.進行gene統計
my.obj <- gene.stats(my.obj, which.data = "main.data") head([email protected][order([email protected]$numberOfCells, decreasing = T),]) # 將數據通過擁有該基因的細胞數按從大到小進行排序 # genes numberOfCells totalNumberOfCells percentOfCells meanExp #30303 TMSB4X 2637 2637 100.00000 38.55948 #3633 B2M 2636 2637 99.96208 45.07327 #14403 MALAT1 2636 2637 99.96208 70.95452 #27191 RPL13A 2635 2637 99.92416 32.29009 #27185 RPL10 2632 2637 99.81039 35.43002 #27190 RPL13 2630 2637 99.73455 32.32106 # SDs condition #30303 7.545968e-15 all #3633 2.893940e+01 all #14403 7.996407e+01 all #27191 2.783799e+01 all #27185 2.599067e+01 all #27190 2.661361e+01 all
12.製作聚類的基因模型
在bulk RNA-seq數據中，用平均表達值meanExp排名前500的基因(設置參數base.mean.rank對樣本進行聚類是非常常見的，主要原因是為了降低雜訊。在scRNA-seq數據中也可以這樣做。並加上ranked.glsf歸一化，會較為擅長處理具有大量零的矩陣。
# 看一下高變基因集 make.gene.model(my.obj, my.out.put = "plot", dispersion.limit = 1.5, base.mean.rank = 500, no.mito.model = T, mark.mito = T, interactive = F, no.cell.cycle = T, out.name = "gene.model") # 將高變基因集注入`my.obj` 中 my.obj <- make.gene.model(my.obj, my.out.put = "data", dispersion.limit = 1.5, base.mean.rank = 500, no.mito.model = T, mark.mito = T, interactive = F, no.cell.cycle = T, out.name = "gene.model") head([email protected]) # 高變基因 # "ACTB" "ACTG1" "ACTR3" "AES" "AIF1" "ALDOA" # get html plot (optional) #make.gene.model(my.obj, my.out.put = "plot", # dispersion.limit = 1.5, # base.mean.rank = 500, # no.mito.model = T, # mark.mito = T, # interactive = T, # out.name = "plot4_gene.model")
13. 進行降維
降維是為了在資訊損失較少的情況下，以比原始特徵數量少得多的特徵（基因）來表示數據，更好地捕捉細胞之間的關係，常見降維的方法有PCA，t-SNE，umap和diffusion map。（Hemberg-lab單細胞轉錄組數據分析（十一）- Scater單細胞表達譜PCA可視化）
PCA分析
my.obj <- run.pca(my.obj, method = "gene.model", gene.list = [email protected],data.type = "main",batch.norm = F) opt.pcs.plot(my.obj) ##下圖為典型的滾石圖，當趨勢越趨於平緩表明其PC值作用越小，下圖PCA取到10。這一步是可選的，可能在一些情況下也是不適用的。# 2 round PCA (to find top genes in the first 10 PCs and re-run PCA for better clustering ## This is optional and might not be good in some cases length([email protected]) # 585 my.obj <- find.dim.genes(my.obj, dims = 1:10,top.pos = 20, top.neg = 20) # (optional) length([email protected]) # 208 # second round PC my.obj <- run.pca(my.obj, method = "gene.model", gene.list = [email protected],data.type = "main",batch.norm = F) [email protected]
Hemberg-lab單細胞轉錄組數據分析（十一）- Scater單細胞表達譜PCA可視化
更多降維方法
# tSNE my.obj <- run.pc.tsne(my.obj, dims = 1:10) # UMAP my.obj <- run.umap(my.obj, dims = 1:10, method = "naive") # or # my.obj <- run.umap(my.obj, dims = 1:10, method = "umap-learn") # diffusion map # this requires python packge phate # pip install --user phate # Install phateR version 2.9 # wget https://cran.r-project.org/src/contrib/Archive/phateR/phateR_0.2.9.tar.gz # install.packages('phateR/', repos = NULL, type="source") library(phateR) my.obj <- run.diffusion.map(my.obj, dims = 1:10, method = "phate")14. 對細胞進行聚類分析
「聚類是針對給定的樣本，依據它們特徵的相似度或距離，將其歸併到若干個「類」或「簇」的數據分析問題。一個類是樣本的一個子集。直觀上，相似的樣本聚集在相同的類，不相似的樣本分散在不同的類。這裡，樣本之間的相似度或距離起著重要作用。」 降維是對feature (gene) 進行操作，而聚類是對sample進行操作，建議在得到高變基因以及降維後操作。(基因表達聚類分析之初探SOM)
建議使用以下默認選項：
my.obj <- run.clustering(my.obj, clust.method = "kmeans", # 聚類方式為k-means dist.method = "euclidean", # 距離計算是歐氏距離，有非常明確的空間距離概念 index.method = "silhouette", max.clust = 25, min.clust = 2, dims = 1:10) # 如果要手動設置cluster數量，而不使用預測的最佳數量，請將最小值和最大值設置為所需的數量: #my.obj <- run.clustering(my.obj, # clust.method = "ward.D", # dist.method = "euclidean", # index.method = "ccc", # max.clust = 8, # min.clust = 8, # dims = 1:10) # more examples #my.obj <- run.clustering(my.obj, # clust.method = "ward.D", # dist.method = "euclidean", # index.method = "kl", # max.clust = 25, # min.clust = 2, # dims = 1:10)
15. 數據可視化
# clusters A= cluster.plot(my.obj,plot.type = "pca",interactive = F) B= cluster.plot(my.obj,plot.type = "umap",interactive = F) C= cluster.plot(my.obj,plot.type = "tsne",interactive = F) D= cluster.plot(my.obj,plot.type = "diffusion",interactive = F) #下圖通過不同降維方式進行展示，其實我們可以看出，在該組數據中還是t-sne降維後的可視化表達方式最好（細胞之間的分離度較好） library(gridExtra) grid.arrange(A,B,C,D) # conditions A= cluster.plot(my.obj,plot.type = "pca",col.by = "conditions",interactive = F) B= cluster.plot(my.obj,plot.type = "umap",col.by = "conditions",interactive = F) C= cluster.plot(my.obj,plot.type = "tsne",col.by = "conditions",interactive = F) D= cluster.plot(my.obj,plot.type = "diffusion",col.by = "conditions",interactive = F) #在下圖中我們可以發現無論哪種降維方式均沒有batch effect不同條件下的細胞重合度較好；library(gridExtra) grid.arrange(A,B,C,D)
16. 三維圖，密度圖和互動式圖
# 2D展示 cluster.plot(my.obj, cell.size = 1, plot.type = "tsne", cell.color = "black", back.col = "white", col.by = "clusters", cell.transparency = 0.5, clust.dim = 2, interactive = F) # 可交互的2D圖 cluster.plot(my.obj, plot.type = "tsne", col.by = "clusters", clust.dim = 2, interactive = T, out.name = "tSNE_2D_clusters") # 可交互的3D圖 cluster.plot(my.obj, plot.type = "tsne", col.by = "clusters", clust.dim = 3, interactive = T, out.name = "tSNE_3D_clusters") # cluster的密度圖 cluster.plot(my.obj, plot.type = "pca", col.by = "clusters", interactive = F, density=T) # Density plot for conditions cluster.plot(my.obj, plot.type = "pca", col.by = "conditions", interactive = F, density=T) #下圖為不同可視化方式，儘管方法多種多樣，但個人感覺還是2D的方式展示是最為合適的。
17. 不同cluster和條件下細胞比例
#如果normalize.ncell = TRUE，它將在三種狀態（conditions）隨機**降取樣（down sample）**，因此所有條件都具有相同數量的細胞，如果為FALSE則輸出原始細胞數。#下圖分別是在不同條件下細胞的比例關係，KO，WT，ctrl都是我們較為常見的細胞分組。# bar plot clust.cond.info(my.obj, plot.type = "bar", normalize.ncell = FALSE, my.out.put = "plot") # Pie chart clust.cond.info(my.obj, plot.type = "pie", normalize.ncell = FALSE, ,my.out.put = "plot") # data my.obj <- clust.cond.info(my.obj, plot.type = "bar", normalize.ncell = F) #head([email protected]) # conditions clusters Freq #1 ctrl 1 199 #2 KO 1 170 #3 WT 1 182 #4 ctrl 2 106 #5 KO 2 116 #6 WT 2 113
到這裡，我們其實已經表現了許多的plot，下期接著介紹吧
