1 亞群合并

當(dāng)有幾類亞群同屬于某類細(xì)胞時(shí)，比如CD4+ T細(xì)胞和CD8+ T細(xì)胞均屬于T細(xì)胞，想要將他們合并在一起時(shí)，可以使用此代碼。

#接seurat標(biāo)準(zhǔn)流程代碼
#命名pbmc為sce，表明其是seurat對(duì)象
sce=pbmc

Idents(sce)#細(xì)胞身份，屬于哪個(gè)亞群
levels(sce)#獲取sce的亞群名
> [1] "Naive CD4 T"  "CD14+ Mono"   "Memory > CD4 T"
> [4] "B"            "CD8 T"        "FCGR3A+ > Mono"
> [7] "NK"           "DC"           "Platelet"
#可以看到Naive CD4 T，Memory CD4T及CD8 T都屬于T細(xì)胞

head(sce@meta.data) #可以看到meta.data中多了seurat_clusters列，由0-9構(gòu)成。
# 方法1 

#根據(jù)levels(sce)可將要合并的亞群寫作如下
new.cluster.ids <- c("T", "Mono", "T", 
                     "B", "T", "Mono",
                     "T", "DC", "Platelet")

names(new.cluster.ids) <- levels(sce)
#更改前
p1<-DimPlot(pbmc, reduction = 'umap', 
            label = TRUE, pt.size = 0.5) + NoLegend()

#重命名亞群
sce <- RenameIdents(sce, new.cluster.ids)
#更改后
p2<-DimPlot(sce, reduction = 'umap', 
        label = TRUE, pt.size = 0.5) + NoLegend()


#方法2
根據(jù)亞群對(duì)應(yīng)關(guān)系，創(chuàng)建一個(gè)向量
cluster2celltype <- c("0"="T",
                  "1"="Mono", 
                  "2"="T", 
                  "3"= "B", 
                  "4"= "T", 
                  "5"= "Mono",
                  "6"= "T", 
                  "7"= "DC", 
                  "8"= "Platelet") 
#unname去掉對(duì)象的名（行名列名）
sce[['cell_type']] = unname(cluster2celltype[sce@meta.data$seurat_clusters])
#繪制散點(diǎn)圖，reduction：使用哪種降維。如果沒(méi)有指定，首先搜索 umap，然后是 tsne，然后是 pca；label:是否標(biāo)記分群
DimPlot(sce, reduction = 'umap', group.by = 'cell_type',
        label = TRUE, pt.size = 0.5) + NoLegend()#刪除圖例

p3<-DimPlot(sce, reduction = 'umap', group.by = 'cell_type',
        label = TRUE, pt.size = 0.5) + NoLegend()#刪除圖例
p1+p2+p3#此函數(shù)需要library(dplyr)

推薦使用第二種方法，這樣不更改原始亞群分群，只是在metadata中增加了一列

2 提取子集

當(dāng)我們想把表達(dá)感興趣基因的細(xì)胞提取出來(lái)單獨(dú)分析時(shí)，可使用此函數(shù)。

p1<-DimPlot(pbmc, reduction = 'umap', group.by ="seurat_clusters",label = TRUE, pt.size = 0.5) + NoLegend()
p2<-FeaturePlot(pbmc, features = c("CD4"))
p3<-DotPlot(sce, group.by = 'seurat_clusters',
           features = unique(genes_to_check)) + RotatedAxis()

p1+p2+p3

比如我們想把某幾個(gè)亞群提取出來(lái)

sce=pbmc
#R 中，提取子集, 需要知道邏輯值，坐標(biāo)及名字
# seurat 取0，2分群的子集
cd4_sce1 = sce[,sce@meta.data$seurat_clusters %in% c(0,2)]
#取Naive CD4 T ,  Memory CD4 T 
cd4_sce2 = sce[, Idents(sce) %in% c( "Naive CD4 T" ,  "Memory CD4 T" )]

取子集后就是新的項(xiàng)目，需要重新標(biāo)準(zhǔn)化，重新用子集走一遍seurat標(biāo)準(zhǔn)流程

3 如何分群是合理的

