（未测试）复现原文（一）：Single-cell RNA sequencing of human kidney（step by step）

https://www.nature.com/articles/s41597-019-0351-8

背景介绍
肾脏是具有许多不同功能的高度复杂的器官，由几个功能和解剖上不连续的部分组成。肾小球和肾小管是肾单位的重要组成部分。足细胞与肾小球内皮细胞一起合成了肾小球基底膜，它是最终的过滤屏障，防止蛋白质损失到尿液中。顶叶上皮细胞（Parietal epithelial cells，PECs）是另一种常见的肾小球细胞类型，可能导致肾小球硬化、新月和假新月形成。近端小管（proximal tubule，PT）通过控制Na+ – H+和HCO3-的转运在调节全身酸碱平衡中起着重要作用，而远曲小管则更多地参与电解质的转运。在先前的研究中，研究人员对肾脏不同组成部分进行了bulk RNA测序（RNA-seq最强综述名词解释&思维导图|关于RNA-seq，你想知道的都在这（续）），为理解不同片段的转录组提供参考。然而，bulk RNA测序不能反映单细胞水平的转录组，只能反映总体平均RNA表达（自从用了这个神器，大规模RNA-seq数据挖掘我也可以）。

正常人肾脏的全面细胞解剖结构对于解决肾脏疾病和肾癌的细胞起源至关重要。一些肾脏疾病可能是细胞类型特异性的，尤其是肾小管细胞。为了研究人肾脏的分类和转录组信息，作者迅速获得了肾脏的单细胞悬液并进行了单细胞RNA测序（scRNA-seq）（重磅综述：三万字长文读懂单细胞RNA测序分析的最佳实践教程（原理、代码和评述））。作者介绍了来自三个人类供体肾脏的23,366个高质量细胞的scRNA-seq数据，并将正常人肾细胞划分为10个clusters。其中，近端肾小管（PT）细胞被分为三个亚型，而集合导管细胞被分为两个亚型。总体而言，该数据为肾细胞生物学和肾脏疾病的研究提供了可靠的参考。

下面我们按照作者的分析思路复现该文章的部分内容：

首先，从GSE131685中下载数据：

里面的文件名要分别改为“barcodes.tsv”、“genes.tsv”和“matrix.mtx”，在Read10X（Hemberg-lab单细胞转录组数据分析（七）- 导入10X和SmartSeq2数据Tabula Muris）时才不会报错。。。

Kidney data loading

library(devtools)
install_github("immunogenomics/harmony")
library(Seurat)
library(magrittr)
library(harmony)
library(dplyr)

#Kidney data loading 并构建seurat object

K1.data <- Read10X(data.dir = "/Users/zhanghu1992/Documents/GSE131685_RAW/kidney1/")
K1 <- CreateSeuratObject(counts = K1.data, project = "kidney1", min.cells = 8, min.features = 200)
K2.data <- Read10X(data.dir = "/Users/zhanghu1992/Documents/GSE131685_RAW/kidney2/")
K2 <- CreateSeuratObject(counts = K2.data, project = "kidney2", min.cells = 6, min.features = 200)
K3.data <- Read10X(data.dir = "/Users/zhanghu1992/Documents/GSE131685_RAW/kidney3/")
K3 <- CreateSeuratObject(counts = K3.data, project = "kidney3", min.cells = 10, min.features = 200)
kid <- merge(x = K1, y = list(K2, K3)) #读取文件并用merge函数进行合并

插一句嘴，我们来看一下华盛顿大学PhD jared.andrews对merge函数的解释：

注意老铁说的“Seurat’s integration method is quite heavy handed in my experience,so if you decide to go the integration route,I’d recommend using the SeuratWrapper around the fastMNN”（单细胞分析Seurat使用相关的10个问题答疑精选！）

# quality control
kid[["percent.mt"]] <- PercentageFeatureSet(kid, pattern = "^MT-") #提取有关线粒体的基因
VlnPlot(kid, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3) #由图可以看出分布还可以

plot1 <- FeatureScatter(kid, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(kid, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))

kid <- subset(kid, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 30) #筛选条件
kid <- NormalizeData(kid, normalization.method = "LogNormalize", scale.factor = 10000)
kid <- NormalizeData(kid) #标准化
kid <- FindVariableFeatures(kid, selection.method = "vst", nfeatures = 2000) #查找高变基因
top10 <- head(VariableFeatures(kid), 10)
plot1 <- VariableFeaturePlot(kid)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
CombinePlots(plots = list(plot1, plot2))

# 计算细胞周期
s.genes <-cc.genes$s.genes
g2m.genes<-cc.genes$g2m.genes
kid <- CellCycleScoring(kid, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)
all.genes <- rownames(kid)
kid <- ScaleData(kid, vars.to.regress = c("S.Score", "G2M.Score"), features = all.genes)

这里我想叨叨几句，据我看到的文献，多数是在进行降维后将细胞周期方面对分群的影响作为一个单独模块去叙述，作者在先期不管细胞周期对聚类是否有影响的情况下就对细胞周期相关基因进行去除也是比较明智的，因为作者并不想让该因素混杂其中影响分群（如何火眼金睛鉴定那些单细胞转录组中的混杂因素）。

#当然我们还是要看是否细胞周期真的有影响，感兴趣的小伙伴可以看一下，确实是有一定影响的！#kid <- ScaleData(kid, features = rownames(kid))
#kid  <- RunPCA(kid , features = c(s.genes, g2m.genes))
#DimPlot(kid)

