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

NGS系列文章包括NGS基础、转录组分析 （Nature重磅综述|关于RNA-seq你想知道的全在这）、ChIP-seq分析 （ChIP-seq基本分析流程）、单细胞测序分析 (重磅综述：三万字长文读懂单细胞RNA测序分析的最佳实践教程 （原理、代码和评述）)、DNA甲基化分析、重测序分析、GEO数据挖掘（典型医学设计实验GEO数据分析 (step-by-step) - Limma差异分析、火山图、功能富集）等内容。

背景介绍

肾脏是具有许多不同功能的高度复杂的器官，由几个功能和解剖上不连续的部分组成。肾小球和肾小管是肾单位的重要组成部分。足细胞与肾小球内皮细胞一起合成了肾小球基底膜，它是最终的过滤屏障，防止蛋白质损失到尿液中。顶叶上皮细胞（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个问题答疑精选！）

QC

# 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)

“harmony”整合不同平台的单细胞数据之旅

鉴定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")

• 单细胞分群后，怎么找到Marker基因定义每一类群？

可视化

#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"))　#　绘制部分基因热图

R语言 - 热图美化

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"))

• R语言 - 箱线图（小提琴图、抖动图、区域散点图）

• 可视化之为什么要使用箱线图？

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")

• Hemberg-lab单细胞转录组数据分析（十二）- Scater单细胞表达谱tSNE可视化

与前面的图是相同的。

提取PT cells进行后续分析

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