我們知道FindClusters函數(shù)中不同的resolution參數(shù)會(huì)帶來(lái)不同的結(jié)果，而且即使某個(gè)亞群內(nèi)部的細(xì)胞也會(huì)有一定的異質(zhì)性，那么到底分群多少是合適的呢？

#先執(zhí)行不同resolution 下的分群
library(Seurat)
library(clustree)
sce <- FindClusters(
  object = sce,
  resolution = c(seq(from=.1,to=1.6,by=.2))
) 
#clustree創(chuàng)建一個(gè)聚類樹(shù)圖，依賴于clustree包，顯示不同分辨率的聚類之間的關(guān)系。prefix指示包含聚類信息的列的字符串
colnames(sce@meta.data)
pz<-clustree(sce@meta.data, prefix = "RNA_snn_res.")

p2=DimPlot(sce, reduction = 'umap', group.by = 'RNA_snn_res.0.5',
           label = TRUE, pt.size = 0.5) + NoLegend()
p3=DimPlot(sce, reduction = 'umap',group.by = 'RNA_snn_res.1.5',
           label = TRUE, pt.size = 0.5) + NoLegend()
library(patchwork)
p1+p2

根據(jù)clustree可確定合適的分群數(shù)，若resolution設(shè)置過(guò)大，會(huì)導(dǎo)致分群過(guò)度，把原來(lái)分群的中應(yīng)有的異質(zhì)性也提煉出來(lái)單獨(dú)作為一群。無(wú)論如何分群，都應(yīng)該結(jié)合FindMarkers函數(shù)，確認(rèn)亞群標(biāo)志基因，結(jié)合生物學(xué)背景來(lái)解釋分群結(jié)果。

4 細(xì)胞亞群等比例抽樣

當(dāng)數(shù)據(jù)特別大時(shí)，可以采用等比例抽樣的方式，降低運(yùn)算時(shí)間。

sce=pbmc

p1<-DoHeatmap(sce , 
          features = features, 
          size = 3
          )
table(Idents(sce))
#使用subset函數(shù)，每個(gè)亞群抽取15個(gè)細(xì)胞，做熱圖
p2<-DoHeatmap(subset(sce, downsample = 15), 
          features = features, size = 3)+
  labs(title = paste("downsample = 15"))
p1+p2

# 每個(gè)細(xì)胞亞群抽10 
allCells=names(Idents(sce))
allType = levels(Idents(sce))

choose_Cells = unlist(lapply(allType, function(x){
  cgCells = allCells[Idents(sce)== x ]
  cg=sample(cgCells,10) #sample 從指定的元素中獲取指定大小的樣本。
  cg
}))
#將抽取出的細(xì)胞組成一個(gè)新的sce對(duì)象
cg_sce = sce[, allCells %in% choose_Cells]
> An object of class Seurat 
> 13714 features across 90 samples within 1 > assay 
> Active assay: RNA (13714 features, 2000 > > variable features)
>  2 dimensional reductions calculated: pca, umap

as.data.frame(table(Idents(cg_sce)))

5 如何查看基因的原始表達(dá)量

# 如何看原始表達(dá)量 slot:要使用的數(shù)據(jù)槽，從“raw.data”、“data”或“scale.data”中選擇;size 顏色條上方文字大小
#不加slot默認(rèn)是從之前2000個(gè)FindVariableFeatures基因作圖
DoHeatmap(subset(sce ), 
          features = features, 
          size = 3, slot='data'
)

6 多個(gè)亞群的多個(gè)標(biāo)記基因可視化

#首先需要將各個(gè)亞群的標(biāo)記基因找出來(lái)
sce.markers <- FindAllMarkers(object = sce, 
               only.pos = TRUE,min.pct = 0.25, thresh.use = 0.25)

#此時(shí)可以創(chuàng)建一個(gè)HTML小部件來(lái)顯示矩形數(shù)據(jù)
DT::datatable(sce.markers) 

可以以交互形式的HTML格式查看各亞群的標(biāo)記基因

#挑選各個(gè)亞群差異基因的前5名
top5 <- sce.markers %>% group_by(cluster) %>% top_n(5,avg_log2FC)

