NGS系列文章包括NGS基礎、轉錄組分析?（Nature重磅綜述|關于RNA-seq你想知道的全在這）、ChIP-seq分析?（ChIP-seq基本分析流程）、單細胞測序分析?(重磅綜述：三萬字長文讀懂單細胞RNA測序分析的最佳實踐教程 （原理、代碼和評述）)、DNA甲基化分析、重測序分析、GEO數(shù)據(jù)挖掘（典型醫(yī)學設計實驗GEO數(shù)據(jù)分析 (step-by-step) - Limma差異分析、火山圖、功能富集）等內容。

在文章Single-cell RNA sequencing of human kidney中，作者介紹了來自三個人類供體腎臟的23,366個高質量細胞的scRNA-seq數(shù)據(jù)，并將正常人腎細胞劃分為10個clusters。其中，近端腎小管（PT）細胞被分為三個亞型，而集合導管細胞被分為兩個亞型。該數(shù)據(jù)為腎細胞生物學和腎臟疾病的研究提供了可靠的參考。

在本周一的文章（復現(xiàn)原文（一）：Single-cell RNA sequencing of human kidney（step by step））中，我們完成了scRNA-seq數(shù)據(jù)的質控（Hemberg-lab單細胞轉錄組數(shù)據(jù)分析（三）- 原始數(shù)據(jù)質控）、批次校正、找Marker基因（單細胞分群后，怎么找到Marker基因定義每一類群？）、UMAP可視化和tSNE可視化（Hemberg-lab單細胞轉錄組數(shù)據(jù)分析（十二）- Scater單細胞表達譜tSNE可視化）。本期我們將提取近端腎小管（PT）細胞來完成下面的后續(xù)分析。

#Select a subset of PT cells PT <- subset(x = kid, idents= c("Proximal convoluted tubule cells","Proximal tubule cells","Proximal straight tubule cells")) saveRDS(PT,file="PT.rds")

數(shù)據(jù)準備

加載monocle軟件R包：

library(monocle)

monocle輸入文件類型有3種類型的格式：

• 表達量文件：exprs基因在所有細胞中的count值矩陣。

• 表型文件：phenoData。

• featureData

3個特征文件轉換成CellDataSet格式：

data <- as(as.matrix(PT@assays$RNA@counts), 'sparseMatrix') pd <- new('AnnotatedDataFrame', data = PT@meta.data) fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data)) fd <- new('AnnotatedDataFrame', data = fData) my_cds <- newCellDataSet(data,phenoData = pd,featureData = fd,lowerDetectionLimit = 0.5,expressionFamily = negbinomial.size())

對monocle對象進行歸一化：

my_cds <- estimateSizeFactors(my_cds) my_cds <- estimateDispersions(my_cds) my_cds <- detectGenes(my_cds, min_expr = 0.1) ##過濾掉低質量的基因。

查看數(shù)據(jù)：

head(fData(my_cds))

head(pData(my_cds))

pData(my_cds)$UMI <- Matrix::colSums(exprs(my_cds)) disp_table <- dispersionTable(my_cds) head(disp_table)

在進行降維聚類之前，先獲得高表達的基因集作為每個聚類用來order的Feature基因。也可以使用所有的基因，但一些表達量特別低的基因提供的聚類信號往往會制造分析噪音，Feature基因的選擇性很多，一種是可以根據(jù)基因的平均表達水平來進行篩選，另外我們也可以選擇細胞間異常變異的基因。這些基因往往能較好地反映不同細胞的狀態(tài)（對一篇單細胞RNA綜述的評述：細胞和基因質控參數(shù)的選擇）。

table(disp_table$mean_expression>=0.1) unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1)#細胞平均表達量大于0.1 my_cds <- setOrderingFilter(my_cds, unsup_clustering_genes$gene_id) plot_ordering_genes(my_cds)

plot_ordering_genes函數(shù)顯示了基因表達的變異性（分散性）是如何取決于整個細胞的平均表達的。

