让你的单细胞数据动起来！|iCellR(一)
生活随笔
收集整理的這篇文章主要介紹了
让你的单细胞数据动起来！|iCellR(一)
小編覺得挺不錯的,現在分享給大家,幫大家做個參考.
今天在翻閱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個samplessample1 <- 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 03. 制作iCellR類的對象。
my.obj <- make.obj(my.data) my.obj ################################### ,--. ,-----. ,--.,--.,------. `--'' .--./ ,---. | || || .--. ' ,--.| | | .-. :| || || '--'.' | |' '--'\ --. | || || | `--' `-----' `----'`--'`--'`--' '--' ################################### An object of class iCellR version: 0.99.0Raw/original data dimentions (rows,columns): 32738,2700Data 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: FALSEClustering performed: FALSE , Number of clusters: 0tSNE performed: FALSE , UMAP performed: FALSE , DiffMap performed: FALSEMain data dimentions (rows,columns): 0 0Normalization 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 0ADT main data dimentions (rows,columns): 0 0ADT 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的表達分布（原則上遵循線性分布）R語言 - 散點圖繪制
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(my.obj@main.data) #[1] 32738 2637 # 過濾掉了63個細胞Hemberg-lab單細胞轉錄組數據分析（九）- Scater包單細胞過濾
Hemberg-lab單細胞轉錄組數據分析（十）- Scater基因評估和過濾
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 normalized9. 進行二次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")R語言 - 箱線圖（小提琴圖、抖動圖、區域散點圖）
10. Scale data
my.obj <- data.scale(my.obj) # 對數據進行中心化和標準化處理11.進行gene統計
my.obj <- gene.stats(my.obj, which.data = "main.data") head(my.obj@gene.data[order(my.obj@gene.data$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 all12.制作聚類的基因模型
在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(my.obj@gene.model) # 高變基因 # "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 = my.obj@gene.model,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(my.obj@gene.model) # 585 my.obj <- find.dim.genes(my.obj, dims = 1:10,top.pos = 20, top.neg = 20) # (optional) length(my.obj@gene.model) # 208 # second round PC my.obj <- run.pca(my.obj, method = "gene.model", gene.list = my.obj@gene.model,data.type = "main",batch.norm = F) my.obj@opt.pcsHemberg-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. 對細胞進行聚類分析
“聚類是針對給定的樣本，依據它們特征的相似度或距離，將其歸并到若干個“類”或“簇”的數據分析問題。一個類是樣本的一個子集。直觀上，相似的樣本聚集在相同的類，不相似的樣本分散在不同的類。這里，樣本之間的相似度或距離起著重要作用?！?br />降維是對feature (gene) 進行操作，而聚類是對sample進行操作，建議在得到高變基因以及降維后操作。(基因表達聚類分析之初探SOM)
建議使用以下默認選項：
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(my.obj@my.freq) # 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，下期接著介紹吧！
單細胞
收藏 北大生信平臺” 單細胞分析、染色質分析” 視頻和PPT分享
Science: 小鼠腎臟單細胞轉錄組+突變分析揭示腎病潛在的細胞靶標
Science：通過單細胞轉錄組測序揭示玉米減數分裂進程 | 很好的單細胞分析案例
Nature 首次對阿爾茨海默病進行單細胞轉錄組分析|詳細解讀
Cell 深度 一套普遍適用于各類單細胞測序數據集的錨定整合方案
骨髓基質在正常和白血病個體中的細胞圖譜 Cell,Nature聯袂解析
癌中之王：基質微環境塑造胰腺癌瘤內結構|Cell
Nature系列 整合單細胞轉錄組學和質譜流式確定類風濕性關節炎滑膜組織中的炎癥細胞狀態 詳細解讀
單細胞轉錄組教程匯總
10X單細胞測序分析軟件:Cell ranger，從拆庫到定量
Hemberg-lab單細胞轉錄組數據分析（一）- 引言
Hemberg-lab單細胞轉錄組數據分析（二）- 實驗平臺
Hemberg-lab單細胞轉錄組數據分析（三）- 原始數據質控
Hemberg-lab單細胞轉錄組數據分析（四）- 文庫拆分和細胞鑒定
Hemberg-lab單細胞轉錄組數據分析（五）- STAR, Kallisto定量
Hemberg-lab單細胞轉錄組數據分析（六）- 構建表達矩陣，UMI介紹
Hemberg-lab單細胞轉錄組數據分析（七）- 導入10X和SmartSeq2數據Tabula Muris
Hemberg-lab單細胞轉錄組數據分析（八）- Scater包輸入導入和存儲
Hemberg-lab單細胞轉錄組數據分析（九）- Scater包單細胞過濾
Hemberg-lab單細胞轉錄組數據分析（十）- Scater基因評估和過濾
Hemberg-lab單細胞轉錄組數據分析（十一）- Scater單細胞表達譜PCA可視化
Hemberg-lab單細胞轉錄組數據分析（十二）- Scater單細胞表達譜tSNE可視化
如何火眼金睛鑒定那些單細胞轉錄組中的混雜因素
什么？你做的差異基因方法不合適？
單細胞分群后，怎么找到Marker基因定義每一類群？
在線平臺如何做單細胞測序分析全套？有它so easy!
植物單細胞轉錄組的春天來了，還不上車？Science, PC, PP, MP, bioRxiv各一個
三人成虎，概率卻不足十分之五？
一文掌握GSEA，超詳細教程
這個只需一步就可做富集分析的網站還未發表就被CNS等引用超過350次
什么，你算出的P-value看上去像齊天大圣變的廟？
GO、GSEA富集分析一網打進
GSEA富集分析 - 界面操作
無需寫代碼的高顏值富集分析神器
去東方，最好用的在線GO富集分析工具
跨物種單細胞分析發現胰腺導管癌中一類有免疫原性的抗原呈遞成纖維細胞
NCB|心咽發育多樣化的單細胞轉錄軌跡分析
七龍珠|召喚一份單細胞數據庫匯總
用了這么多年的PCA可視化竟然是錯的！！！
單細胞預測Doublets軟件包匯總-過渡態細胞是真的嗎？
Seurat亮點之細胞周期評分和回歸
cellassign：用于腫瘤微環境分析的單細胞注釋工具（9月Nature）
Nature重磅綜述 |關于RNA-seq，你想知道的都在這
高顏值免費在線繪圖
往期精品(點擊圖片直達文字對應教程)
機器學習
后臺回復“生信寶典福利第一波”或點擊閱讀原文獲取教程合集
總結
以上是生活随笔為你收集整理的让你的单细胞数据动起来！|iCellR(一)的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 癌症精准医疗上市公司泛生子基因 - 内推
- 下一篇： Nat Commun｜单细胞ATAC-s