head(top5)
#檢測(cè)并去除重復(fù)值 ！表示“與、或、非”中的“非”的意思
top5=top5[!duplicated(top5$gene),] 
select_genes_all=split(top5$gene,top5$cluster)#將表拆分成列表
select_genes_all
DotPlot(object = sce, features=select_genes_all, 
        assay = "RNA") + theme(axis.text.x = element_text(angle= 45,  vjust = 0.5, hjust=0.5))

不同亞群多個(gè)差異基因的可視化

7 對(duì)任意細(xì)胞亞群做差異分析

#若沒(méi)有運(yùn)行過(guò)FindAllMarkers，運(yùn)行FindAllMarkers函數(shù)，查找各亞群標(biāo)記基因
sce.markers <- FindAllMarkers(object = sce, 
               only.pos = TRUE,min.pct = 0.25, thresh.use = 0.25)

#table使用交叉分類的因素建立一個(gè)列聯(lián)表,計(jì)數(shù)在每個(gè)組合的因素水平。
table(sce.markers$cluster)

Naive CD4 T   CD14+ Mono Memory CD4 T            B        CD8 T 
         164          391          178          148          144 
FCGR3A+ Mono           NK           DC     Platelet 
         608          369          633          242 

挑選亞群和對(duì)應(yīng)基因

# 挑選細(xì)胞，若挑選多個(gè)亞群細(xì)胞，可采用正則表達(dá)式選擇，比如挑選 Mono和T細(xì)胞：grepl(("Mono|T") , Idents(sce ))，這里只挑選Mono類亞群
kp=grepl(("Mono") , Idents(sce ))   
table(kp)
sce=sce[,kp]
sce
table( Idents(sce ))

# 挑選對(duì)應(yīng)的標(biāo)記基因
kp=grepl('Mono',sce.markers$cluster)#返回匹配與否的邏輯值
table(kp)
cg_sce.markers = sce.markers [ kp ,]
#這里認(rèn)為vg_log2FC>2的才是標(biāo)記基因
cg_sce.markers=cg_sce.markers[cg_sce.markers$avg_log2FC>2,]

查找兩群的差異基因

markers_df <- FindMarkers(object = sce, 
                 ident.1 = 'FCGR3A+ Mono',#ident 聚類名
                ident.2 = 'CD14+ Mono',
                 #logfc.threshold = 0,
                  min.pct = 0.25)
                  
head(markers_df)

篩選差異基因

#選取avg_log2FC絕對(duì)值大于1基因
cg_markers_df=markers_df[abs(markers_df$avg_log2FC) >1,]

dim(cg_markers_df)

DoHeatmap(subset(sce, downsample = 50),
          slot = 'counts',
          unique(rownames(cg_markers_df)),size=3) 
#取兩個(gè)基因集的交集（基因在某個(gè)群中即高表達(dá)，又有差異）
intersect(rownames(cg_markers_df) ,
          cg_sce.markers$gene)

8 人為劃分亞群挑選差異基因

當(dāng)我們對(duì)表達(dá)某類基因的細(xì)胞感興趣時(shí)，可以把表達(dá)這類基因的的細(xì)胞挑選出來(lái)，單獨(dú)劃分為一個(gè)群。

highCells= colnames(subset(x = sce, subset = FCGR3A > 1,slot = 'counts')) 

#a %in% table，a值是否包含于table中，為真輸出TURE，否者輸出FALSE
highORlow=ifelse(colnames(sce) %in% highCells,'high','low')

table(highORlow)
table(Idents(sce))
table(Idents(sce),highORlow)

#將highORlow添加到meta.data中
sce@meta.data$highORlow<-highORlow

#查找劃分出的亞群與其他亞群的差異基因
markers <- FindMarkers(sce, ident.1 = "high", group.by = 'highORlow' )
head(markers)

cg_markers=markers[abs(markers$avg_log2FC) >1,]

dim(cg_markers)

DoHeatmap(subset(sce, downsample = 15),
          rownames(cg_markers) ,size=3 ) 

9 查看任意文獻(xiàn)的單細(xì)胞亞群標(biāo)記基因

可以直接使用DotPlot函數(shù)以文獻(xiàn)中標(biāo)記基因作圖

