GRCh38版本染色体位置转换GRCh37(hg19)
目录
- 方法 1:使用 Ensembl REST API(推荐,适用于少量位点查询)
- 方法 2:使用 UCSC API
- 方法 3:使用 NCBI API 并转换坐标(需要额外步骤)
- 方法 4:使用本地数据库(最可靠,适合批量查询)
- 1、下载 GRCh37 基因注释文件:
- 2、创建 Python 查询脚本:
- 其他
- 注意事项
- Ensembl API 限制:
- 基因组版本差异:
- 基因位置表示:
数据准备:
要获取 GRCh37(hg19)版本的染色体位置,我们需要使用专门针对该基因组版本的数据源。
方法 1:使用 Ensembl REST API(推荐,适用于少量位点查询)
import requests
import timedef get_gene_location_hg19(gene_symbol, max_retries=3):"""获取基因在 GRCh37/hg19 版本的位置"""base_url = "http://grch37.rest.ensembl.org"endpoint = f"/lookup/symbol/homo_sapiens/{gene_symbol}"url = base_url + endpointparams = {"expand": 0,"utr": 0,"content-type": "application/json"}for attempt in range(max_retries):try:response = requests.get(url, params=params, timeout=30)response.raise_for_status()data = response.json()chrom = data['seq_region_name']start = data['start']end = data['end']strand = data['strand']# 返回位置信息(添加chr前缀)return f"chr{chrom}:{start}-{end}"except requests.exceptions.RequestException as e:if attempt < max_retries - 1:wait_time = 2 ** attempt # 指数退避print(f"Attempt {attempt+1} failed. Retrying in {wait_time} seconds...")time.sleep(wait_time)else:raise Exception(f"Failed to get location for {gene_symbol}: {str(e)}")# 测试
print(get_gene_location_hg19("TP53")) # 输出: chr17:7571719-7590868
print(get_gene_location_hg19("BRCA1")) # 输出: chr17:41196312-41277500
方法 2:使用 UCSC API
import requests
import xml.etree.ElementTree as ETdef get_gene_location_ucsc_hg19(gene_symbol):"""使用 UCSC API 获取 hg19 位置"""base_url = "https://api.genome.ucsc.edu"endpoint = f"/getData/track"params = {"genome": "hg19","track": "refGene","gene": gene_symbol,"maxItemsOutput": 1}try:response = requests.get(base_url + endpoint, params=params, timeout=30)response.raise_for_status()data = response.json()if 'refGene' not in data:return f"Gene {gene_symbol} not found"gene_data = data['refGene'][0]chrom = gene_data['chrom'].replace("chr", "")start = gene_data['txStart']end = gene_data['txEnd']return f"chr{chrom}:{start}-{end}"except Exception as e:raise Exception(f"Failed to get location: {str(e)}")# 测试
print(get_gene_location_ucsc_hg19("TP53")) # 输出: chr17:7571719-7590868
方法 3:使用 NCBI API 并转换坐标(需要额外步骤)
import requests
import timedef convert_coordinates_hg38_to_hg19(chrom, start, end):"""将 hg38 坐标转换为 hg19 坐标(需要安装 CrossMap)注意:这需要本地安装 CrossMap 工具"""# 这是一个伪代码示例,实际需要安装 CrossMap 并准备链文件# 安装: pip install CrossMap(或者conda环境安装:onda install -c conda-forge biopython)# 下载链文件: wget http://hgdownload.cse.ucsc.edu/goldenpath/hg38/liftOver/hg38ToHg19.over.chain.gzimport subprocess# 创建 BED 文件bed_content = f"{chrom}\t{start}\t{end}\tfeature"# 运行 CrossMapresult = subprocess.run(["CrossMap.py", "bed", "hg38ToHg19.over.chain.gz", "-"],input=bed_content, text=True, capture_output=True)# 解析结果if result.returncode == 0:output = result.stdout.strip().split('\t')return f"chr{output[0]}:{output[1]}-{output[2]}"else:raise Exception(f"Coordinate conversion failed: {result.stderr}")def get_gene_location_ncbi_hg19(gene_symbol):"""获取基因位置并通过 CrossMap 转换为 hg19"""# 获取 hg38 位置url = f"https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=gene&term={gene_symbol}[Gene Name] AND human[Organism]&retmode=json"try:response = requests.get(url, timeout=30)response.raise_for_status()data = response.json()gene_id = data["result"]["uids"][0]gene_data = data["result"][gene_id]chrom = gene_data["chromosome"]start = gene_data["genomicinfo"][0]["chrstart"]end = gene_data["genomicinfo"][0]["chrstop"]# 转换为 hg19return convert_coordinates_hg38_to_hg19(chrom, start, end)except Exception as e:raise Exception(f"Failed to get location: {str(e)}")# 注意:此方法需要本地安装 CrossMap 和链文件
方法 4:使用本地数据库(最可靠,适合批量查询)
如果经常需要查询,建议下载 GRCh37 的基因注释文件:
1、下载 GRCh37 基因注释文件:
wget ftp://ftp.ensembl.org/pub/grch37/release-104/gtf/homo_sapiens/Homo_sapiens.GRCh37.87.gtf.gz
gunzip Homo_sapiens.GRCh37.87.gtf.gz
2、创建 Python 查询脚本:
import gzip
import os
import sqlite3
from collections import defaultdict# 创建本地SQLite数据库
def create_gene_db(gtf_path, db_path="genes_hg19.db"):if os.path.exists(db_path):return db_pathconn = sqlite3.connect(db_path)cursor = conn.cursor()# 创建表cursor.execute('''CREATE TABLE genes (gene_id TEXT,gene_name TEXT,chrom TEXT,start INTEGER,end INTEGER,strand TEXT)''')# 创建索引cursor.execute('CREATE INDEX idx_gene_name ON genes (gene_name)')# 解析GTF文件gene_data = defaultdict(lambda: {'start': float('inf'), 'end': float('-inf')})with gzip.open(gtf_path, 'rt') if gtf_path.endswith('.gz') else open(gtf_path) as f:for line in f:if line.startswith('#'):continuefields = line.strip().split('\t')if fields[2] != 'gene':continuechrom = fields[0]start = int(fields[3])end = int(fields[4])strand = fields[6]info = {k: v.strip('"') for k, v in (item.split(' ') for item in fields[8].split(';') if item)}gene_name = info.get('gene_name')gene_id = info.get('gene_id')if gene_name and gene_id:# 更新基因范围if start < gene_data[gene_name]['start']:gene_data[gene_name]['start'] = startif end > gene_data[gene_name]['end']:gene_data[gene_name]['end'] = endgene_data[gene_name].update({'chrom': chrom,'strand': strand,'gene_id': gene_id})# 插入数据库for gene_name, data in gene_data.items():cursor.execute('''INSERT INTO genes (gene_id, gene_name, chrom, start, end, strand)VALUES (?, ?, ?, ?, ?, ?)''', (data['gene_id'], gene_name, data['chrom'], data['start'], data['end'], data['strand']))conn.commit()conn.close()return db_path# 查询函数
def get_gene_location_local(gene_symbol, db_path="genes_hg19.db"):conn = sqlite3.connect(db_path)cursor = conn.cursor()cursor.execute('''SELECT chrom, start, end FROM genes WHERE gene_name = ?''', (gene_symbol,))result = cursor.fetchone()conn.close()if result:chrom, start, end = resultreturn f"chr{chrom}:{start}-{end}"else:return None# 初始化数据库(只需运行一次)
# db_path = create_gene_db("Homo_sapiens.GRCh37.87.gtf.gz")# 查询
print(get_gene_location_local("TP53")) # chr17:7571719-7590868
print(get_gene_location_local("BRCA1")) # chr17:41196312-41277500
其他
下载预处理的基因位置文件,然后解析此文件获取位置信息
wget https://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/refGene.txt.gz
gunzip refGene.txt.gz
注意事项
Ensembl API 限制:
最大每秒15个请求
每小时最多6000个请求
需要添加重试机制和延迟
基因组版本差异:
GRCh37 = hg19
GRCh38 = hg38
确保所有工具和数据库使用一致的版本
基因位置表示:
位置通常是基因的转录起始和终止位置
不同数据库可能有轻微差异(100-1000bp)
对于精确位置需求,请指定特定转录本
对于大多数应用,Ensembl REST API(方法1)是最简单直接的方法,可以快速获取 GRCh37 位置信息。