### 使用 Python 从 GEO 数据库下载数据并进行分析 #### 准备工作 为了实现从 GEO 数据库下载数据并进行分析，可以借助 `biopython` 库中的模块来完成这一任务。Biopython 提供了一个名为 `Bio.GeoUtils` 的工具，用于解析 GEO 数据库的内容[^2]。 首先需要确保已安装 Biopython 和其他必要的依赖项。可以通过以下命令安装所需的库： ```bash pip install biopython pandas numpy requests beautifulsoup4 ``` --- #### 下载 GEO 数据 通过 Python 编写脚本可以从 GEO 数据库获取特定的 GSE/GPL/GSM 文件。以下是具体方法： 1. **使用 Entrez 工具** NCBI 提供了强大的 API 接口——Entrez，允许开发者查询和下载公共数据库中的内容。下面是一个简单的例子展示如何使用 Biopython 中的 `Entrez` 模块检索 GEO 数据。 ```python from Bio import Entrez def fetch_geo_data(gse_id, email): """Fetch data using the given GSE ID.""" Entrez.email = email # Always tell NCBI who you are handle = Entrez.esearch(db="gds", term=gse_id) record = Entrez.read(handle) gse_uid = record['IdList'][0] geo_handle = Entrez.efetch(db="gds", id=gse_uid, retmode="text") with open(f"{gse_id}_data.txt", "w") as file: file.write(geo_handle.read()) fetch_geo_data("GSE123456", "your_email@example.com") # 替换为实际的 GSE ID 和邮箱地址 ``` 2. **直接抓取矩阵文件** 如果目标是下载 `.txt` 或 `.gz` 格式的矩阵文件，则可以直接访问 GEO 官方网站上的链接，并通过 HTTP 请求将其保存到本地磁盘中。例如： ```python import os import urllib.request def download_matrix_file(url, output_dir="."): filename = url.split("/")[-1] filepath = os.path.join(output_dir, filename) urllib.request.urlretrieve(url, filepath) print(f"File downloaded to {filepath}") matrix_url = "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE123nnn/GSE123456/matrix/GSE123456_series_matrix.txt.gz" download_matrix_file(matrix_url) ``` --- #### 解析与预处理数据 一旦成功下载了原始数据文件（通常是压缩后的 TXT），下一步就是加载它以便进一步操作。这里介绍两种常见的场景及其对应的解决方案： 1. **读入 Series Matrix 文本文件** 这类文件通常包含了关于实验设计、样本属性以及表达量数值等重要信息。我们可以利用 Pandas 来简化导入过程。 ```python import gzip import pandas as pd def parse_series_matrix(file_path): with gzip.open(file_path, 'rt') as f: lines = f.readlines() header_end_index = next(i for i, line in enumerate(lines) if "!dataset_table_begin" in line.lower()) + 1 table_lines = [line.strip().split("\t") for line in lines[header_end_index:] if not line.startswith("!")] df = pd.DataFrame(table_lines[1:], columns=table_lines[0]) return df parsed_df = parse_series_matrix("GSE123456_series_matrix.txt.gz") print(parsed_df.head()) ``` 2. **基因探针映射至官方名称** 很多时候我们需要把平台定义下的唯一标识符（如 Affymetrix Probe Set IDs）转化为更通用的标准命名法（比如 HGNC 符号）。此步骤可能涉及额外资源的支持，像 MyGene.info Web Service 就非常适合此类需求。 ```python import mygene mg = mygene.MyGeneInfo() probe_ids = ["1728_at", "3969_s_at"] # Example list of probesets results = mg.querymany(probe_ids, scopes='affy_hg_u133a', fields='symbol,name') mapping_dict = {} for res in results: key = res.get('query') value = { 'Symbol': res.get('symbol'), 'Name': res.get('name') } mapping_dict[key] = value print(mapping_dict) ``` --- #### 数据可视化与统计学检验 最后一步是对整理好的数据执行探索性和验证性的计算。Matplotlib 及 Seaborn 是绘制图表的好帮手；而 SciPy 则提供了丰富的算法支持来进行差异显著性测试等工作。 示例代码如下所示： ```python import matplotlib.pyplot as plt import seaborn as sns from scipy.stats import ttest_ind def plot_expression_levels(df, gene_symbol): subset = df[df["ID_REF"].isin([k for k,v in mapping_dict.items() if v['Symbol'] == gene_symbol])] control_values = subset.loc[:,subset.columns.str.contains("Control")].values.flatten() treatment_values = subset.loc[:,~subset.columns.str.contains("Control") & ~subset.columns.str.match("ID_REF|Unnamed")].values.flatten() _, p_value = ttest_ind(control_values, treatment_values) fig, ax = plt.subplots(figsize=(8,6)) sns.boxplot(data=[control_values,treatment_values], palette=["lightblue","orange"], showmeans=True, meanprops={"marker":"o", "markerfacecolor":"white", "markersize":8, "markeredgewidth":2}) ax.set_xticklabels(['Control','Treatment']) ax.set_title(f'Expression Levels ({gene_symbol})\nP-value={p_value:.3f}') plt.show() plot_expression_levels(parsed_df, "TP53") ``` ---