（單細胞-SingleCell）Seurat流程文獻複現——單細胞實戰分析流程

單細胞項目：來自于30個病人的49個組織樣品，跨越3個治療階段

Therapy-Induced Evolution of Human Lung Cancer Revealed by Single-Cell RNA Sequencing

這篇教程我将分為四個階段完整的闡述單細胞的主流的下遊分析流程

1. 資料預處理（資料準備階段）
2. seurat 基礎分析
3. 免疫細胞識别
4. inferCNV 的實作

先來第一部分資料預處理

這裡的資料預處理很自由，其中一部分是必須按照嚴格的标準進行比如計算外源基因，線粒體， 核糖體等等

# 導入包
library(Seurat)
library(dplyr)
library(Matrix)
library(magrittr)
library(stringr)
library(tidyverse)
options(stringsAsFactors = F)
options(as.is = T)
setwd('/Users/yuansh/Desktop/單細胞/scell_lung_adenocarcinoma-master')
           

文章的話一共是有 2 個表達矩陣需要處理

我們這裡的話要很明确的知道表達譜的預處理到底預處理什麼東西

1. 看一下 count 矩陣是否是 hista 比對結果，如果是需要剔除 hista 比對的注釋資訊。（有的也是已經處理完了，但是也得檢查）
2. 确定一下 cell 名的組成，建構好 metadata（行名） 和 count（列名） 矩陣能夠對應
3. 過濾掉外源 RNA （ERCC）
4. 合并成為 seurat 對象
5. 計算線粒體（MT），核糖體的 RNA（^RP[SL]） 比例并篩選資料

我們先處理第一個表達矩陣

# 這個是治療前的第一個表達矩陣
# 導入資料
raw.data = read.csv('csv_files/S01_datafinal.csv', header=T, row.names = 1)
metadata = read.csv("csv_files/S01_metacells.csv", row.names=1, header=T)
dim(raw.data)
dim(metadata)

           

判斷一下是否是 hista 比對的結果并且沒有處理，如果是的話後面 5 行是沒用的

然後我們發現并不是

繼續觀察發現,raw.data 是 count 矩陣,metadata 是資訊矩陣

count 矩陣的列名是 meta 矩陣的 well 和 plate 合并的

稍微對資料進行處理一下,然後儲存，確定 count 矩陣的列名和 metadata 的行名一樣

這樣的話下次讀取就很快

> tail(raw.data[,1:4])
A10_1001000329 A10_1001000362 A10_1001000366 A10_1001000367
ZXDC                0              0              0              0
ZYG11A              0              0              0              0
ZYG11B              0              0              0              0
ZYX               292              0              0              0
ZZEF1              66              0              0              0
ZZZ3                0              0              0              0

> raw.data[1:4,1:4]
A10_1001000329 A10_1001000362 A10_1001000366 A10_1001000367
A1BG                  3              0              0              0
A1BG-AS1              0              0              0              0
A1CF                  0              0              0              0
A2M                  68              0              0              0

> metadata[1:4,1:4]
well      plate        cell_id sample_name
1  A10 1001000329 A10_1001000329      LT_S07
2  A10 1001000362 A10_1001000362      LT_S12
3  A10 1001000366 A10_1001000366      LT_S13
4  A10 1001000367 A10_1001000367      LT_S13
           

paste0(metadata$well[1],'_',metadata$plate[1]) # 在進行大幅度修改的時候要判斷代碼有沒有寫錯
colnames(raw.data)[1]
# 修改一下 metadata 的行名
rownames(metadata) = paste0(metadata$well,'_',metadata$plate)
# 确定了行名metadata 的行名和 count 矩陣 raw.data 的列名相等
identical(colnames(raw.data),rownames(metadata))
           

接下來就是基因的處理，

按照上面說的步驟：

1. 剔除外源基因，然後建構 seurat 對象
2. 計算線粒體基因（有沒有都算就行了）
3. 計算核糖體基因

# 這裡要過濾 ERCC(External RNA Control Consortium)
# ERCC 是外源 RNA 主要的作用就是用來質控的
erccs = grep('^ERCC-', x= rownames(x=raw.data),value = T) # value = T 擷取名字
# 計算 ercc 比例
percent.ercc = Matrix::colSums(raw.data[erccs, ])/Matrix::colSums(raw.data)
fivenum(percent.ercc) 
ercc.index = grep(pattern = "^ERCC-", x = rownames(x = raw.data), value = FALSE) # 擷取 index
# 删除 Count 矩陣中的 ERCC
raw.data = raw.data[-ercc.index,]
dim(raw.data) 
# [1] 26485 27489
           

建構 seurat 對象

# 建構對象
# min.cells 某一個基因至少在多少個基因中表達
# min.features 某個細胞至少表達多少個基因
min.cells = 0 # 3 因為文章沒有篩選，是以這裡設為0
min.features = 0 # 150 # 同上
sce = CreateSeuratObject(counts = raw.data,metadata = metadata,min.cells =min.cells,min.features =min.features)
sce
           

我們先看一下 sce 對象中的 meta.data,這個是建構 seurat 對象自動生成的，缺少了很多資訊，是以要把 metdata 添加上去