利用harmony算法去除批次效应并细胞分类

#Eliminate batch effects with harmony and cell classification
kid <- RunPCA(kid, pc.genes = kid@var.genes, npcs = 20, verbose = FALSE)
options(repr.plot.height = 2.5, repr.plot.width = 6)
kid <- kid %>%
    RunHarmony("orig.ident", plot_convergence = TRUE)　#等候时间较长，请溜达溜达吧
harmony_embeddings <- Embeddings(kid, 'harmony')
harmony_embeddings[1:5, 1:5]
kid <- kid %>%
    RunUMAP(reduction = "harmony", dims = 1:20) %>%
    FindNeighbors(reduction = "harmony", dims = 1:20) %>%
    FindClusters(resolution = 0.25) %>%
    identity()
new.cluster.ids <- c(0,1, 2, 3, 4, 5, 6, 7,8,9,10)
names(new.cluster.ids) <- levels(kid)
kid <- RenameIdents(kid, new.cluster.ids)

鉴定Marker基因

#Calculating differentially expressed genes (DEGs) and Save rds file
kid.markers <- FindAllMarkers(kid, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25)#寻找高变基因
write.table(kid.markers,sep="\t",file="/home/yuzhenyuan/Seurat/0.2_20.xls")
saveRDS(kid,file="/home/yuzhenyuan/kid/har/0.25_20.rds")

可视化

#Some visual figure generation
DimPlot(kid, reduction = "umap", group.by = "orig.ident", pt.size = .1, split.by = 'orig.ident')

DimPlot(kid, reduction = “umap”, group.by = “Phase”, pt.size = .1)　#按照细胞周期进行划分

DimPlot(kid, reduction = "umap", label = TRUE, pt.size = .1)　#注意作者在用同样参数设置后分为10个clusters，其实无关紧要，都需要通过marker重新贴现。

根据作者提供的marker对细胞亚群进行贴现，如下图所示：

其实部分marker并不是特异性marker，所以在进行区分的时候一定要好好甄别。

与以下原文图基本相同，个人感觉tSNE是不是也有什么随机种子的东东，感觉总会略有不同：

DoHeatmap(kid, features = c("SLC13A3","SLC34A1","GPX3","DCXR","SLC17A3","SLC22A8","SLC22A7","GNLY","NKG7","CD3D","CD3E","LYZ","CD14","KRT8","KRT18","CD24","VCAM1","UMOD","DEFB1","CLDN8","AQP2","CD79A","CD79B","ATP6V1G3","ATP6V0D2","TMEM213"))　#　绘制部分基因热图

VlnPlot(kid, pt.size =0, idents= c(1,2,3), features = c(“GPX3”, “DCXR”,”SLC13A3″,”SLC34A1″,”SLC22A8″,”SLC22A7″))
VlnPlot(kid, idents= c(8,10), features = c(“AQP2”, “ATP6V1B1″,”ATP6V0D2″,”ATP6V1G3”))

tSNE plot

##tSNE Plot
kid <-RunTSNE(kid, reduction = "harmony", dims = 1:20)
TSNEPlot(kid, do.label = T, label = TRUE, do.return = T, pt.size = 1)
TSNEPlot(kid, do.return = T, pt.size = 1, group.by = "orig.ident", split.by = 'orig.ident')
TSNEPlot(kid, do.return = T, pt.size = 1, group.by = "Phase")

与前面的图是相同的。

提取PT cells进行后续分析

#Select a subset of PT cells（近端小管）
PT <- SubsetData(kid, ident.use = c(0,1,2), subset.raw = T)

在文章Single-cell RNA sequencing of human kidney中，作者介绍了来自三个人类供体肾脏的23,366个高质量细胞的scRNA-seq数据，并将正常人肾细胞划分为10个clusters。其中，近端肾小管（PT）细胞被分为三个亚型，而集合导管细胞被分为两个亚型。该数据为肾细胞生物学和肾脏疾病的研究提供了可靠的参考。

在本周一的文章（复现原文（一）：Single-cell RNA sequencing of human kidney（step by step））中，我们完成了scRNA-seq数据的质控（Hemberg-lab单细胞转录组数据分析（三）- 原始数据质控）、批次校正、找Marker基因（单细胞分群后，怎么找到Marker基因定义每一类群？）、UMAP可视化和tSNE可视化（Hemberg-lab单细胞转录组数据分析（十二）- Scater单细胞表达谱tSNE可视化）。本期我们将提取近端肾小管（PT）细胞来完成下面的后续分析。

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

数据准备

加载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) ##过滤掉低质量的基因。

查看数据：

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基因的选择性很多，一种是可以根据基因的平均表达水平来进行筛选，另外我们也可以选择细胞间异常变异的基因。这些基因往往能较好地反映不同细胞的状态（对一篇单细胞RNA综述的评述：细胞和基因质控参数的选择）。

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函数显示了基因表达的变异性（分散性）是如何取决于整个细胞的平均表达的。

expressed_genes <- row.names(subset(fData(my_cds), num_cells_expressed >= 10)) #细胞表达个数大于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代码）

选择定义细胞发展的基因

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) #将细胞按照伪时间排序

伪时序轨迹按不同方面绘图

注意参数color_by。（NBT|45种单细胞轨迹推断方法比较，110个实际数据集和229个合成数据集）

plot_cell_trajectory(my_cds_subset, color_by = "State")

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

这里由于我们在前期分类的时候使用的是相同的参数，作者分为了3个cluster，我们分成了4个cluster，本质没有什么太大的区别。

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,])