1.写在前面

最近在中国知网刷到一篇关于网络药理学的文章，文章的逻辑和使用的技术还挺简单。记录一下学习流程，为了避免不必要的麻烦，就不用原文举例了。

为了简化接下来的教程，我就以关键词“肝癌”和“中药”在知网搜索，查看有哪些中药可以治疗肝癌。搜索到这篇文章

Screen Shot 2023-12-31 at 08.57.00.png

链接如下：https://kns.cnki.net/kcms2/article/abstract?v=OHSidZhWKUkSVAcEooQbXpIa97MtfCTTNozL3raYBCkQtL2jMZtpJzuXggfz3MR0MjMc87DUBdV6DSf-WGGq6fZEe4rbC8GCIqB_FKeVzdR3ELZ2dLBN7rFND6cc-o-8&uniplatform=NZKPT&flag=copy

接下来的研究对象就有了：
中药：山豆根
疾病：肝癌

2.数据准备

1.TCMSP：https://old.tcmsp-e.com/tcmsp.php

这个数据库主要用于搜索中药成分和中药成分对应的靶点，参考山西医科大学发布的TCMSP教程《7--手把手教你从TCMSP数据库获取中药有效成分和中药靶点》：http://www.sxmu.edu.cn/bdcd/info/1110/1291.htm

1.打开官网，搜索“山豆根”；

Screen Shot 2023-12-31 at 09.06.04.png

2.点击“Sophorae Tonkinensis Radix Et Rhizome”；

Screen Shot 2023-12-31 at 09.07.44.png

3.按照山西医科学发布的教程，选择）OB>30%，DL>0.18，筛选中药中的有效成分。

Screen Shot 2023-12-31 at 09.10.10.png

4.请注意Mol ID这一列，我们要根据这一列下载Related Targets。

TCSMP不好的一点是不能直接下载搜索到的结果，需要一页页的复制到excel里。注意这里的Mol ID是无序的，下载起来会很麻烦，点一Mol ID列可以根据升序进行排列。

接下来下载Related Targets，将Mol ID按照生序排列，然后根据刚刚下载的有效成分的Mol ID下载这个有效成分的Related Targets，依然是一页页的复制到excel里。

素质使然就不要用爬虫去爬了。

Screen Shot 2023-12-31 at 09.20.12.png

5. 结果：我们一共筛选到山豆根的21个有效成分，这个21个有效成分又有394个靶基因。

注意：不是每个有效成分都有靶基因，如果找不到某个有效成分的靶基因直接跳过这个有效成分就可以了。

3.准备肝癌基因相关的数据集

1. genecars: https://www.genecards.org

Screen Shot 2023-12-31 at 09.42.18.png

2.在右上角搜索框以“liver cancer”搜索

Screen Shot 2023-12-31 at 09.44.18.png

3.点击export导出搜索结果。操作此步骤需要注册genecards账户。

Screen Shot 2023-12-31 at 09.50.29.png

到目前为止，所需要的数据文件我们已经准备的差不多了，但需要对下载的文件进行进一步的处理。

3.数据处理

打开中药成分的靶基因文件，在Target name这一列，可以看到显示的是基因的全名，而在我们从genecards导出的肝癌数据显示的是基因的Symbol，因此，我们需要把Target name这一列转换成基因的Symbol。

1.打开uniprot网站：https://www.uniprot.org

Screen Shot 2023-12-31 at 09.58.55.png

2.点击proteins，然后搜索，按照箭头指示，最后导出搜索到的数据就可以了。

Screen Shot 2023-12-31 at 10.01.59.png

3.最终把下载的文件整理成以下格式即可。

Screen Shot 2023-12-31 at 10.21.34.png

4. 在靶基因文件中添加Symbol列，使用R语言完成。

library(tidyverse)
library(dplyr)
library(readxl)

# 读取靶基因文件
Related_Targets <- read_excel("Related Targets.xlsx")
head(Related_Targets)

# 读取基因注释文件
gene_annotation <- read_excel("gene_annotation.xlsx")
head(gene_annotation)

data = merge(Related_Targets, gene_annotation, by.x = "Target name", by.y = "Protein names")
head(data)

# 保存文件
write.table(data, "ingredient_target.txt",row.names = F,sep = "/t", quote = F)

最终整理成的数据格式如下所示：

Screen Shot 2023-12-31 at 10.41.01.png

5.筛选genecards下载的"liver cancer"数据集