> [email protected][1:5,1:5]
orig.ident nCount_RNA nFeature_RNA sample_name patient_id
A10_1001000329 SeuratProject     681122         2996      LT_S07      TH103
A10_1001000362 SeuratProject         30           17      LT_S12      TH172
A10_1001000366 SeuratProject         29           19      LT_S13      TH169
A10_1001000367 SeuratProject        286            3      LT_S13      TH169
A10_1001000372 SeuratProject        599          240      LT_S14      TH153
           

sce = AddMetaData(object = sce, metadata = metadata)
sce = AddMetaData(object = sce, percent.ercc, col.name = "percent.ercc")
# Head to check
head([email protected])
           

然後就是資料清洗

資料清洗原則：

1. 不符合品質的：這裡一般指細胞中基因數太少或者某些基因隻在兩三個細胞中表達
2. 過濾線粒體 DNA占比過高的

3. 過濾核糖體占比過高的

   在過濾前先計算

table(grepl("^MT-",rownames(sce)))
table(grepl("^RP[SL][[:digit:]]",rownames(sce)))
sce[["percent.mt"]]  = PercentageFeatureSet(sce, pattern = "^MT-")
sce[["percent.rp"]]  = PercentageFeatureSet(sce, pattern = "^RP[SL][[:digit:]]")
dim(sce)
summary(sce[["nCount_RNA"]])
summary(sce[["nFeature_RNA"]])
summary(sce[["percent.mt"]]  )
summary(sce[["percent.ercc"]] )
summary(sce[["percent.rp"]]  )
           

文章中是啥都沒去，就正常的過濾了一下基因和表達過低的細胞

# 過濾篩選
# 一般根據實際情況進行篩選，這次的篩選主要還是根據文章裡面的流程
# 這個過濾蠻自由的，是以随便搞一下就行
myData1 = subset(x=sce, subset = nCount_RNA > 50000 & nFeature_RNA > 500)
myData1
save(myData1,file='myData1.Rdata')
           

第一階段的處理流程就是這樣

接着處理第二個表達矩陣

# 這個是治療前的第一個表達矩陣
# 導入資料
raw.data = read.csv('Data_input/csv_files/neo-osi_rawdata.csv', header=T, row.names = 1)
metadata = read.csv("Data_input/csv_files/neo-osi_metadata.csv", row.names=1, header=T)
raw.data[1:4,1:4]
metadata[1:4,1:4]
           

判斷一下是否是 hista 比對的結果并且沒有處理，如果是的話後面 5 行是沒用的

然後發現果然是

繼續觀察發現,第二套資料的命名規則和第一套有點不太一樣,多了_S10.homo,是以要把他去掉

稍微對資料進行處理一下,然後儲存，確定 count 矩陣的列名和 metadata 的行名一樣

這樣的話下次讀取就很快

> tail(raw.data[,1:4])
A10_B000561_S10.homo A10_B000563_S10.homo A10_B000568_S10.homo A10_B001544_S10.homo
ZZZ3                                      0                    0                    0                    0
__no_feature                          31084                92098               131276              1197071
__ambiguous                             353                 2231                  570                21247
__too_low_aQual                           0                    0                    0                    0
__not_aligned                             0                    0                    0                    0
__alignment_not_unique                67598                44778                47959               500993
> raw.data[1:4,1:4]
A10_B000561_S10.homo A10_B000563_S10.homo A10_B000568_S10.homo A10_B001544_S10.homo
A1BG                        0                    0                    0                    0
A1BG-AS1                    0                    0                    0                    0
A1CF                        0                    0                    0                    0
A2M                         0                    0                    0                    0
> metadata[1:4,1:4]
plate sample_name patient_id sample_type
1 B001567       AZ_01      AZ003       tumor
2 B001544       AZ_01      AZ003       tumor
3 B001546       AZ_01      AZ003       tumor
4 B001481       AZ_01      AZ003       tumor
           

然後我們進一步的發現，這套資料好像是沒有 well 資訊的

是以又要進一步處理一下，確定count 矩陣的列名和 metadata 矩陣的行名可以對應

gsub('_S[[:digit:]]*.homo','',colnames(raw.data))[1:5]
colnames(raw.data) = gsub('_S[[:digit:]]*.homo','',colnames(raw.data))
raw.data = raw.data[-grep('__',rownames(raw.data)),]
tail(raw.data[,1:4])
metadata$well[1] # 發現這套資料中的 metadata 中沒有 well 資訊
metadata$plate[1] # 好在有 plate資訊
           

是以稍微進行處理一下

# 将 count 矩陣的名字拆開自己建構 well 資訊矩陣
library(stringr)
wellInfo = str_split(colnames(raw.data),'_',simplify = T)%>% as.data.frame()
colnames(wellInfo) = c("well", "plate")
rownames(wellInfo) = paste0(wellInfo[,1],'_',wellInfo[,2]) 
# 修改一下 metadata 的行名
wellInfo[1:4,]
metadata = left_join(wellInfo, metadata, by = "plate")
metadata[1:4,]
wellInfo[1:4,]
rownames(metadata) = rownames(wellInfo)
# 判斷一下是否一緻
identical(colnames(raw.data),rownames(metadata))
           