expressed_genes <- row.names(subset(fData(my_cds), num_cells_expressed >= 10)) #細胞表達個數(shù)大于10 my_cds_subset <- my_cds[expressed_genes, ] my_cds_subset head(pData(my_cds_subset)) my_cds_subset <- detectGenes(my_cds_subset, min_expr = 0.1) fData(my_cds_subset)$use_for_ordering <- fData(my_cds_subset)$num_cells_expressed > 0.05 * ncol(my_cds_subset) table(fData(my_cds_subset)$use_for_ordering) plot_pc_variance_explained(my_cds_subset, return_all = FALSE)

下面進行降維與聚類

my_cds_subset <- reduceDimension(my_cds_subset,max_components = 2,norm_method = 'log',num_dim = 10,reduction_method = 'tSNE',verbose = TRUE) my_cds_subset <- clusterCells(my_cds_subset, verbose = FALSE) plot_rho_delta(my_cds_subset, rho_threshold = 2, delta_threshold = 10) my_cds_subset <- clusterCells(my_cds_subset,rho_threshold = 2,delta_threshold = 10,skip_rho_sigma = T,verbose = FALSE) table(pData(my_cds_subset)$Cluster) plot_cell_clusters(my_cds_subset)

從上圖可以看出一種聚成9個clusters。（還在用PCA降維？快學學大牛最愛的t-SNE算法吧, 附Python/R代碼）

選擇定義細胞發(fā)展的基因

head(pData(my_cds_subset)) clustering_DEG_genes <- differentialGeneTest(my_cds_subset,fullModelFormulaStr = '~Cluster',cores = 22) dim(clustering_DEG_genes)

將細胞按照偽時間排序

library(dplyr) clustering_DEG_genes %>% arrange(qval) %>% head() my_ordering_genes <- row.names(clustering_DEG_genes)[order(clustering_DEG_genes$qval)][1:1000] my_cds_subset <- setOrderingFilter(my_cds_subset, ordering_genes = my_ordering_genes) my_cds_subset <- reduceDimension(my_cds_subset, method = 'DDRTree') #降維 my_cds_subset <- orderCells(my_cds_subset) #將細胞按照偽時間排序

偽時序軌跡按不同方面繪圖

注意參數(shù)color_by。（NBT|45種單細胞軌跡推斷方法比較，110個實際數(shù)據(jù)集和229個合成數(shù)據(jù)集）

plot_cell_trajectory(my_cds_subset, color_by = "State")

plot_cell_trajectory(my_cds_subset, color_by = "RNA_snn_res.0.25")

結果圖

原文圖

這里由于我們在前期分類的時候使用的是相同的參數(shù)，作者分為了3個cluster，我們分成了4個cluster，本質沒有什么太大的區(qū)別。

結果圖

原文圖

head(pData(my_cds_subset)) my_pseudotime_de <- differentialGeneTest(my_cds_subset,fullModelFormulaStr = "~sm.ns(Pseudotime)",cores = 22) my_pseudotime_de %>% arrange(qval) %>% head() my_pseudotime_de %>% arrange(qval) %>% head() %>% select(gene_short_name) -> my_pseudotime_gene plot_cell_trajectory(my_cds_subset, color_by = "Pseudotime")

結果圖

原文圖

影響fate decision的gene變化

#"A" stand for top 6 genes of affecting the fate decisions A=c("AKR1A1","PDZK1","AKR7A3","AKR7A2","FABP3","GADD45A") my_pseudotime_gene <-A plot_genes_in_pseudotime(my_cds_subset[my_pseudotime_gene,])

結果圖

原文圖

從上面我得到的圖和原文圖的比較來看，結果可能存在一定的差異（…

…），但基本模塊是相似的。

整個文章也并沒有對以上的偽時序分析進行描述，其最重要的部分應該還是作為人體腎臟單細胞的resource。大家也可以試試呀！

撰文：張虎
校對：生信寶典

參考文獻

Liao J, Yu Z, Chen Y, Bao M, Zou C, Zhang H, Liu D, Li T, Zhang Q, Li J, Cheng J, Mo Z. Single-cell RNA sequencing of human kidney. Sci Data. 2020 Jan 2;7(1):4.doi: 10.1038

總結

以上是生活随笔為你收集整理的复现原文（二）：Single-cell RNA sequencing of human的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內容還不錯，歡迎將生活随笔推薦給好友。