一般选择Relevance score大于中间值的基因就可以了，也可以自定义筛选标准。

# 读取genecards数据集
liver_data = read.csv("GeneCards-liver cancer.csv",header = T)
dim(liver_data) # 4493    8

filter_liver_data = filter(liver_data, Relevance.score > median(liver_data$Relevance.score))
dim(filter_liver_data) # 2241    8
# 过滤后一共有2241个基因

# 保存文件
write.table(filter_liver_data, "filter_liver_data.txt", row.names = F,sep = "/t", quote = F)

4.到此为止，我们已经准备好了所有的数据，接下来可以进行分析了。

分析所需的文件：ingredient_target.txt，filter_liver_data.txt

1.首先查看药物成分-靶基因与肝癌基因集中有多少个交集基因，绘制venn图。

# 清空变量
rm(list = ls())
# 读取数据
target = read.table("ingredient_target.txt", header = T, sep = "/t")
liver = read.table("filter_liver_data.txt", header = T, sep = "/t")
# 查看数据集
head(target)
head(liver)

# 查看交集基因
inter_gene = intersect(target$Gene.Names, liver$Gene.Symbol)
length(inter_gene) # 一个有98个基因

# 可视化
library(ggvenn)
gene = list(Ingredient_target_gene = target$Gene.Names, liver_cancaer_gene = liver$Gene.Symbol)
ggvenn(gene,fill_color = c("red","blue"))# 根据自己的喜好更改颜色就可以了

image.png

2.将交集中的基因从Ingredient-targets中挑选出来，使用cytoscape软件绘图

我电脑没装cytoscape软件，就不展示这部分结果了，cytoscape软件使用很简单。

select_ingredient_target = subset(target,Gene.Names %in% inter_gene)
dim(select_ingredient_target) # 165   6
write.table(select_ingredient_target,"select_ingredient_target.txt", row.names = F,sep = "/t", quote = F)
head(select_ingredient_target)
# Target.name    Mol.ID Molecule.name   Source    status Gene.Names
# 7  72 kDa type IV collagenase MOL005944       matrine DrugBank validated       MMP2
# 8  72 kDa type IV collagenase MOL000098     quercetin DrugBank validated       MMP2
# 9    Acetyl-CoA carboxylase 1 MOL000098     quercetin DrugBank validated      ACACA
# 33    Apoptosis regulator BAX MOL000422    kaempferol      N/A validated        BAX
# 34    Apoptosis regulator BAX MOL000098     quercetin      N/A validated        BAX
# 35  Apoptosis regulator Bcl-2 MOL000098     quercetin DrugBank validated       BCL2

3.富集分析

就是将筛选出来的98个基因使用clusterprofiler包进行GO和KEGG富集分析。

GO富集分析

可以将GO富集分析的结果保存下来，挑选合适的GO term展示。
在这里随机展示GO富集分析的结果。

 富集分析
library(clusterProfiler)
library(enrichplot)
library(org.Hs.eg.db)

# 基因ID转换
gene_id = bitr(inter_gene, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = "org.Hs.eg.db")$ENTREZID
head(gene_id)

# go富集分析
go = enrichGO(gene_id,
              OrgDb = "org.Hs.eg.db",
              keyType = "ENTREZID",
              ont = "all",
              pvalueCutoff = 0.05,
              pAdjustMethod = "BH",
              qvalueCutoff = 0.2)
barplot(go,
        showCategory = 5,
        label_format = 50,
        split = "ONTOLOGY")+
  facet_grid(ONTOLOGY~., scale = "free")
# 保存go富集分析结果
write.csv(go,"go.csv")

image.png

KEGG富集分析

以可以挑选合适KEGG富集分析的结果用于展示，选择自己感兴趣的信号通路即可。这里随机展示KEGG富集分析的结果。

kegg = enrichKEGG(gene_id,
                  organism = "hsa",
                  keyType = "kegg",
                  pvalueCutoff = 0.05,
                  pAdjustMethod = "BH",
                  qvalueCutoff = 0.2)
barplot(kegg)
# write.csv(kegg,"kegg.csv")

image.png

写在最后：

经过以上分析，我们根据自己的研究背景应该已经可以得出自己想要研究的药物成分和目的基因了，接下来就是进行分子对接，我电脑里目前还没装这些分子对接所需的软件。这几天我把软件装好以后就写接下来如何进行分子对接。

（我就是看到别人的水文时突发奇想写了这个教程，也不知道算不算教程，思路有点乱，等不忙了再修改吧。）