https://www.nature.com/articles/s41597-019-0351-8
NGS系列文章包括NGS基礎(chǔ)狼钮、轉(zhuǎn)錄組分析 （Nature重磅綜述|關(guān)于RNA-seq你想知道的全在這）约急、ChIP-seq分析 （ChIP-seq基本分析流程）标沪、單細(xì)胞測(cè)序分析 (重磅綜述：三萬字長(zhǎng)文讀懂單細(xì)胞RNA測(cè)序分析的最佳實(shí)踐教程 （原理乌妙、代碼和評(píng)述）)、DNA甲基化分析腐魂、重測(cè)序分析谍咆、GEO數(shù)據(jù)挖掘（典型醫(yī)學(xué)設(shè)計(jì)實(shí)驗(yàn)GEO數(shù)據(jù)分析 (step-by-step) - Limma差異分析纵揍、火山圖、功能富集）等內(nèi)容垛玻。
背景介紹
腎臟是具有許多不同功能的高度復(fù)雜的器官割捅，由幾個(gè)功能和解剖上不連續(xù)的部分組成。腎小球和腎小管是腎單位的重要組成部分帚桩。足細(xì)胞與腎小球內(nèi)皮細(xì)胞一起合成了腎小球基底膜亿驾，它是最終的過濾屏障，防止蛋白質(zhì)損失到尿液中账嚎。頂葉上皮細(xì)胞（Parietal epithelial cells莫瞬，PECs）是另一種常見的腎小球細(xì)胞類型参淹，可能導(dǎo)致腎小球硬化、新月和假新月形成乏悄。近端小管（proximal tubule浙值，PT）通過控制Na+ - H+和HCO3-的轉(zhuǎn)運(yùn)在調(diào)節(jié)全身酸堿平衡中起著重要作用，而遠(yuǎn)曲小管則更多地參與電解質(zhì)的轉(zhuǎn)運(yùn)檩小。在先前的研究中开呐，研究人員對(duì)腎臟不同組成部分進(jìn)行了bulk RNA測(cè)序（RNA-seq最強(qiáng)綜述名詞解釋&思維導(dǎo)圖|關(guān)于RNA-seq，你想知道的都在這（續(xù)））规求，為理解不同片段的轉(zhuǎn)錄組提供參考筐付。然而，bulk RNA測(cè)序不能反映單細(xì)胞水平的轉(zhuǎn)錄組阻肿，只能反映總體平均RNA表達(dá)（自從用了這個(gè)神器瓦戚，大規(guī)模RNA-seq數(shù)據(jù)挖掘我也可以）。
正常人腎臟的全面細(xì)胞解剖結(jié)構(gòu)對(duì)于解決腎臟疾病和腎癌的細(xì)胞起源至關(guān)重要丛塌。一些腎臟疾病可能是細(xì)胞類型特異性的较解，尤其是腎小管細(xì)胞。為了研究人腎臟的分類和轉(zhuǎn)錄組信息赴邻，作者迅速獲得了腎臟的單細(xì)胞懸液并進(jìn)行了單細(xì)胞RNA測(cè)序（scRNA-seq）（重磅綜述：三萬字長(zhǎng)文讀懂單細(xì)胞RNA測(cè)序分析的最佳實(shí)踐教程 （原理印衔、代碼和評(píng)述））。作者介紹了來自三個(gè)人類供體腎臟的23,366個(gè)高質(zhì)量細(xì)胞的scRNA-seq數(shù)據(jù)姥敛，并將正常人腎細(xì)胞劃分為10個(gè)clusters奸焙。其中，近端腎小管（PT）細(xì)胞被分為三個(gè)亞型彤敛，而集合導(dǎo)管細(xì)胞被分為兩個(gè)亞型与帆。總體而言墨榄，該數(shù)據(jù)為腎細(xì)胞生物學(xué)和腎臟疾病的研究提供了可靠的參考玄糟。
下面我們按照作者的分析思路復(fù)現(xiàn)該文章的部分內(nèi)容：
首先，從GSE131685中下載數(shù)據(jù)：
里面的文件名要分別改為“barcodes.tsv”渠概、“genes.tsv”和“matrix.mtx”茶凳，在Read10X（Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（七）- 導(dǎo)入10X和SmartSeq2數(shù)據(jù)Tabula Muris）時(shí)才不會(huì)報(bào)錯(cuò)。播揪。贮喧。
Kidney data loading
library(devtools)
install_github("immunogenomics/harmony")
library(Seurat)
library(magrittr)
library(harmony)
library(dplyr)
#Kidney data loading 并構(gòu)建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函數(shù)進(jìn)行合并
插一句嘴，我們來看一下華盛頓大學(xué)PhD jared.andrews對(duì)merge函數(shù)的解釋：
注意老鐵說的“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”（單細(xì)胞分析Seurat使用相關(guān)的10個(gè)問題答疑精選猪狈！）
QC
# quality control
kid[["percent.mt"]] <- PercentageFeatureSet(kid, pattern = "^MT-") #提取有關(guān)線粒體的基因
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) #標(biāo)準(zhǔn)化
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))
# 計(jì)算細(xì)胞周期
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)
這里我想叨叨幾句箱沦，據(jù)我看到的文獻(xiàn)，多數(shù)是在進(jìn)行降維后將細(xì)胞周期方面對(duì)分群的影響作為一個(gè)單獨(dú)模塊去敘述雇庙，作者在先期不管細(xì)胞周期對(duì)聚類是否有影響的情況下就對(duì)細(xì)胞周期相關(guān)基因進(jìn)行去除也是比較明智的谓形，因?yàn)樽髡卟⒉幌胱屧撘蛩鼗祀s其中影響分群（如何火眼金睛鑒定那些單細(xì)胞轉(zhuǎn)錄組中的混雜因素）灶伊。
#當(dāng)然我們還是要看是否細(xì)胞周期真的有影響，感興趣的小伙伴可以看一下寒跳，確實(shí)是有一定影響的聘萨！#kid <- ScaleData(kid, features = rownames(kid))
#kid <- RunPCA(kid , features = c(s.genes, g2m.genes))
#DimPlot(kid)
利用harmony算法去除批次效應(yīng)并細(xì)胞分類
#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) #等候時(shí)間較長(zhǎng)，請(qǐng)溜達(dá)溜達(dá)吧
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”整合不同平臺(tái)的單細(xì)胞數(shù)據(jù)之旅
鑒定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) #按照細(xì)胞周期進(jìn)行劃分
DimPlot(kid, reduction = "umap", label = TRUE, pt.size = .1) #注意作者在用同樣參數(shù)設(shè)置后分為10個(gè)clusters，其實(shí)無關(guān)緊要书释，都需要通過marker重新貼現(xiàn)翘贮。
根據(jù)作者提供的marker對(duì)細(xì)胞亞群進(jìn)行貼現(xiàn)，如下圖所示：
其實(shí)部分marker并不是特異性marker爆惧，所以在進(jìn)行區(qū)分的時(shí)候一定要好好甄別狸页。
與以下原文圖基本相同，個(gè)人感覺tSNE是不是也有什么隨機(jī)種子的東東扯再，感覺總會(huì)略有不同：
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進(jìn)行后續(xù)分析
#Select a subset of PT cells（近端小管）
PT <- SubsetData(kid, ident.use = c(0,1,2), subset.raw = T)