> metadata[1:5,1:5]
well   plate sample_name patient_id sample_type
A10_B000561  A10 B000561       AZ_04  AZ008_NAT         NAT
A10_B000563  A10 B000563       AZ_04  AZ008_NAT         NAT
A10_B000568  A10 B000568       AZ_05      AZ008       tumor
A10_B001544  A10 B001544       AZ_01      AZ003       tumor
A10_B001546  A10 B001546       AZ_01      AZ003       tumor
           

接下來的步驟和上面完全一樣

# ERCC 是外源 RNA 主要的作用就是用來質控的
erccs = grep('^ERCC-', x= rownames(x=raw.data),value = T) # value = T 擷取名字
# 計算 ercc 比例
percent.ercc = Matrix::colSums(raw.data[erccs, ])/Matrix::colSums(raw.data)
fivenum(percent.ercc) # 檢視前 5 個結果
ercc.index = grep(pattern = "^ERCC-", x = rownames(x = raw.data), value = FALSE) # 擷取 index
# 删除 Count 矩陣中的 ERCC
raw.data = raw.data[-ercc.index,]
dim(raw.data) 
# [1] 26485  3507
# 建構對象
# 建構對象
# min.cells 某一個基因至少在多少個基因中表達
# min.features 某個細胞至少表達多少個基因
min.cells = 0 # 3 
min.features = 0 # 150 # 
sce = CreateSeuratObject(counts = raw.data,metadata = metadata,min.cells =min.cells,min.features =min.features)
# 添加metadata 和 外源基因資訊
sce = AddMetaData(object = sce,metadata = metadata)  
sce = AddMetaData(object = sce, percent.ercc, col.name = "percent.ercc") 
# 2.資料清洗
# 資料清洗原則 
# 1. 不符合品質的
# 2. 過濾線粒體 DNA
# 3. 過濾外原DNA
# 在過濾前先計算
table(grepl("^MT-",rownames(sce)))
table(grepl("^RP[SL][[:digit:]]",rownames(sce)))
sce[["percent.mt"]]  = PercentageFeatureSet(sce, pattern = "^MT-")
sce[["percent.rp"]]  = PercentageFeatureSet(sce, pattern = "^RP[SL][[:digit:]]")
dim(sce)
summary(sce[["nCount_RNA"]])
summary(sce[["nFeature_RNA"]])
summary(sce[["percent.mt"]]  )
summary(sce[["percent.ercc"]] )
summary(sce[["percent.rp"]]  )

# 過濾篩選
# 一般根據實際情況進行篩選，這次的篩選主要還是根據文章裡面的流程
myData2 = subset(x=sce, subset = nCount_RNA > 50000 & nFeature_RNA > 500)
myData2
# 26485 features across 2070 samples within 1 assay
save(myData2,file='myData2.Rdata')
           

第一部分結束，準備進入第二部分

第二部分 seurat 标準流程

小提琴圖觀察資料分布情況

rm(list = ls())
# 小提琴圖可視化
load("myData1.Rdata")
load("myData2.Rdata")

allData = merge(x = myData1, y = myData2)
table([email protected]$sample_name) %>% length()
VlnPlot(allData, features = 'nFeature_RNA', group.by = 'sample_name')
           

接下來就是走流程

先觀察一下資料的分布情況

rm(list = ls())
# 小提琴圖可視化
load("allData.Rdata")
raw_sce = allData
rm(allData)
plot1 = FeatureScatter(raw_sce, feature1 = "nCount_RNA", feature2 = "percent.ercc")
plot2 = FeatureScatter(raw_sce, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))

           

