前言
目前基于RNA做分析的文章中几乎都有 GSEA 的分析内容，并且都是出现在正文，那么这个也是表达基因筛选的一种重要方式，下面我将整个流程梳理一下，供大家参考。
GSEA（Gene Set EnrichmentAnalysis），即基因集富集分析，它的基本思想是使用预定义的基因，将基因按照在两类样本中的差异表达程度排序，然后检验预先设定的基因集合是否在这个排序表的顶端或者底端富集。
GSEA和GO、KEGG pathway不同的地方在于，后两者会提前设定一个阈值，只关注差异变化大的基因（相当于重点班）。这样子容易遗漏部分差异表达不显著却有重要生物学意义的基因（成绩一般，但是很有天赋）。所以GSEA分析比较适用于，传统分析方法筛选后样本过少的数据集。
GSEA数据库收集了很多分子标记数据，有9大分类的基因，如下：
九大分类如下：
实例解读
1. 数据读取
数据的读取我们仍然使用的是 TCGA-COAD 的数据集，表达数据的读取以及临床信息分组的获得我们上期已经提过，我们使用的是edgeR 软件包计算出来的差异表达结果，提取上调基因 2832 的 ENSEMBL 号，
###########基因列表
DEG=read.table("DEG-resdata.xls",sep="\t",check.names=F,header = T)
geneList<-DEG[DEG$sig=="Up",]$Row.names
table(DEG$sig)##
## Down Up
## 1296 2832直接用之前的转录组差异分析后的数据来演示，数据格式如下：
head(DEG[,1:7])## Row.names logFC logCPM LR PValue
## 1 ENSG00000142959 -5.924830 3.20085238 1457.1137 8.159939e-319
## 2 ENSG00000163815 -4.223787 2.59792097 1145.7948 3.679815e-251
## 3 ENSG00000107611 -5.131288 1.71477711 1122.0643 5.289173e-246
## 4 ENSG00000162461 -4.101967 1.51480635 1085.9423 3.752089e-238
## 5 ENSG00000163959 -4.295806 3.39558390 1080.7407 5.067873e-237
## 6 ENSG00000144410 -6.284258 -0.02616284 916.3497 2.739233e-201
## FDR sig
## 1 2.592984e-314 Down
## 2 5.846674e-247 Down
## 3 5.602468e-242 Down
## 4 2.980753e-234 Down
## 5 3.220836e-233 Down
## 6 1.450744e-197 Down因为GSEA只需要SYMBOL（基因名）和foldchange （或logFC)两列，所以我们提取这两列，如下：
geneList<-DEG[,1:2]
head(geneList)## Row.names logFC
## 1 ENSG00000142959 -5.924830
## 2 ENSG00000163815 -4.223787
## 3 ENSG00000107611 -5.131288
## 4 ENSG00000162461 -4.101967
## 5 ENSG00000163959 -4.295806
## 6 ENSG00000144410 -6.2842582. GSEA富集结果
首先我们同样需要安装软件包并加载，这里面主程序就是 clusterProfiler 软件包，如下：
########
if(!require(clusterProfiler)){
BiocManager::install("clusterProfiler")
}
if(!require(org.Hs.eg.db)){
BiocManager::install("org.Hs.eg.db")
}
if(!require(DOSE)){
BiocManager::install("DOSE")
}
if(!require(topGO)){
BiocManager::install("topGO")
}
if(!require(pathview)){
BiocManager::install("pathview")
}
if(!require(KEGG.db)){
BiocManager::install("KEGG.db")
}
if(!require(enrichplot)){
BiocManager::install("enrichplot")
}
library(org.Hs.eg.db)
library(clusterProfiler)
library(DOSE)
library(topGO)
library(pathview)
library(KEGG.db)
library(enrichplot)如基因名是symbol，需要将基因ID转换为Entrez ID格式。Entrez ID实际上是指的Entrez gene ID，是对应于染色体上一个gene location的。每一个发现的基因都会被编制一个统一的编号，而Entrez ID是指的来自于NCBI旗下的Entrez gene数据库所使用的编号。因为Entrez ID具有特异性，所以后续分析更适合用Entrez ID
eg <- bitr(geneList$Row.names,
fromType="ENSEMBL",
toType=c("ENTREZID","ENSEMBL",'SYMBOL'),
OrgDb="org.Hs.eg.db")## 'select()' returned 1:many mapping between keys and columns## Warning in bitr(geneList$Row.names, fromType = "ENSEMBL", toType =
## c("ENTREZID", : 30.89% of input gene IDs are fail to map...head(eg)## ENSEMBL ENTREZID SYMBOL
## 1 ENSG00000142959 266675 BEST4
## 2 ENSG00000163815 7123 CLEC3B
## 3 ENSG00000107611 8029 CUBN
## 4 ENSG00000162461 284723 SLC25A34
## 5 ENSG00000163959 200931 SLC51A
## 6 ENSG00000144410 130749 CPO合并基因列表，就是我们需要输入的文件了，如下：
mergedata<-merge(eg,geneList,by.x="ENSEMBL",by.y="Row.names")
head(mergedata)## ENSEMBL ENTREZID SYMBOL logFC
## 1 ENSG00000003249 79007 DBNDD1 2.474017
## 2 ENSG00000004776 126393 HSPB6 -3.543461
## 3 ENSG00000004799 5166 PDK4 -3.281233
## 4 ENSG00000004846 340273 ABCB5 -3.084616
## 5 ENSG00000005001 64063 PRSS22 4.095319
## 6 ENSG00000005981 51666 ASB4 2.237660mergedata_sort <- mergedata[order(mergedata$logFC, decreasing = T),]#先按照logFC降序排序
gene_fc = mergedata_sort$logFC #把foldchange按照从大到小提取出来
names(gene_fc) <- mergedata_sort$ENTREZID #给上面提取的foldchange加上对应上ENTREZID
head(gene_fc)## 4102 26831 26832 5653 646960 4105
## 10.010124 9.817905 9.807069 9.553685 9.448116 9.440482GSEA 分析
需要gmt文件，http://www.gsea-msigdb.org/gsea/downloads.jsp 路径下载，选择合适的，我们这里是结直肠癌中癌和癌旁的差异比较分析，所以我们需要找到致癌基因，选择 c5: Ontology gene sets 即 c5.all.v7.5.1.entrez.gmt 文件，然后进行 GSEA 富集分析，如下：
library(stats)
kegg_gmt <- read.gmt("c5.all.v7.5.1.entrez.gmt") #读gmt文件
gsea <- GSEA(gene_fc,
TERM2GENE = kegg_gmt) #GSEA分析
head(gsea)## ID
## GOBP_ADAPTIVE_IMMUNE_RESPONSE GOBP_ADAPTIVE_IMMUNE_RESPONSE
## GOBP_CATION_TRANSPORT GOBP_CATION_TRANSPORT
## GOBP_IMMUNE_RESPONSE_REGULATING_SIGNALING_PATHWAY GOBP_IMMUNE_RESPONSE_REGULATING_SIGNALING_PATHWAY
## GOBP_IMMUNOGLOBULIN_PRODUCTION GOBP_IMMUNOGLOBULIN_PRODUCTION
## GOBP_ION_TRANSPORT GOBP_ION_TRANSPORT
## GOBP_LIPID_METABOLIC_PROCESS GOBP_LIPID_METABOLIC_PROCESS
## Description
## GOBP_ADAPTIVE_IMMUNE_RESPONSE GOBP_ADAPTIVE_IMMUNE_RESPONSE
## GOBP_CATION_TRANSPORT GOBP_CATION_TRANSPORT
## GOBP_IMMUNE_RESPONSE_REGULATING_SIGNALING_PATHWAY GOBP_IMMUNE_RESPONSE_REGULATING_SIGNALING_PATHWAY
## GOBP_IMMUNOGLOBULIN_PRODUCTION GOBP_IMMUNOGLOBULIN_PRODUCTION
## GOBP_ION_TRANSPORT GOBP_ION_TRANSPORT
## GOBP_LIPID_METABOLIC_PROCESS GOBP_LIPID_METABOLIC_PROCESS
## setSize enrichmentScore
## GOBP_ADAPTIVE_IMMUNE_RESPONSE 111 -0.4743354
## GOBP_CATION_TRANSPORT 185 -0.2838816
## GOBP_IMMUNE_RESPONSE_REGULATING_SIGNALING_PATHWAY 52 -0.4982248
## GOBP_IMMUNOGLOBULIN_PRODUCTION 51 -0.6152482
## GOBP_ION_TRANSPORT 252 -0.2763134
## GOBP_LIPID_METABOLIC_PROCESS 163 -0.2922133
## NES pvalue
## GOBP_ADAPTIVE_IMMUNE_RESPONSE -4.455485 1e-10
## GOBP_CATION_TRANSPORT -2.957025 1e-10
## GOBP_IMMUNE_RESPONSE_REGULATING_SIGNALING_PATHWAY -3.636516 1e-10
## GOBP_IMMUNOGLOBULIN_PRODUCTION -4.430346 1e-10
## GOBP_ION_TRANSPORT -2.791438 1e-10
## GOBP_LIPID_METABOLIC_PROCESS -2.998976 1e-10
## p.adjust
## GOBP_ADAPTIVE_IMMUNE_RESPONSE 1.548421e-08
## GOBP_CATION_TRANSPORT 1.548421e-08
## GOBP_IMMUNE_RESPONSE_REGULATING_SIGNALING_PATHWAY 1.548421e-08
## GOBP_IMMUNOGLOBULIN_PRODUCTION 1.548421e-08
## GOBP_ION_TRANSPORT 1.548421e-08
## GOBP_LIPID_METABOLIC_PROCESS 1.548421e-08
## qvalues rank
## GOBP_ADAPTIVE_IMMUNE_RESPONSE 1.040997e-08 1137
## GOBP_CATION_TRANSPORT 1.040997e-08 1082
## GOBP_IMMUNE_RESPONSE_REGULATING_SIGNALING_PATHWAY 1.040997e-08 1130
## GOBP_IMMUNOGLOBULIN_PRODUCTION 1.040997e-08 1137
## GOBP_ION_TRANSPORT 1.040997e-08 852
## GOBP_LIPID_METABOLIC_PROCESS 1.040997e-08 916
## leading_edge
## GOBP_ADAPTIVE_IMMUNE_RESPONSE tags=86%, list=40%, signal=54%
## GOBP_CATION_TRANSPORT tags=64%, list=38%, signal=42%
## GOBP_IMMUNE_RESPONSE_REGULATING_SIGNALING_PATHWAY tags=88%, list=39%, signal=55%
## GOBP_IMMUNOGLOBULIN_PRODUCTION tags=100%, list=40%, signal=61%
## GOBP_ION_TRANSPORT tags=51%, list=30%, signal=40%
## GOBP_LIPID_METABOLIC_PROCESS tags=58%, list=32%, signal=42%
## core_enrichment
## GOBP_ADAPTIVE_IMMUNE_RESPONSE 28885/28432/28420/28943/28834/28908/6977/28803/28876/28785/28820/28903/28781/4063/28442/729230/28784/28815/974/28892/3514/28775/28409/28388/28774/11126/28900/28941/28518/6363/28823/28449/28448/28923/28412/660/5579/28804/28392/30835/28938/1191/28811/28912/28423/634/28935/3662/28786/28444/28464/28450/973/28821/28937/28875/28874/28915/1380/28921/10990/3493/28385/28401/10462/28906/57379/28408/3570/28893/28772/28817/28776/28780/2208/10332/28896/28930/28457/339390/28516/28883/83417/28410/28940/28881/28870/608/79908/23495/28877/730/3494/28919/3512/28452
## GOBP_CATION_TRANSPORT 3773/79901/1756/9472/6328/493/2316/6549/729230/285195/6263/6330/2852/114780/3738/387849/2668/245972/170850/3782/9058/130497/6363/492/57113/5023/5579/482/170572/2824/143425/2903/534/4129/3359/6355/1179/5174/169026/931/6568/6616/84239/84502/844/1272/1393/53826/5649/6387/3741/146395/1261/148/9992/6446/6572/6543/3777/6564/57030/3778/57282/162514/287/57158/6860/9963/341359/140738/5348/3779/4842/2273/846/4638/285242/5350/1908/3736/2641/11280/6335/10050/3769/22871/1258/7432/1268/3745/6327/115584/22953/6865/6550/159963/60482/117/644139/2823/6338/140803/845/6340/8671/6332/286133/477/1804/201140/246213/5354/340024/6561/55532/347741/6555/92736
## GOBP_IMMUNE_RESPONSE_REGULATING_SIGNALING_PATHWAY 28432/28420/2206/28834/7098/28442/974/3514/28409/81793/28388/11126/26228/28449/28448/28412/660/5579/28392/148823/948/23228/28423/634/933/640/3662/28444/28464/28450/973/931/1380/3493/28385/28401/28408/64922/10917/11148/28457/28410/79908/3494/9971/28452
## GOBP_IMMUNOGLOBULIN_PRODUCTION 28943/51237/28908/28803/28876/28785/28820/28903/28781/28784/28815/28892/3514/28775/28774/28900/28941/28518/28823/28923/29802/28804/28938/28811/28912/933/28935/28786/28821/28937/28875/28874/28915/28921/28906/57379/28893/28772/28817/28776/28780/28896/28930/28516/28883/28940/28881/28870/28877/28919
## GOBP_ION_TRANSPORT 57113/220963/5023/5579/2350/482/170572/2824/5322/143425/2903/948/534/634/4129/3359/6355/1179/5028/5174/169026/931/6568/6616/84239/84502/844/1272/1393/53826/5649/2893/6387/3741/146395/387700/1261/148/9992/64805/28965/6446/55065/2565/6572/2901/6543/3777/1188/6564/2566/57030/3778/57282/162514/391013/287/57158/6860/9963/203859/341359/140738/5348/3779/4842/2273/846/4638/285242/5350/1908/3736/2641/11280/6335/10050/3769/22871/553/5030/2899/1258/7432/353189/1268/3745/6327/115584/22953/6865/54831/6550/159963/60482/117/5593/644139/11136/2823/123264/1836/6338/140803/200931/845/6340/8671/1811/6332/760/286133/477/9429/1804/8647/201140/762/6514/246213/5354/22802/340024/6561/766/55532/266675/347741/6555/92736
## GOBP_LIPID_METABOLIC_PROCESS 6817/5288/124/1562/6363/5959/283985/5408/5023/660/389396/5333/2494/5322/1191/3158/948/23228/79154/634/3291/345557/1565/129807/84803/10690/54988/54657/653/30815/10924/10317/64805/28965/23584/5105/159296/137872/2168/3248/1580/54511/80168/5742/1907/6822/54658/51046/13/391013/57016/339221/374569/84830/27306/2169/7366/8708/10170/54857/5166/53345/54979/9227/3294/133121/1268/185/54575/9971/6799/7148/80157/7367/126/124872/51059/6338/2938/54576/125/10351/1103/57733/1576/8647/8029/5354/3284/2538/338/335/345/337-
ID ：信号通路
-
Description ：信号通路的描述
-
setSize ：富集到该信号通路的基因个数
-
enrichmentScore ：富集分数，也就是ES
-
NES ：标准化以后的ES，全称normalized enrichment score
-
pvalue：富集的P值
-
p.adjust ：校正后的P值
-
qvalues ：FDR （false discovery rate）错误发现率
-
rank ：排名
-
core_enrichment：富集到该通路的基因列表。
绘制气泡图
气泡图解读需要说明一下，富集程度通过Gene ratio、Pvalue和富集到此基因集上的基因个数来衡量。
-
横坐标是Gene ratio，数值越大表示富集程度越大。Count 位于该基因集下的差异表达基因数
-
纵坐标是富集程度较高的基因集（一般选取富集最显著的20条进行展示，不足20条则全部列出）。
-
Pvalue取值范围[0, 1]，以颜色表示，越红表示Pvalue越小，说明富集越明显。
-
点的大小表示该基因集下差异基因的个数，点越大表示基因数越多**。**
气泡图绘制，如下：
dotplot(gsea)
GO 富集分析
# Run GO and KEGG enrichment analysis
#GO富集
GO <- gseGO(
gene_fc, #gene_fc
ont = "BP",# "BP"、"MF"和"CC"或"ALL"
OrgDb = org.Hs.eg.db,#人类注释基因
keyType = "ENTREZID",
pvalueCutoff = 0.05,
pAdjustMethod = "BH",#p值校正方法
)
head(GO[,1:4])## ID Description
## GO:0002250 GO:0002250 adaptive immune response
## GO:0002253 GO:0002253 activation of immune response
## GO:0002377 GO:0002377 immunoglobulin production
## GO:0002764 GO:0002764 immune response-regulating signaling pathway
## GO:0006139 GO:0006139 nucleobase-containing compound metabolic process
## GO:0006396 GO:0006396 RNA processing
## setSize enrichmentScore
## GO:0002250 110 -0.4726934
## GO:0002253 50 -0.5141143
## GO:0002377 51 -0.6152482
## GO:0002764 52 -0.4982248
## GO:0006139 443 0.3195185
## GO:0006396 105 0.4952210可以进行一些调整以接近文献
1）修改GSEA线条颜色
2）添加P值的table
3）展示指定的GO或pathway
4）展示多个GO或pathway
5）只展示上两部
展示指定的通路，绘制GSEA 图形，这里选择我们注释的一个结果 "GO:0002250"，如下：
gseaplot2(GO, "GO:0002250", color = "firebrick", rel_heights=c(1, .2, .6))
展示多个GO， 展示多个GSEA - GO 富集结果，添加P值的table，如下：
go<-GO$ID[1:4]
gseaplot2(GO,
go,
pvalue_table = TRUE,
color = colorspace::rainbow_hcl(4),
base_size = 10)
不添加P值的table，如下：
gseaplot2(GO,
go,
pvalue_table = FALSE,
color = colorspace::rainbow_hcl(4),
base_size = 10)
只展示上两部，如下：
gseaplot2(GO,go,
color = colorspace::rainbow_hcl(4),
subplots=c(1,2),
pvalue_table = FALSE)
KEGG 富集分析
#KEGG富集
KEGG<-gseKEGG(
gene_fc,
organism = "hsa",
keyType = "kegg",
exponent = 1,
minGSSize = 10,
maxGSSize = 500,
eps = 1e-10,
pvalueCutoff = 0.05,
pAdjustMethod = "BH",
verbose = TRUE,
use_internal_data = FALSE,
seed = FALSE,
by = "fgsea"
)
sortKEGG<-KEGG[order(KEGG$enrichmentScore, decreasing = T),]#按照enrichment score从高到低排序
head(sortKEGG)## ID Description setSize
## hsa04613 hsa04613 Neutrophil extracellular trap formation 36
## hsa05322 hsa05322 Systemic lupus erythematosus 35
## hsa05203 hsa05203 Viral carcinogenesis 25
## hsa05034 hsa05034 Alcoholism 48
## hsa04080 hsa04080 Neuroactive ligand-receptor interaction 81
## hsa01100 hsa01100 Metabolic pathways 143
## enrichmentScore NES pvalue p.adjust qvalues
## hsa04613 0.5819289 2.796643 6.536810e-08 3.819955e-06 2.637779e-06
## hsa05322 0.5712915 2.717399 3.160108e-07 9.547559e-06 6.592841e-06
## hsa05203 0.5117522 2.204931 3.377795e-04 3.247880e-03 2.242746e-03
## hsa05034 0.4554528 2.402982 1.454689e-05 2.272951e-04 1.569533e-04
## hsa04080 -0.2327878 -1.935807 1.679454e-03 1.208257e-02 8.343336e-03
## hsa01100 -0.3054455 -3.049852 3.470165e-10 4.337706e-08 2.995300e-08
## rank leading_edge
## hsa04613 731 tags=72%, list=25%, signal=55%
## hsa05322 731 tags=71%, list=25%, signal=54%
## hsa05203 520 tags=52%, list=18%, signal=43%
## hsa05034 865 tags=62%, list=30%, signal=44%
## hsa04080 792 tags=48%, list=28%, signal=36%
## hsa01100 933 tags=59%, list=32%, signal=42%
## core_enrichment
## hsa04613 8368/8346/8364/3018/8366/8360/8335/85235/8354/317772/8356/8331/8332/8367/8348/8342/8359/8344/8968/8340/8343/5582/2244/8336/8350/440689
## hsa05322 8368/8346/8364/3018/8366/8360/8335/85235/8354/317772/8356/8331/8332/8367/8348/8342/8359/8344/8968/8340/8343/2904/8336/8350/440689
## hsa05203 8368/8346/8364/3018/8366/8360/8367/8348/8342/8359/8344/8340/8343
## hsa05034 8368/8346/8364/3018/8366/8360/8335/85235/8354/317772/8356/8331/8332/8367/8348/2906/8342/8359/8344/8968/8340/8343/2904/8336/8350/440689/7054/2792/1813/3013
## hsa04080 3363/2903/6752/7068/155/5028/66004/1816/2893/9934/148/84539/2565/2901/3360/1907/2566/5179/1511/1908/2641/2587/553/5030/1129/2899/4852/7432/1268/185/22953/6865/9340/117/5539/4887/6750/5697/10022
## hsa01100 51301/79153/5288/124/1562/771/5959/5408/5333/2878/377841/5322/7498/3158/534/4129/3291/129807/84803/10690/54988/11181/54657/5152/957/30815/5138/10317/8639/23382/57452/83539/5105/54511/80168/5742/25834/2949/111/54658/2875/197257/391013/57016/339221/27306/7366/4842/64902/79799/6476/60495/8708/10170/9227/3294/133121/5145/7166/229/54575/2651/290/6898/443/5837/270/7367/126/2938/54576/125/248/150572/760/57733/1576/762/8972/2939/1610/3284/766/2538/759展示指定通路，如下：
gseaplot2(KEGG, "hsa04613", color = "firebrick", rel_heights=c(1, .2, .6))
展示指定多个通路，如下：
paths <- c("hsa04613", "hsa05322", "hsa05203", "hsa05034")#选取你需要展示的通路ID
gseaplot2(KEGG,
paths,
pvalue_table = TRUE,
color = colorspace::rainbow_hcl(4),
base_size = 20)
不添加Pvalue表格
gseaplot2(KEGG,paths,
color = colorspace::rainbow_hcl(4),
pvalue_table = FALSE)
只展示上面两部分图形，如下：
gseaplot2(KEGG,paths,
color = colorspace::rainbow_hcl(4),
subplots=c(1,2),
pvalue_table = TRUE)
可以选择点图 ES_geom = "dot"，如下：
gseaplot2(KEGG,paths,
color = colorspace::rainbow_hcl(4),
pvalue_table = FALSE)
GSEA的结果解读
1、图最上面部分展示的是富集分数（ES,enrichment score）值计算过程，从左至右每到一个基因，计算出一个ES值，连成线。在最左侧或最右侧有一个特别明显的峰值就是基因集表型上的ES值。图中间部分每一条线代表基因集中的一个基因，及其在基因列表中的排序位置。
2、最下面部分展示的是基因与表型关联的矩阵，红色为与第一个表型(class A)正相关，在class A中表达高，蓝色与第二个表型(class B)正相关，在class B中表达高。
3、Leading-edge subset 对富集得分贡献最大的基因成员。若富集得分为正值，则是峰左侧的基因；若富集得分为负值，则是峰右侧的基因。
4、FDR GSEA默认提供所有的分析结果，并且设定FDR<0.25为可信的富集，最可能获得有功能研究价值的结果。但如果样品数目少，而且选择了gene_set作为Permumation type则需要使用更为严格的标准，比如FDR<0.05。
个人感觉总结的已经很全面了，不足的地方请大家指正。
关注公众号，每日更新，扫码进群交流不停歇，马上就出视频版，关注我，您最佳的选择！
桓峰基因 生物信息分析，SCI文章撰写及生物信息基础知识学习：R语言学习，perl基础编程，linux系统命令，Python遇见更好的你 47篇原创内容 --> 公众号
References:
-
Zhuang W, Sun H, Zhang S, et al. An immunogenomic signature for molecular classification in hepatocellular carcinoma. Mol Ther Nucleic Acids. 2021;25:105-115. Published 2021 Jul 2. doi:10.1016/j.omtn.2021.06.024
-
Cao J, Zhang C, Jiang GQ, et al. Identification of hepatocellular carcinoma-related genes associated with macrophage differentiation based on bioinformatics analyses. Bioengineered. 2021;12(1):296-309. doi:10.1080/21655979.2020.1868119
-
Ju A, Tang J, Chen S, Fu Y, Luo Y. Pyroptosis-Related Gene Signatures Can Robustly Diagnose Skin Cutaneous Melanoma and Predict the Prognosis. Front Oncol. 2021;11:709077. Published 2021 Jul 13. doi:10.3389/fonc.2021.709077
-
Tan L, Xu Q, Shi R, Zhang G. Bioinformatics analysis reveals the landscape of immune cell infiltration and immune-related pathways participating in the progression of carotid atherosclerotic plaques. Artif Cells Nanomed Biotechnol. 2021;49(1):96-107. doi:10.1080/21691401.2021.1873798
网友评论