#若想查看下述marker_genes在T細(xì)胞中的表達(dá)，首先需要將T細(xì)胞相關(guān)亞群提取出來(lái)

sce=pbmc
table(Idents(sce))
#提取單細(xì)胞相關(guān)亞群
t_sce = sce[, Idents(sce) %in% c("Naive CD4 T" ,  "Memory CD4 T" , "CD8 T","NK")]  

# 然后進(jìn)行標(biāo)準(zhǔn)的降維聚類分群,代碼不要變動(dòng)
sce=t_sce
sce <- NormalizeData(sce, normalization.method = "LogNormalize",
                     scale.factor = 1e4) 
sce <- FindVariableFeatures(sce, selection.method = 'vst',
                            nfeatures = 2000)
sce <- ScaleData(sce, vars.to.regress = "percent.mt")
sce <- RunPCA(sce, features = VariableFeatures(object = sce)) 
sce <- FindNeighbors(sce, dims = 1:10)
sce <- FindClusters(sce, resolution = 1 )
head(Idents(sce), 5)
table(sce$seurat_clusters) 
sce <- RunUMAP(sce, dims = 1:10)
DimPlot(sce, reduction = 'umap')

查看文獻(xiàn)中提到的基因

marker_genes=c("LEF1","TCF7","SELL","IL7R","CD40LG","ANXA1","FOS","JUN","FOXP3","SAT1","IL2RA","CTLA4","PDCD1","CXCL13","CD200","TNFRSF18","CCR7","NELL2","CD55","KLF2","TOB1","ZNF683","CCL5","GZMK","EOMES","ITM2C","CX3CR1","GNLY","GZMH","GZMB","LAG3","CCL4L2","FCGR3A","FGFBP2","TYROBP","AREG","XCL1","KLRC1","TRDV2","TRGV9","MTRNR2L8","KLRD1","TRDV1","KLRC3","CTSW","CD7","MKI67","STMN1","TUBA1B","HIST1H4C" )

DotPlot(sce, features = marker_genes,
             assay='RNA'  )  + coord_flip() +
  theme(axis.text.x = element_text(angle = 90))

10 多個(gè)單細(xì)胞亞群各自標(biāo)記基因聯(lián)合進(jìn)行GO及kegg數(shù)據(jù)庫(kù)注釋

#加載需要的包
library(Seurat)
library(gplots)
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)

#加載上述cg_sce.markers
pro='all'
load('sce.markers.all_10_celltype.Rdata')
table(sce.markers$cluster)

#基因ID轉(zhuǎn)換bitr(geneID, fromType, toType, OrgDb, drop = TRUE)
ids=bitr(sce.markers$gene,'SYMBOL','ENTREZID','org.Hs.eg.db')

#合并數(shù)據(jù)
sce.markers=merge(sce.markers,ids,by.x='gene',by.y='SYMBOL')

#拆分id和分群，將ENTREZID拆分成cluster定義的組 
gcSample=split(sce.markers$ENTREZID, sce.markers$cluster)

#比較亞群功能
xx <- compareCluster(gcSample, fun="enrichKEGG",organism="hsa", pvalueCutoff=0.05)
p=dotplot(xx) 
p+ theme(axis.text.x = element_text(angle = 45,  vjust = 0.5, hjust=0.5))
p

對(duì)相關(guān)基因進(jìn)行富集分析

