公开数据库汇总及下载3-、ClinVar、dbVar、dbSNP、RefSeqoncoKB
3.3 肿瘤驱动基因
OncoKB
from playwright.sync_api import Playwright, sync_playwright, expect
import time
import random
def run(playwright: Playwright) -> None:
# 基于列表
finished = []
finished_num = 0
with open("finished.txt", "r", encoding="utf-8") as f:
for line in f:
finished.append(line.strip())
finished_num += 1
print("finished_num", finished_num)
genelist = []
with open("cancerGeneList.tsv", "r", encoding="utf-8") as g:
for line in g:
if not line.startswith("Hugo Symbol"):
lineList = line.strip().split("\t")
if not lineList[0] in finished:
genelist.append(lineList[0])
print("downlaod gene list", genelist)
# 随机打乱列表
random.shuffle(genelist)
browser = playwright.chromium.launch(headless=True)
context = browser.new_context()
# 注入反反爬脚本
page = context.new_page()
page.add_init_script(path="steal.js")
# 基因遍历
for gene in genelist:
page.wait_for_load_state("domcontentloaded")
page.goto("https://www.oncokb.org/gene/{gene}".format(gene=gene), timeout=30000) # 设置一个合理的超时时间，比如30秒
time.sleep(random.randint(12, 20))
try:
page.get_by_role("tab", name="Annotated Alterations").click()
time.sleep(2)
html = page.content()
with open("html/gene/" + gene + ".alterations.html", "w", encoding="utf-8") as f:
f.write(html)
except:
pass
try:
page.get_by_role("tab", name="Therapeutic").click()
time.sleep(2)
html = page.content()
with open("html/gene/" + gene + ".therapeutic.html", "w", encoding="utf-8") as f:
f.write(html)
except:
pass
try:
page.get_by_role("tab", name="Diagnostic").click()
time.sleep(2)
html = page.content()
with open("html/gene/" + gene + ".diagnostic.html", "w", encoding="utf-8") as f:
f.write(html)
except:
pass
try:
page.get_by_role("tab", name="Prognostic").click()
time.sleep(2)
html = page.content()
with open("html/gene/" + gene + ".prognostic.html", "w", encoding="utf-8") as f:
f.write(html)
except:
pass
try:
page.get_by_role("tab", name="FDA-Recognized Content").click()
time.sleep(2)
html = page.content()
with open("html/gene/" + gene + ".fda.html", "w", encoding="utf-8") as f:
f.write(html)
except:
pass
html = page.content()
with open("html/gene/" + gene + ".empty.html", "w", encoding="utf-8") as f:
f.write(html)
time.sleep(0.3)
print("【已完成】", gene)
with open("finished.txt", "a", encoding="utf-8") as f:
f.write(gene + "\n")
# 加大爬取间隔，防止被禁
time.sleep(random.randint(20, 30))
page.close()
context.close()
browser.close()
with sync_playwright() as playwright:
run(playwright)
import os
from bs4 import BeautifulSoup
# 基因基本信息
def gene_base_info(gene):
htmlFile = open("html/gene/{gene}.empty.html".format(gene=gene), "r", encoding="utf-8")
htmlHandle = htmlFile.read()
soup = BeautifulSoup(htmlHandle, "html.parser")
div_ele = soup.find_all("div", role="alert")
negative = False
for d in div_ele:
if "We do not have any information for this gene" in d.text:
negative = True
if negative:
outputList = [gene, "无此基因记录", "", "", "", "", ""]
else:
geneInfoTable = soup.find_all("table")[0]
tds = geneInfoTable.find_all("td")
tdDict = {}
for n in range(len(tds)):
if n % 2 == 0:
tdDict[tds[n].text] = tds[n+1].text
ncbi = tdDict["NCBI Gene"]
if "Ensembl Gene" in tdDict:
embl = tdDict["Ensembl Gene"].split(" (")[0]
else:
embl = "-"
grch37 = grch38 = "-"
if "Location" in tdDict:
if "GRch37" in tdDict["Location"]:
grch37 = tdDict["Location"].split(" (GRch37)")[0].replace("Chr", "chr")
if "GRch38" in tdDict["Location"]:
if "GRch37" in tdDict["Location"]:
grch38 = tdDict["Location"].split(" (GRch37)")[1].rstrip(" (GRch38)").replace("Chr", "chr")
else:
grch38 = tdDict["Location"].split(" (GRch38)")[0].replace("Chr", "chr")
if "Ensembl Transcript" in tdDict:
embl_t = tdDict["Ensembl Transcript"].split(" (")[0]
else:
embl_t = "-"
if "RefSeq" in tdDict:
refseq = tdDict["RefSeq"].split(" (")[0]
else:
refseq = "-"
outputList = [gene, ncbi, embl, grch37, grch38, embl_t, refseq]
htmlFile.close()
return outputList
# alterations
def alterations_info(gene):
htmlFile = open("html/gene/{gene}.alterations.html".format(gene=gene), "r", encoding="utf-8")
htmlHandle = htmlFile.read()
soup = BeautifulSoup(htmlHandle, "html.parser")
rows = soup.find_all("div", role="row")
outputList = []
for i in rows[1:]:
divs = i.find_all("div")
outputList.append([gene, divs[0].text, divs[1].text, divs[2].text])
htmlFile.close()
return outputList
# diagnostic
def diagnostic_info(gene):
htmlFile = open("html/gene/{gene}.diagnostic.html".format(gene=gene), "r", encoding="utf-8")
htmlHandle = htmlFile.read()
soup = BeautifulSoup(htmlHandle, "html.parser")
rows = soup.find_all("div", role="row")
outputList = []
for i in rows[1:]:
divs = i.find_all("div")
level = divs[0].span.i.get("class")[-1].lstrip("level-")
output = [gene, level, divs[2].text, divs[3].text]
outputList.append(output)
htmlFile.close()
return outputList
# therapeutic
def therapeutic_info(gene):
htmlFile = open("html/gene/{gene}.therapeutic.html".format(gene=gene), "r", encoding="utf-8")
htmlHandle = htmlFile.read()
soup = BeautifulSoup(htmlHandle, "html.parser")
rows = soup.find_all("div", role="row")
outputList = []
for i in rows[1:]:
divs = i.find_all("div")
level = divs[0].span.i.get("class")[-1].lstrip("level-")
output = [gene, level, divs[2].text, divs[3].text, divs[4].text]
outputList.append(output)
htmlFile.close()
return outputList
# prognostic
def prognostic_info(gene):
htmlFile = open("html/gene/{gene}.prognostic.html".format(gene=gene), "r", encoding="utf-8")
htmlHandle = htmlFile.read()
soup = BeautifulSoup(htmlHandle, "html.parser")
rows = soup.find_all("div", role="row")
outputList = []
for i in rows[1:]:
divs = i.find_all("div")
level = divs[0].span.i.get("class")[-1].lstrip("level-")
output = [gene, level, divs[2].text, divs[3].text]
outputList.append(output)
htmlFile.close()
return outputList
# fda
def fda_info(gene):
htmlFile = open("html/gene/{gene}.fda.html".format(gene=gene), "r", encoding="utf-8")
htmlHandle = htmlFile.read()
soup = BeautifulSoup(htmlHandle, "html.parser")
rows = soup.find_all("div", role="row")
outputList = []
for i in rows[1:]:
divs = i.find_all("div")
output = [gene, divs[1].text, divs[2].text, divs[3].text]
outputList.append(output)
return outputList
####################
geneList = []
with open("oncokb.all.txt", "r", encoding="utf-8") as g:
for line in g:
geneList.append(line.strip())
output1 = open("output_all/output.summary.txt", "w", encoding="utf-8")
output2 = open("output_all/output.anno.txt", "w", encoding="utf-8")
output3 = open("output_all/output.thera.txt", "w", encoding="utf-8")
output4 = open("output_all/output.diagn.txt", "w", encoding="utf-8")
output5 = open("output_all/output.progn.txt", "w", encoding="utf-8")
output6 = open("output_all/output.fda.txt", "w", encoding="utf-8")
htmlList = os.listdir("html/gene")
for gene in geneList:
print(gene)
for i in ["empty", "alterations", "diagnostic", "therapeutic", "prognostic", "fda"]:
if (gene + "." + i + ".html") in htmlList:
if i == "empty":
output1.write("\t".join(gene_base_info(gene)) + "\n")
elif i == "alterations":
for j in alterations_info(gene):
output2.write("\t".join(j) + "\n")
elif i == "therapeutic":
for j in therapeutic_info(gene):
output3.write("\t".join(j) + "\n")
elif i == "diagnostic":
for j in diagnostic_info(gene):
output4.write("\t".join(j) + "\n")
elif i == "prognostic":
try:
for j in prognostic_info(gene):
output5.write("\t".join(j) + "\n")
except:
continue
elif i == "fda":
try:
for j in fda_info(gene):
output6.write("\t".join(j) + "\n")
except:
continue
output1.close()
output2.close()
output3.close()
output4.close()
output5.close()
output6.close()
3.4 遗传变异数据库
ClinVar
它是一个专门用于收集和存储人类基因变异信息的公共数据库。这些基因变异信息包括变异的类型、变异的临床意义（如是否与疾病相关或对药物的反应）、变异所在的基因等。ClinVar数据库的目标是帮助科研人员和医生更好地理解基因变异对人类健康的影响。
ClinVar数据库的内容主要包括三个部分：变异信息、临床信息和相关证据。变异信息是指变异的类型、变异所在的基因等基本信息。临床信息是指变异的临床意义，包括是否与疾病相关、是否对药物反应等。相关证据是指支持临床信息的研究数据或文献。
Clinvar核心是与临床相关的突变，最重要的应用是从个人基因组的海量变异位点中寻找致病或可能致病（P/LP）的位点。**OMIM核心主要集中在疾病与基因层面，**适合在网站上对单个疾病或基因进行检索。
数据库有多种类型xml、vcf等。我们主要用在各个注释软件（VEP、annovar）中对varidct等软件拷出的snv、indel进行注释，标注致病性等信息。
下载链接
wget ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh37/clinvar_20240825.vcf.gz
wget ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh37/clinvar_20240825.vcf.gz.md5
wget ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh37/clinvar_20240825.vcf.gz.tbi
# 具体可都说明文档
wget https://ftp.ncbi.nlm.nih.gov/pub/clinvar/README_VCF.txt
参考文件：
ClinVar数据库介绍
人类遗传变异神库 | ClinVar数据库详解
dbVar
数据库简介
DBVar数据库是一个用于存储和管理人类基因组中的结构变异和重复序列变异信息的数据库。数据库将这些变异信息进行整合、注释和分类，并提供丰富的功能和工具，以便用户能够方便地浏览、搜索和分析这些变异数据。
染色体结构变异structural variation(SV), 被定义为1kb以上范围的DNA结构变化，通常包括缺失，重复，倒位，易位，当然也包含拷贝数变异(CNV), 染色体的结构变异来源于内因或者外因的作用，很多疾病就是由于染色体结构变异导致的，环境因素也可以诱导产生结构变异。
在该数据库中，数据分为以下3个层级:
study,用std表示, 代表一个研究项目，该项目包含了一组结构变异分析结果，每个study有一个唯一的编号，如果来自NCBI，则以nstd开头，如果来自EBI,则以estd开头
variant regions, 用sv表示，代表存在结构变异的基因组区域，同样有一个唯一的编号，如果来自NCBI，则以nsv开头，如果来自EBI,则以esv开头
variant calls, 用ssv表示，代表一个具体的结构变异事件，包含了缺失，重复，CNV等多种类型，同样有一个唯一的编号，如果来自NCBI，则以nssv开头，如果来自EBI,则以essv开头
简单使用：利用基因、位置等信息，查询该位置真实世界的SV信息
下载方法：
wget -r -np -nH --reject "index.html*" -e robots=off ftp://ftp.ncbi.nlm.nih.gov/pub/dbVar/data/Homo_sapiens/
GRCh37.variant_call.all.vcf.gz:
包含了所有类型的变异调用（variant calls），无论其临床意义如何。这是一个综合性的变异数据集，涵盖了 dbVar 数据库中所有结构变异的调用。
GRCh37.variant_call.clinical.benign_or_likely_benign.vcf.gz:
包含所有被认为是良性或可能良性的临床变异。变异被标记为“良性”或“可能良性”意味着这些变异预计不会对健康产生不利影响。
GRCh37.variant_call.clinical.conflicting_or_other.vcf.gz:
包含临床意义存在冲突或其他未归类的变异。这些变异的注释来自不同的来源，它们可能对于相同的变异有不同的临床意义解释。
GRCh37.variant_call.clinical.pathogenic_or_likely_pathogenic.vcf.gz:
**包含所有被认为是致病或可能致病的临床变异**。这些变异可能与某些疾病或健康状况有关。
GRCh37.variant_call.clinical.uncertain_significance.vcf.gz:
包含临床意义不确定的变异。即，这些变异的健康或疾病影响尚不明确，通常需要进一步的研究来确定其临床意义。
GRCh37.variant_call.somatic.vcf.gz:
包含所有体细胞变异（somatic variants）。体细胞变异是发生在个体体细胞中的变异，而不是生殖细胞变异，它们通常与癌症或其他疾病有关。
GRCh37.variant_region.all.vcf.gz:
包含所有类型的变异区域数据。与“variant call”文件相比，变异区域数据通常描述了变异所在的染色体区段，而不是具体的变异事件。
GRCh37.variant_region.somatic.vcf.gz:
包含所有体细胞变异的区域数据。描述了体细胞变异所在的染色体区段，可能与癌症或其他体细胞疾病相关。
参考文件：
dbvar:染色体结构变异数据库
dbSNP
dbsnp是NCBI于1998年建立的主要存储单核苷酸多态性（SNP）的免费公共数据库。该数据库包含多种模式生物。虽然其名称为dbSNP，但该数据库实际上包括多种分子变异：
- 单核苷酸多态性 SNP
- 短缺失和插入多态性 short deletion and insertion polymorphisms (indels/DIPs)
- 微卫星标记或短串联重复 microsatellite markers or short tandem repeats (STRs)
- 多核苷酸多态性 multinucleotide polymorphisms (MNPs)
- 杂合序列 heterozygous sequences
- 命名变体 named variants
首先是snp的基础信息，包括物种，位置，等位基因，变异类型，频率等等
通过切换不同tab可以看到与该snp相关的其他详细信息
下载数据
全部下载的命令行：
wget -r -np -nH --reject "index.html*" -e robots=off ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh37p13/
提供了common和All 两种，All包含所有的SNP位点，common只包含了MAF大于0.01的生殖细胞变异位点,通常-All.vcf.gz。
注意把对应的md5和tbi文件一起下载，md5用于检测下载文件的完整性；tbi文件是vcf文件的索引，方便gatk等程序读取。
此文件通常在对bam文件Realn、BQSR时使用，可通过panel的bed对下载的文件进行过滤，只保留bed范围内的snv或indel突变。
在对突变注释、过滤过程中，由于在正常人中出现的高频突变通常不会是致病位点，因此我们一般注释并过滤掉在dbSNP中已经收录的高频突变位点，保留未被dbSNP收录的突变位点作为候选致病位点，再进行下一步分析。
四大类！外显子组测序中的注释数据库该怎么选？ | 外显子组数据那些事
哪里有人类SNP位点的数据库？
RefSeqGene(人类基因特异性参考基因组序列的集合)
其实，常用的gene注释有不同的来源，这个来源一般是某一个组织通过一定的方法来确定下来的参考gene的相关注释信息。比如常用的有：
RefSeq Gene注释，对gene的不同转录本进行注释，1个转录本对应1个编号成为RefSeq id，例如对于可以翻译成蛋白的转录本，都会以NM_开头如NM_015658；对于不能翻译的转录本，都会以NR_开头如NR_027055；
Ensembl注释；对gene的不同转录本进行注释，以ENSG开头的表示Ensembl gene_id如ENSG00000227232，以ENST开头的表示Ensembl transcript id如ENST00000438504
UCSC gene注释；对gene的不同转录本进行注释，一般是类似uc004cpf这样的名称。
针对1个参考基因组版本，比如human的hg19参考基因组版本，会有来自各种不同组织，不同方法的基因注释。这些基因注释各有优劣，没有绝对的好与绝对的坏，只有掌握他们的注释规则，才能找到最适合我们课题的基因注释。所以，我们今天先来介绍一下RefSeq Gene的注释规则与原理。
RefSeq = The NCBI Reference Sequence
这个计划是由NCBI提出的，意图是为所有常见生物提供非冗余，人工选择过的参考序列。一个物种的RefSeq注释通常包含：参考基因组，参考转录组，参考蛋白序列，参考SNP信息，参考CNV信息等等。
RefSeq基本上我们最常用的一个gene注释版本了，因为经过了人工挑选，挑选出来的gene或者是转录本都十分可靠。
1.人类基因组各基因的序列
# NG_034012.1 Homo sapiens midasin AAA ATPase 1 (MDN1), RefSeqGene on chromosome 6
https://ftp.ncbi.nih.gov/refseq/H_sapiens/RefSeqGene/
2.gtf
# NCBI不断更新人类基因组版本,里面有注释文件
https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.25_GRCh37.p13/
参考文件：
生物信息学100个基础问题 —— 第26题 什么是RefSeq Gene?
基因组注释文件(GFF,GTF)下载的五种方法