[外鍊圖檔轉存失敗,源站可能有防盜鍊機制,建議将圖檔儲存下來直接上傳(img-bex3RwOY-1614138377387)(https://tva1.sinaimg.cn/large/0081Kckwly1gkir1ibb7cj31jm0lc43b.jpg)]

流程：

1. log : NormalizeData
2. 找特征 : FindVariableFeatures
3. 标準化 : ScaleData
4. pca : RunPCA
5. 建構圖 : FindNeighbors
6. 聚類 : FindClusters
7. tsne /umap : RunTSNE RunUMAP
8. 差異基因 : FindAllMarkers / FindMarkers

VlnPlot(raw_sce, features = c("percent.rp", "percent.ercc"), ncol = 2)
           

# 開始走流程
# 1.log
sce = raw_sce
rm(raw_sce)
sce = NormalizeData(object = sce,normalization.method =  "LogNormalize",  scale.factor = 1e6)
# 2. 特征提取
sce = FindVariableFeatures(object = sce,selection.method = "vst", nfeatures = 2000)
# 所有特征
VariableFeatures(sce)
top20=head(VariableFeatures(sce),20)# 提取差異最大的 top20 基因
plot1= VariableFeaturePlot(sce)
plot2=LabelPoints(plot = plot1, points = top20, repel = TRUE)
plot2
           

# 3.标準化
sce = ScaleData(object = sce)
# 4. PCA
sce = RunPCA(object = sce, do.print = FALSE)
# 可視化
ElbowPlot(sce)
VizDimLoadings(sce, dims = 1:5, reduction = "pca", nfeatures = 10)
DimHeatmap(sce, dims = 1:10, cells = 100, balanced = TRUE)
           

[外鍊圖檔轉存失敗,源站可能有防盜鍊機制,建議将圖檔儲存下來直接上傳(img-MW8tmYr0-1614138377390)(https://tva1.sinaimg.cn/large/0081Kckwly1gkir227e5yj30u00z9dr5.jpg)]

# 5.建構圖
#細胞聚類
#首先我們在PCA空間裡根據歐氏距離建構一個KNN圖，并在任意兩個細胞間要确定它們的邊緣權重，這個過程#用FindNeighbors函數，這裡的input就是前面我們定義的dataset的維數。
# 選 20 個主成分
sce= FindNeighbors(sce, dims = 1:20)
           

# 6. 聚類
#再用FindClusters函數，該函數有一個“分辨率”的參數，該參數設定下遊聚類的“粒度”，值
#越高，得到的聚類數越多。這個參數設定在0.4-1.2之間，
#對于3千個左右的單細胞資料通常會得到
#比較好的結果。對于較大的資料集，最佳分辨率通常會增加。
sce = FindClusters(sce, resolution = 0.5) 
           

# 小劇場
# 不同的 resolution 分群的結果不一樣
# 下面介紹兩種可視化方法
# 首先先用兩種不同的resolution 進行計算
sce = FindClusters(sce, resolution = 0.2)
sce = FindClusters(sce, resolution = 0.8)

library(gplots)
tab.1=table([email protected]$RNA_snn_res.0.2,[email protected]$RNA_snn_res.0.8) 
# 這個是第一幅圖
balloonplot(tab.1)

# Check clustering stability at given resolution  
# Set different resolutions 
res.used = c(0.2,0.8)
res.used
# Loop over and perform clustering of different resolutions 
for(i in res.used){
  sce = FindClusters(object = sce, verbose = T, resolution = res.used)
}
# Make plot 
library(clustree)
clustree(sce) +
  theme(legend.position = "bottom") + 
  scale_color_brewer(palette = "Set1") +
  scale_edge_color_continuous(low = "grey80", high = "red")
# 把資料删除
[email protected] = [email protected][,-grep("snn_res.",ids,value = T)]
# 然後用文章給的 0.5
sce = FindClusters(sce, resolution = 0.5)
           

這兩張圖顯示了不同的分辨率情況下的分類結果(這兩張圖要是看不懂就真的極其的過分)

[外鍊圖檔轉存失敗,源站可能有防盜鍊機制,建議将圖檔儲存下來直接上傳(img-07YIOHvF-1614138377392)(https://tva1.sinaimg.cn/large/0081Kckwly1gkir387hzrj316a0u07gs.jpg)]

resolution 等于0.2的時候分了 19 類，等于 0.8 的時候分成了 31 類

# 7.tsne
#Seurat提供了幾種非線性的降維技術，如tSNE和UMAP。這些算法的目标是在低維空間中将相似的#細胞放在一起。上面所确定的基于圖的叢集中的單元應該在這些降維圖上共同定位。作為
#UMAP和tSNE的input，我們建議使用與聚類分析的輸入相同的PCs作為輸入。
sce=RunTSNE(sce,dims.use = 1:20)  ##tsne降維
sce=RunUMAP(sce,dims = 1:20)  ##umap降維
##可視化
p1=DimPlot(sce, reduction = "tsne")
p1
p2=DimPlot(sce, reduction = "umap")
p2
CombinePlots(plots =list(p1, p2))

           

[外鍊圖檔轉存失敗,源站可能有防盜鍊機制,建議将圖檔儲存下來直接上傳(img-IHPqaumn-1614138377393)(https://tva1.sinaimg.cn/large/0081Kckwly1gkir3g66vjj31gi0u01kx.jpg)]

# 用 ggplot 可視化
phe=data.frame(cell=rownames([email protected]),
               cluster [email protected]$seurat_clusters)
head(phe)
table(phe$cluster)
tsne_pos=Embeddings(sce,'tsne') 
head(phe)
table(phe$cluster)
head(tsne_pos) 
dat=cbind(tsne_pos,phe)
library(ggplot2)
p=ggplot(dat,aes(x=tSNE_1,y=tSNE_2,color=cluster))+geom_point(size=0.95)
p=p+stat_ellipse(data=dat,aes(x=tSNE_1,y=tSNE_2,fill=cluster,color=cluster),
                 geom = "polygon",alpha=0.2,level=0.9,type="t",linetype = 2,show.legend = F)+coord_fixed()
print(p) 
theme= theme(panel.grid =element_blank()) +   ## 删去網格
  theme(panel.border = element_blank(),panel.background = element_blank()) +   ## 删去外層邊框
  theme(axis.line = element_line(size=1, colour = "black")) 
p=p+theme+guides(colour = guide_legend(override.aes = list(size=5)))
print(p)
# umpa 也是一樣的就不畫了
           

# 8.差異基因
# 尋找某幾個簇之間的差異基因
#cluster5.markers = FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
# 找全部的差異基因（要非常久）
sce.markers <- FindAllMarkers(object = sce, only.pos = TRUE, min.pct = 0.25, 
                              thresh.use = 0.25)

DT::datatable(sce.markers)
write.csv(sce.markers,file=paste0(pro,'_sce.markers.csv'))
library(dplyr) 
top10 <- sce.markers %>% group_by(cluster) %>% top_n(10, avg_logFC)
# 這個圖要畫太久，是以隻選 15 個
DoHeatmap(sce,top10$gene[1:15],size=3)
ggsave(filename=paste0(pro,'_sce.markers_heatmap.pdf'))
sce.first=sce
save(sce.first,sce.markers,file = 'first_sce.Rdata')
rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
load(file = 'first_sce.Rdata')
sce = sce.first
rm(sce.first)

# 這個圖要畫非常就，所一就畫一幅基于可以了
# 這個圖是看合個組分之間的差異
for( i in 1 ){
  markers_df <- FindMarkers(object = sce, ident.1 = i, min.pct = 0.25)
  print(x = head(markers_df))
  markers_genes =  rownames(head(x = markers_df, n = 5))
  VlnPlot(object = sce, features =markers_genes,log =T )
  FeaturePlot(object = sce, features=markers_genes )
}

           

第二部分結束，準備進入第三部分

第三部分免疫細胞識别

免疫細胞的識别，一般情況下分為兩種

第一種就是你自己一直知道某些基因和自己要找的特定的免疫細胞相關

這裡你可以參考文獻，也可以自己定義，非常的自由

這種方法不止針對免疫細胞，而且針對所有你自己想 diy 識别的細胞，隻要能夠找到靶dna

例如以下這些代表基因：

epithelial/cancer (EpCAM+,EPCAM),

immune (CD45+,PTPRC),

stromal (CD10+,MME,fibo or CD31+,PECAM1,endo)

rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(ggplot2)
load(file = 'first_sce.Rdata')
sce=sce.first
rm(sce.first)

genes_to_check = c("PTPRC","EPCAM",'PECAM1','MME',"CD3G","CD3E", "CD79A")
p = DotPlot(sce, features = genes_to_check,
             assay='RNA' )  

p = p$data
p = p[which((p$pct.exp>50) & p$avg.exp.scaled>-0.5),] 

imm = p[which(p$features.plot == 'PTPRC'),]$id %>% as.character()
epi = p[which(p$features.plot == 'EPCAM'),]$id %>% as.character()
stromal=setdiff(as.character(0:27),c(imm,epi));

           

例如按照文章的定義

genes_to_check = c("PTPRC","EPCAM","CD3G","CD3E", "CD79A", "BLNK","MS4A1", "CD68", "CSF1R", 
                   "MARCO", "CD207", "PMEL", "ALB", "C1QB", "CLDN5", "FCGR3B", "COL1A1")
# All on Dotplot 
p = DotPlot(sce, features = genes_to_check) + coord_flip()
p
#然後根據自定義規則區分細胞
           

[外鍊圖檔轉存失敗,源站可能有防盜鍊機制,建議将圖檔儲存下來直接上傳(img-k0GADJ7e-1614138377395)(https://tva1.sinaimg.cn/large/0081Kckwly1gkir4ahlx2j31hf0u0tf5.jpg)]

[email protected]$immune_annotation = ifelse([email protected]$seurat_clusters  %in% imm ,'immune',
                                         ifelse([email protected]$seurat_clusters  %in% epi ,'epi','stromal') )
# MAke a table 
table([email protected]$immune_annotation)

p = TSNEPlot(object = sce, group.by = 'immune_annotation')
p 

genes_to_check = c("PTPRC","EPCAM","CD3G","CD3E", "CD79A", "BLNK","MS4A1", "CD68", "CSF1R", 
                   "MARCO", "CD207", "PMEL", "ALB", "C1QB", "CLDN5", "FCGR3B", "COL1A1")
# All on Dotplot 
p = DotPlot(sce, features = genes_to_check,group.by = 'immune_annotation') + coord_flip()
p
# Generate immune and non-immune cell lists
[email protected]
table(phe$immune_annotation)
save(phe,file = 'phe-of-first-anno.Rdata')
           

[外鍊圖檔轉存失敗,源站可能有防盜鍊機制,建議将圖檔儲存下來直接上傳(img-5KHP0V24-1614138377396)(https://tva1.sinaimg.cn/large/0081Kckwly1gkir4m54zej30u00yb7wh.jpg)]

其實，細胞注釋到這裡已經結束了，其實不難發現，細胞注釋這個流程，沒有什麼難的，隻需要做兩件事情

1. 定義基因，根據基因集合畫點圖 DotPlot + coord_flip()
2. 根據點圖自己定義什麼是什麼細胞

那麼到這裡的話基本上面免疫細胞大群體已經區分開了，下一步就是直接對大群體内部細分

rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(ggplot2)
load(file = 'first_sce.Rdata')  
load(file = 'phe-of-first-anno.Rdata')
sce=sce.first
table(phe$immune_annotation)
[email protected]=phe 
# 既然是對免疫群體内部分類，則隻需要提取免疫細胞就行
cells.use = row.names([email protected])[which(phe$immune_annotation=='immune')]
length(cells.use)
sce =subset(sce, cells=cells.use)  
sce 
# 走标準流程，一直走到 TSNE
#1.Normalize data
require(Seurat)
sce = NormalizeData(object = sce, normalization.method =  "LogNormalize",scale.factor = 1e6)

#2.Find variable genes
sce = FindVariableFeatures(object = sce, selection.method = "vst", nfeatures = 2000)

#3.Scale data
sce = ScaleData(object = sce)

#4.Perform PCA
sce = RunPCA(object = sce, pc.genes = VariableFeatures(sce))
ElbowPlot(sce)
VizDimLoadings(sce, dims = 1:5, reduction = "pca", nfeatures = 10)
DimHeatmap(sce, dims = 1:10, cells = 100, balanced = TRUE)
pca.obj = sce@reductions$pca
pc.coords = [email protected]
df1 = [email protected][,c("nFeature_RNA","nCount_RNA","percent.rp")]
df2 = pc.coords[,c(1:10)]
cordf12 = cor(df1,df2)
# Make a correlation plot
library(corrplot)
# 補一張相關性圖
corrplot(cordf12, method = "number", main="Correlation of PCs and metadata")

#5. Construct Neighbor graph  
sce = FindNeighbors(object = sce, dims = 1:15)

#6. 分辨率
tiss_subset = FindClusters(object = tiss_subset, verbose = T, resolution = 1)

#7. Run and project TSNEs
sce = RunTSNE(sce, dims = 1:15)
# 可視化
DimPlot(sce, reduction = "tsne", label = T)
DimPlot(sce,reduction = "tsne",label=T, group.by = "patient_id")
TSNEPlot(object = sce, do.label = F, group.by = "driver_gene")
TSNEPlot(object = sce, do.label = F, group.by = "patient_id")
tab=table([email protected]$seurat_clusters,[email protected]$patient_id)
balloonplot(tab)
           

根據參考資料庫自動化注釋

sce_for_SingleR = GetAssayData(sce, slot="data")
sce_for_SingleR
library(SingleR)
# 根據資料庫進行自動注釋
hpca.se = HumanPrimaryCellAtlasData() # 這個要搞非常久
hpca.se
[email protected]$seurat_clusters
clusters[1:5]
# test 稀疏矩陣
# ref 參考資料庫
table( hpca.se$label.main)

pred.hesc = SingleR(test = sce_for_SingleR, ref = hpca.se, labels = hpca.se$label.main,
                     method = "cluster", clusters = clusters, 
                     assay.type.test = "logcounts", assay.type.ref = "logcounts")
table(pred.hesc$labels)
celltype = data.frame(ClusterID=rownames(pred.hesc), celltype=pred.hesc$labels, stringsAsFactors = F) 
[email protected]$singleR=celltype[match(clusters,celltype$ClusterID),'celltype']
DimPlot(sce, reduction = "tsne", group.by = "singleR")
[email protected]
table(phe$singleR)
save(phe,file = 'phe-of-subtypes-Immune-by-singleR.Rdata')
           

[外鍊圖檔轉存失敗,源站可能有防盜鍊機制,建議将圖檔儲存下來直接上傳(img-eV5QDrSR-1614138377396)(https://tva1.sinaimg.cn/large/0081Kckwly1gkir4we5j5j31ir0u01kx.jpg)]

小結

那麼這裡再一次的簡單總結一下細胞注釋的流程（并不隻是免疫細胞，任何的都可以）:

首先我們需要找出和要注釋的細胞的相關的基因

graph TB
找出注釋細胞相關基因-->繪制點圖
繪制點圖-->根據自定義門檻值确定細胞類型
根據自定義門檻值确定細胞類型-->提取感興趣的類
提取感興趣的類-->走标準流程1-7
走标準流程1-7-->資料庫自動化注釋

           

不同的注釋庫的功能不太一樣（這些都是大佬整理好的，講實話我很好奇這些資料庫都是哪兒找的）

?MonacoImmuneData

?BlueprintEncodeData

?DatabaseImmuneCellExpressionData

?NovershternHematopoieticData

?MonacoImmuneData

?ImmGenData

?MouseRNAseqData

?HumanPrimaryCellAtlasData

其他的細胞類型注釋流程幾乎完全一樣，如果後續需要個性化分析，隻需要提取出自己識别好的類，然後在這些類中繼續走标準流程 1-7，之後還是使用點圖進行驗證，單純的隻是多了個 group.by 參數而已

第三部分結束，準備進入第四部分

第四部分癌細胞識别

這個癌細胞識别略微小複雜，通常的流程是在已知的細胞類别上繼續進行分類。

這種分類方式大部分也都是屬于人工自己篩選的

例如在上皮細胞中識别癌細胞和非癌細胞

rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(ggplot2)
load(file = 'first_sce.Rdata')  
load(file = 'phe-of-first-anno.Rdata')
sce=sce.first
rm(sce.first)
table(phe$immune_annotation)
[email protected]=phe
# 那麼第一步就是提取上皮細胞中的癌細胞
cells.use = row.names([email protected])[which(phe$immune_annotation=='epi')]
length(cells.use)
sce = subset(sce, cells=cells.use)  
sce
# 走标準流程，一直走到 TSNE
#1.Normalize data
require(Seurat)
sce = NormalizeData(object = sce, normalization.method =  "LogNormalize",scale.factor = 1e6)

#2.Find variable genes
sce = FindVariableFeatures(object = sce, selection.method = "vst", nfeatures = 2000)

#3.Scale data
sce = ScaleData(object = sce)

#4.Perform PCA
sce = RunPCA(object = sce, pc.genes = VariableFeatures(sce))
ElbowPlot(sce)
VizDimLoadings(sce, dims = 1:5, reduction = "pca", nfeatures = 10)
DimHeatmap(sce, dims = 1:10, cells = 100, balanced = TRUE)
pca.obj = [email protected]$pca
pc.coords = [email protected]
df1 = [email protected][,c("nFeature_RNA","nCount_RNA","percent.rp")]
df2 = pc.coords[,c(1:10)]
cordf12 = cor(df1,df2)
# Make a correlation plot
library(corrplot)
# 補一張相關性圖
corrplot(cordf12, method = "number", main="Correlation of PCs and metadata")

#5. Construct Neighbor graph  
sce = FindNeighbors(object = sce, dims = 1:15)

#6. 分辨率
tiss_subset = FindClusters(object = tiss_subset, verbose = T, resolution = 1)

#7. Run and project TSNEs
sce = RunTSNE(sce, dims = 1:15)
# 可視化
DimPlot(sce, reduction = "tsne", label = T)
DimPlot(sce,reduction = "tsne",label=T, group.by = "patient_id")
TSNEPlot(object = sce, do.label = F, group.by = "driver_gene")
TSNEPlot(object = sce, do.label = F, group.by = "patient_id")
tab=table([email protected]$seurat_clusters,[email protected]$patient_id)
balloonplot(tab)
           

這裡的話就存在一個前提假設，如果細胞會跨患者聚類則這些細胞不是癌細胞

[外鍊圖檔轉存失敗,源站可能有防盜鍊機制,建議将圖檔儲存下來直接上傳(img-JUUMdTTy-1614138377397)(https://tva1.sinaimg.cn/large/0081Kckwly1gkir56lncoj31hm0u01kx.jpg)]

tab = tab %>% as.matrix()
ids = apply(tab,1,function(x){ x/sum(x) })
nonCancer = which(apply(ids,2,max)<0.85) %>% colnames(ids)[.]

[email protected]$cancer <-ifelse([email protected]$seurat_clusters %in% nonCancer ,'non-cancer','cancer')
# MAke a table 
table([email protected]$cancer)
[email protected]

# cancer marker CEACAM6, cell-cycle-related gene CCND2, and apoptosis-related gene BAX 
genes_to_check = c("CEACAM6","CCND2","BAX")
# All on Dotplot 
p <- DotPlot(sce, features = genes_to_check,group.by = 'cancer') + coord_flip()
p 

save(phe,file = 'phe-of-cancer-or-not.Rdata')
           

其實到這裡癌細胞注釋和上一步的免疫細胞注釋幾乎一樣

是以下面介紹 inferCNV 的算法

在使用 infercnv其實隻有 2 行代碼而已

重要的是引導檔案的建構

1. 第一個引導檔案是 count 矩陣
2. 第二個引導檔案是 label 資訊 (第一列是細胞 ID；第二列是細對應的細胞類型，這個資料是由前面的分析結果得到的)
3. 第三個引導檔案是樣本染色體或基因資訊（這個擷取方法有很多種，1.抄大佬代碼（推薦），2.自己下載下傳）

rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(ggplot2)
load(file = 'first_sce.Rdata')  
load(file = 'phe-of-first-anno.Rdata')
sce=sce.first
table(phe$immune_annotation)
[email protected]=phe 
table(phe$immune_annotation,phe$seurat_clusters) 
# BiocManager::install("infercnv")
library(infercnv)
##### inferCNV引導檔案的建構
# 上皮細胞 count 矩陣
epi.cells  = row.names([email protected])[which(phe$immune_annotation=='epi')]
epiMat=as.data.frame(GetAssayData(subset(sce,
                                         cells=epi.cells), slot='counts'))
# 間質細胞
cells.use = row.names([email protected])[which(phe$immune_annotation=='stromal')]
sce =subset(sce, cells=cells.use)  
load(file = 'phe-of-subtypes-stromal.Rdata')
fib.cells  =  row.names([email protected])[phe$singleR=='Fibroblasts']
endo.cells  =  row.names([email protected])[phe$singleR=='Endothelial_cells']

# 成纖細胞 count 矩陣
fibMat=as.data.frame(GetAssayData(subset(sce,
                                         cells=fib.cells),  slot='counts'))
# 内皮細胞 count 矩陣
endoMat=as.data.frame(GetAssayData(subset(sce,
                                          cells=endo.cells),  slot='counts'))
#合并
dat=cbind(epiMat,fibMat,endoMat)
ids1 = data.frame(cellId = colnames(epiMat))
ids2 = data.frame(cellId = colnames(fibMat))
ids3 = data.frame(cellId = colnames(endoMat))
ids1$cellType = 'epi'
ids2$cellType = 'Fibroblasts'
ids3$cellType = 'Endothelial_cells'
data.frame(ids1,ids2,ids3)
groupFiles = rbind(ids1,ids2) %>% rbind(.,ids3) %>% as.data.frame
rm(epiMat,fibMat,endoMat,ids1,ids2,ids3)

library(AnnoProbe)
geneInfor=annoGene(rownames(dat),"SYMBOL",'human')
geneInfor=geneInfor[with(geneInfor, order(chr, start)),c(1,4:6)]
geneInfor=geneInfor[!duplicated(geneInfor[,1]),]
## 這裡可以去除性染色體
# 也可以把染色體排序方式改變
dat=dat[rownames(dat) %in% geneInfor[,1],]
dat=dat[match( geneInfor[,1], rownames(dat) ),] 
dim(dat)
expFile='expFile.txt'
write.table(dat,file = expFile,sep = '\t',quote = F)
groupFiles='groupFiles.txt'
head(groupinfo)
write.table(groupinfo,file = groupFiles,sep = '\t',quote = F,col.names = F,row.names = F)
head(geneInfor)
geneFile='geneFile.txt'
write.table(geneInfor,file = geneFile,sep = '\t',quote = F,col.names = F,row.names = F)

table(groupinfo[,2])
           

以上就是建構三個引導檔案的代碼

最終引導檔案分别張這樣：

> dat[1:5,1:5] #count 矩陣
          A10_1001000407 A10_1001000408 A10_1001000412 A10_B000863 A10_B001470
DDX11L1                0              0              0           0           0
WASH7P                 0              0              0           0           0
MIR6859-1              0              0              0           0           0
MIR1302-2              0              0              0           0           0
FAM138A                0              0              0           0           0

> groupinfo[1:5,]# 分組矩陣
          cellId cellType
1 A10_1001000407      epi
2 A10_1001000408      epi
3 A10_1001000412      epi
4    A10_B000863      epi
5    A10_B001470      epi

> geneInfor[1:5,]#染色體或基因資訊
     SYMBOL  chr start   end
1   DDX11L1 chr1 11869 14409
2    WASH7P chr1 14404 29570
3 MIR6859-1 chr1 17369 17436
5 MIR1302-2 chr1 30366 30503
6   FAM138A chr1 34554 36081
           

然後就直接跑流程就可以

# 這個注意一下這個流程跑起來非常慢，不知道是我的電腦的問題還是怎麼回事
# 大概是我中午睡了一覺起來還沒好
rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(ggplot2)
library(infercnv)
expFile='expFile.txt' 
groupFiles='groupFiles.txt'  
geneFile='geneFile.txt'
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=expFile,
                                    annotations_file=groupFiles,
                                    delim="\t",
                                    gene_order_file= geneFile,
                                    ref_group_names=c('Fibroblasts','Endothelial_cells'))  ## 這個是 normalcell 的 grouplabel


##  文獻的代碼：
infercnv_obj2 = infercnv::run(infercnv_obj,
                              cutoff=1,  
                              out_dir=  'plot_out/inferCNV_output2' , 
                              cluster_by_groups=F,   # cluster
                              hclust_method="ward.D2", plot_steps=F) 

### 可以添加 denoise=TRUE以及 HMM=TRUE
           

最後就是 inferCNV 的結果解讀了

link1

link2

全文總結

其實，單細胞從頭到尾的流程就那麼幾個代碼，感覺有 70%以上的都是重複的

1. 資料預處理，確定 count 資料能夠和 metadata 一一對應
2. 計算外源RNA的百分比，添加到 metadata；删除 count 矩陣中的外源 RNA 資訊
3. 建構 seurat 對象，計算 MT，RP 百分比資訊
4. 根據相應的資訊自主選擇過濾，一般都是過濾基因，細胞，MT 百分比，RP 百分比
5. 第一輪基礎流程，然後儲存資料
6. 後續的所謂的個性化分析都是基于特定的細胞群體
7. 細胞群體的識别方法按這個流程的話隻有一種，就是找到特定的基因，畫點圖然後提取
8. 反複的根據提取的子細胞群落用第 5 步存下來的資料一直跑基礎流程，一直畫點圖，TSNE
9. 如果做 inferCNV 隻要做好引導檔案即可，除了修改

   ref_group_names

   以外不懂的情況下一個字别改

如果您覺得我的文章對您有幫助，請點贊+關注，可以的話打個賞獎勵一杯星巴克(～￣(OO)￣)ブ

Best Regards,  
Yuan.SH;
School of Basic Medical Sciences,  
Fujian Medical University,  
Fuzhou, Fujian, China.  
please contact with me via the following ways:  
(a) e-mail :[email protected]