#將上升下降的差異基因挑選出來(lái)
deg=FindMarkers(object = sce, ident.1 = 'FCGR3A+ Mono',ident.2 = 'CD14+ Mono', #logfc.threshold = 0,min.pct = 0.25)

gene_up=rownames(deg[deg$avg_log2FC>0,])
gene_down=rownames(deg[deg$avg_log2FC < 0,])

library(org.Hs.eg.db)
#把SYMBOL改為ENTREZID，這里會(huì)損失一部分無(wú)法匹配到的
gene_up=as.character(na.omit(AnnotationDbi::select(org.Hs.eg.db, keys = gene_up,columns = 'ENTREZID',keytype = 'SYMBOL')[,2]))

gene_down=as.character(na.omit(AnnotationDbi::select(org.Hs.eg.db,keys = gene_down,columns = 'ENTREZID',keytype = 'SYMBOL')[,2]))

library(clusterProfiler)

這里用到了生信技能樹(shù)-健明老師的一鍵分析代碼，感謝健明老師的無(wú)私分享

## KEGG pathway analysis
### 做KEGG數(shù)據(jù)集超幾何分布檢驗(yàn)分析，重點(diǎn)在結(jié)果的可視化及生物學(xué)意義的理解。
run_kegg <- function(gene_up,gene_down,pro='test'){
  library(ggplot2)
  gene_up=unique(gene_up)
  gene_down=unique(gene_down)
  gene_diff=unique(c(gene_up,gene_down))
  ###   over-representation test
  # 下面把3個(gè)基因集分開(kāi)做超幾何分布檢驗(yàn)
  # 首先是上調(diào)基因集。
  kk.up <- enrichKEGG(gene         = gene_up,
                      organism     = 'hsa',
                      #universe     = gene_all,
                      pvalueCutoff = 0.9,
                      qvalueCutoff =0.9)
  head(kk.up)[,1:6]
  kk=kk.up
  dotplot(kk)
  kk=DOSE::setReadable(kk, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
  write.csv(kk@result,paste0(pro,'_kk.up.csv'))
  
  # 首先是下調(diào)基因集。
  kk.down <- enrichKEGG(gene         =  gene_down,
                        organism     = 'hsa',
                        #universe     = gene_all,
                        pvalueCutoff = 0.9,
                        qvalueCutoff =0.9)
  head(kk.down)[,1:6]
  kk=kk.down
  dotplot(kk)
  kk=DOSE::setReadable(kk, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
  write.csv(kk@result,paste0(pro,'_kk.down.csv'))
  
  # 最后是上下調(diào)合并后的基因集。
  kk.diff <- enrichKEGG(gene         = gene_diff,
                        organism     = 'hsa',
                        pvalueCutoff = 0.05)
  head(kk.diff)[,1:6]
  kk=kk.diff
  dotplot(kk)
  kk=DOSE::setReadable(kk, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
  write.csv(kk@result,paste0(pro,'_kk.diff.csv'))
  
  
  kegg_diff_dt <- as.data.frame(kk.diff)
  kegg_down_dt <- as.data.frame(kk.down)
  kegg_up_dt <- as.data.frame(kk.up)
  down_kegg<-kegg_down_dt[kegg_down_dt$pvalue<0.01,];down_kegg$group=-1
  up_kegg<-kegg_up_dt[kegg_up_dt$pvalue<0.01,];up_kegg$group=1
  #畫圖設(shè)置, 這個(gè)圖很丑，大家可以自行修改。
  g_kegg=kegg_plot(up_kegg,down_kegg)
  print(g_kegg)
  
  ggsave(g_kegg,filename = paste0(pro,'_kegg_up_down.png') )
  
  if(F){
    ###  GSEA 
    ## GSEA算法跟上面的使用差異基因集做超幾何分布檢驗(yàn)不一樣。
    kk_gse <- gseKEGG(geneList     = geneList,
                      organism     = 'hsa',
                      nPerm        = 1000,
                      minGSSize    = 20,
                      pvalueCutoff = 0.9,
                      verbose      = FALSE)
    head(kk_gse)[,1:6]
    gseaplot(kk_gse, geneSetID = rownames(kk_gse[1,]))
    gseaplot(kk_gse, 'hsa04110',title = 'Cell cycle') 
    kk=DOSE::setReadable(kk_gse, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
    tmp=kk@result
    write.csv(kk@result,paste0(pro,'_kegg.gsea.csv'))
    
    
    # 這里找不到顯著下調(diào)的通路，可以選擇調(diào)整閾值，或者其它。
    down_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore < 0,];down_kegg$group=-1
    up_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore > 0,];up_kegg$group=1
    
    g_kegg=kegg_plot(up_kegg,down_kegg)
    print(g_kegg)
    ggsave(g_kegg,filename = paste0(pro,'_kegg_gsea.png'))
    # 
  }
  
}

### GO database analysis 
### 做GO數(shù)據(jù)集超幾何分布檢驗(yàn)分析，重點(diǎn)在結(jié)果的可視化及生物學(xué)意義的理解。
run_go <- function(gene_up,gene_down,pro='test'){
  gene_up=unique(gene_up)
  gene_down=unique(gene_down)
  gene_diff=unique(c(gene_up,gene_down))
  g_list=list(gene_up=gene_up,
              gene_down=gene_down,
              gene_diff=gene_diff)
  # 因?yàn)間o數(shù)據(jù)庫(kù)非常多基因集，所以運(yùn)行速度會(huì)很慢。
  if(T){
    go_enrich_results <- lapply( g_list , function(gene) {
      lapply( c('BP','MF','CC') , function(ont) {
        cat(paste('Now process ',ont ))
        ego <- enrichGO(gene          = gene,
                        #universe      = gene_all,
                        OrgDb         = org.Hs.eg.db,
                        ont           = ont ,
                        pAdjustMethod = "BH",
                        pvalueCutoff  = 0.99,
                        qvalueCutoff  = 0.99,
                        readable      = TRUE)
        
        print( head(ego) )
        return(ego)
      })
    })
    save(go_enrich_results,file =paste0(pro, '_go_enrich_results.Rdata'))
    
  }
  # 只有第一次需要運(yùn)行，就保存結(jié)果哈，下次需要探索結(jié)果，就載入即可。
  load(file=paste0(pro, '_go_enrich_results.Rdata'))
  
  n1= c('gene_up','gene_down','gene_diff')
  n2= c('BP','MF','CC') 
  for (i in 1:3){
    for (j in 1:3){
      fn=paste0(pro, '_dotplot_',n1[i],'_',n2[j],'.png')
      cat(paste0(fn,'\n'))
      png(fn,res=150,width = 1080)
      print( dotplot(go_enrich_results[[i]][[j]] ))
      dev.off()
    }
  }
  
  
}

kegg_plot <- function(up_kegg,down_kegg){
  dat=rbind(up_kegg,down_kegg)
  colnames(dat)
  dat$pvalue = -log10(dat$pvalue)
  dat$pvalue=dat$pvalue*dat$group 
  
  dat=dat[order(dat$pvalue,decreasing = F),]
  
  g_kegg<- ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)), y=pvalue, fill=group)) + 
    geom_bar(stat="identity") + 
    scale_fill_gradient(low="blue",high="red",guide = FALSE) + 
    scale_x_discrete(name ="Pathway names") +
    scale_y_continuous(name ="log10P-value") +
    coord_flip() + theme_bw()+theme(plot.title = element_text(hjust = 0.5))+
    ggtitle("Pathway Enrichment") 
}

運(yùn)行以上代碼后，就可以使用run_kegg，run_go和kegg_plot函數(shù)了

#直接利用腳本中的代碼對(duì)正、負(fù)相關(guān)基因進(jìn)行kegg富集分析
run_kegg(gene_up,gene_down,pro='test')

# 下面的 run_go 會(huì)比較慢，根據(jù)需要
# run_go(gene_up,gene_down,pro='shRNA')

#對(duì)正相關(guān)基因進(jìn)行富集分析
go <- enrichGO(gene_up, OrgDb = "org.Hs.eg.db", ont="all") 
library(ggplot2)
library(stringr)
p1<-barplot(go, split="ONTOLOGY")+ facet_grid(ONTOLOGY~., scale="free")
go=DOSE::setReadable(go, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
write.csv( go@result,file = 'gene_up_GO_all_barplot.csv')

#對(duì)負(fù)相關(guān)基因進(jìn)行富集分析
go <- enrichGO(gene_down, OrgDb = "org.Hs.eg.db", ont="all") 
p2<-barplot(go, split="ONTOLOGY")+ facet_grid(ONTOLOGY~., scale="free") 
 
go=DOSE::setReadable(go, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
write.csv( go@result,file = 'gene_down_GO_all_barplot.csv')

p1+p2

參考來(lái)源
公眾號(hào)：生信技能樹(shù)
#section 3已更新#「生信技能樹(shù)」單細(xì)胞公開(kāi)課2021_嗶哩嗶哩_bilibili

致謝
I thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

THE END