Files
usher-exploring/README.md
gbanyan dc36730cb4 docs: update gitignore to track report results, add README cross-links
- Revised .gitignore to ignore raw data/cache but track data/report/
  (candidates TSV/Parquet, plots, reproducibility metadata)
- Added zh↔en cross-links between README.md and README.en.md

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
2026-02-13 05:57:27 +08:00

549 lines
22 KiB
Markdown
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
# Usher Cilia Candidate Gene Discovery Pipeline
> **[English version (README.en.md)](README.en.md)**
一套可重現的生物資訊分析管線,用於系統性篩選與 Usher 症候群及纖毛病變ciliopathies相關的候選基因。
本管線對人類約 22,600 個蛋白質編碼基因,透過六個獨立的證據層面進行評分與排序,最終產出分層候選基因清單,供後續實驗驗證參考。
---
## 目錄
- [研究背景](#研究背景)
- [管線總覽](#管線總覽)
- [安裝與環境設定](#安裝與環境設定)
- [執行管線](#執行管線)
- [六大證據層面詳解](#六大證據層面詳解)
- [1. gnomAD 約束性分析](#1-gnomad-約束性分析)
- [2. 基因功能註釋](#2-基因功能註釋)
- [3. 組織表達特異性](#3-組織表達特異性)
- [4. 亞細胞定位](#4-亞細胞定位)
- [5. 動物模型表型](#5-動物模型表型)
- [6. 文獻探勘](#6-文獻探勘)
- [綜合評分與分層](#綜合評分與分層)
- [驗證機制](#驗證機制)
- [產出檔案說明](#產出檔案說明)
- [設定檔調整](#設定檔調整)
- [已知限制](#已知限制)
- [引用資料庫](#引用資料庫)
---
## 研究背景
**Usher 症候群**是最常見的遺傳性聾盲症候群,患者同時出現感音神經性聽力損失與視網膜色素變性。目前已知的致病基因(如 MYO7A、USH2A、CDH23 等)約 10 個,但臨床上仍有部分患者無法找到已知基因的致病變異,暗示可能存在尚未被發現的致病或修飾基因。
Usher 蛋白在纖毛cilia與纖毛相關結構中扮演關鍵角色特別是視網膜光受器細胞的連接纖毛connecting cilium及內耳毛細胞的靜纖毛stereocilia。因此本管線以**纖毛生物學**為核心,整合多面向的基因體與功能體學資料,系統性地搜尋可能被忽略的候選基因。
### 核心設計理念
- **缺失資料 ≠ 零分**若某基因在特定證據層無資料NULL不會被扣分僅以有資料的層面計算加權平均。這避免了對研究不足基因的系統性懲罰。
- **多面向正交驗證**:六個證據層面分別從不同角度衡量基因與纖毛/Usher 的關聯性,降低單一資料來源偏差的影響。
- **可重現性**:所有資料版本、參數設定、分析步驟均有完整記錄。
---
## 管線總覽
```
┌─────────────────────────────────────────────────────┐
│ Step 1: 建立基因宇宙 (Gene Universe) │
│ 透過 mygene API 取得 ~22,600 個人類蛋白質編碼基因 │
└────────────────────────┬────────────────────────────┘
┌─────────────────────────────────────────────────────┐
│ Step 2: 六大證據層面 (並行可獨立執行) │
│ │
│ ┌───────────┐ ┌───────────┐ ┌────────────────┐ │
│ │ gnomAD │ │ 功能註釋 │ │ 組織表達特異性 │ │
│ │ 約束性 │ │ │ │ │ │
│ └───────────┘ └───────────┘ └────────────────┘ │
│ ┌───────────┐ ┌───────────┐ ┌────────────────┐ │
│ │ 亞細胞定位 │ │ 動物模型 │ │ 文獻探勘 │ │
│ │ │ │ 表型 │ │ │ │
│ └───────────┘ └───────────┘ └────────────────┘ │
└────────────────────────┬────────────────────────────┘
┌─────────────────────────────────────────────────────┐
│ Step 3: 綜合加權評分 + 信心分層 │
│ NULL-aware weighted average → HIGH/MEDIUM/LOW │
└────────────────────────┬────────────────────────────┘
┌─────────────────────────────────────────────────────┐
│ Step 4: 驗證 + 報告產出 │
│ 已知 Usher/cilia 基因排名驗證 → TSV + Parquet + 圖表 │
└─────────────────────────────────────────────────────┘
```
---
## 安裝與環境設定
### 系統需求
- Python 3.11 以上
- 約 5 GB 磁碟空間(用於快取下載的資料庫檔案)
- 網路連線(首次執行需下載外部資料)
### 安裝步驟
打開終端機Terminal進入專案資料夾後執行
```bash
# 1. 建立虛擬環境(建議)
python -m venv .venv
source .venv/bin/activate # macOS / Linux
# .venv\Scripts\activate # Windows
# 2. 安裝管線
pip install -e ".[dev]"
```
安裝完成後,可以用以下指令確認是否成功:
```bash
usher-pipeline info
```
如果看到版本號與設定摘要,表示安裝成功。
### NCBI API Key文獻探勘層需要
文獻探勘層使用 NCBI PubMed E-utilities API。此 API 需要提供電子郵件,並建議申請 API Key 以提高查詢速度3 次/秒 → 10 次/秒)。
申請方式:前往 https://www.ncbi.nlm.nih.gov/account/ 註冊 NCBI 帳號,在帳號設定中產生 API Key。
---
## 執行管線
以下是完整的執行流程。每一步都是獨立的指令,按順序執行即可。
> **提示**:每個指令執行完畢後會顯示摘要統計。如果中途失敗,修正問題後重新執行同一步即可(所有步驟皆為冪等操作,重複執行不會產生重複資料)。
### Step 1建立基因宇宙
```bash
usher-pipeline setup
```
這一步會:
- 透過 mygene API 查詢所有人類蛋白質編碼基因(約 22,600 個)
- 建立 Ensembl Gene ID ↔ HGNC Symbol ↔ UniProt Accession 的對應關係
- 結果儲存至本地 DuckDB 資料庫(`data/pipeline.duckdb`
### Step 2執行六大證據層面
六個層面可以按任意順序執行,彼此之間互相獨立。但請**逐一執行**不要同時開兩個DuckDB 僅支援單一寫入者)。
```bash
# 2a. gnomAD 約束性指標
usher-pipeline evidence gnomad
# 2b. 基因功能註釋
usher-pipeline evidence annotation
# 2c. 組織表達特異性
usher-pipeline evidence expression
# 2d. 亞細胞定位
usher-pipeline evidence localization
# 2e. 動物模型表型
usher-pipeline evidence animal-models
# 2f. 文獻探勘(需要 email建議提供 API key
usher-pipeline evidence literature --email your@email.com --api-key YOUR_KEY
```
> **注意**:文獻探勘層需要查詢 22,600 個基因的 PubMed 記錄,速度較慢(約 8 基因/分鐘完整執行可能需要較長時間。此層支援中斷後續跑checkpoint-restart如果中途斷線重新執行同一指令即可從上次進度繼續。
### Step 3綜合評分
```bash
usher-pipeline score
```
此步驟會:
- 以 LEFT JOIN 合併六個證據層的分數
- 計算 NULL-aware 加權綜合分數
- 對已知 Usher/纖毛基因進行排名驗證(正控制)
- 將基因分為 HIGH / MEDIUM / LOW 三個信心層級
### Step 4產出報告
```bash
usher-pipeline report
```
產出檔案會放在 `data/report/` 目錄,包含:
- `candidates.tsv` — 候選基因清單(可用 Excel 開啟)
- `candidates.parquet` — 同上,但為高效能二進位格式(適合 R / Python 後續分析)
- `score_distribution.png` — 分數分佈圖
- `evidence_coverage.png` — 各證據層覆蓋率圖
- `tier_distribution.png` — 信心層級分佈圓餅圖
- `reproducibility.md` — 可重現性報告(含所有參數與資料版本)
### 選用:驗證報告
```bash
usher-pipeline validate
```
驗證已知 Usher 基因與 SYSCILIA 纖毛基因是否排名在前 25%,確認評分系統的有效性。
---
## 六大證據層面詳解
### 1. gnomAD 約束性分析
**生物學問題**此基因是否在人類族群中顯示強烈的功能喪失Loss-of-Function選擇約束
**資料來源**gnomAD v4.1 約束性指標
**科學依據**
如果一個基因對正常生理功能是必要的例如感覺功能那麼在人類族群中應該很少觀察到該基因的功能喪失型變異LoF variants因為帶有這些變異的個體會有降低的適應性。gnomAD 資料庫透過比較「觀察到的 LoF 變異數量」與「預期的 LoF 變異數量」,計算出約束性指標。
**關鍵指標**
- **LOEUF**Loss-of-function Observed/Expected Upper bound Fraction越低代表約束越強
- **pLI**Probability of LoF Intolerance越高代表越不容忍 LoF 變異
**評分方式**
將 LOEUF 值反轉並正規化至 01 範圍:
```
loeuf_normalized = (LOEUF_max LOEUF) / (LOEUF_max LOEUF_min)
```
分數越高 = 約束性越強 = 該基因對正常功能越重要。
**品質控管**:要求平均定序深度 ≥30x 且 CDS 涵蓋率 ≥90%,低於此閾值的基因標記為 `incomplete_coverage`
---
### 2. 基因功能註釋
**生物學問題**:此基因的功能註釋完整程度如何?
**資料來源**
- Gene OntologyGO— 生物過程、分子功能、細胞組分
- UniProt — 蛋白質註釋品質分數05
- KEGG / Reactome — 代謝與訊號傳遞路徑
**科學依據**
功能註釋的完整程度反映了一個基因被研究的深度。注意:此層面與「新穎性」呈**負相關**——註釋越少的基因可能代表尚未被發現的生物學。因此在綜合評分中,此層僅佔 15% 權重,讓缺乏註釋但有其他證據支持的新穎候選基因仍能獲得高排名。
**評分方式**
```
annotation_score = 0.5 × GO 組分 + 0.3 × UniProt 組分 + 0.2 × Pathway 組分
其中:
GO 組分 = log₂(GO term 數量 + 1) / log₂(最大 GO term 數量 + 1)
UniProt 組分 = UniProt 分數 / 5.0
Pathway 組分 = 若有任何 pathway 註釋則為 1否則為 0
```
---
### 3. 組織表達特異性
**生物學問題**:此基因是否特異性地表達在 Usher 症候群相關的組織(視網膜、內耳)?
**資料來源**
- **HPA**Human Protein Atlasv23 — 組織層級 RNA 表達量TPM
- **GTEx** v8 — 54 個人體組織的大量 RNA 定序資料
- **CellxGene** — 單細胞 RNA 定序(光受器細胞、毛細胞群體)
**科學依據**
Usher 症候群影響視網膜光受器細胞與耳蝸毛細胞。若一個基因在這些組織中高度富集表達,暗示其在這些組織具有特化功能,是值得關注的候選基因。
**關鍵指標**
1. **Tau 組織特異性指數**(τ):衡量一個基因在各組織間的表達均勻程度
- τ = 0遍在表達housekeeping gene
- τ = 1高度組織特異性
- 公式:`τ = Σ(1 xᵢ/x_max) / (n 1)`
2. **Usher 組織富集度**:目標組織(視網膜、小腦、光受器、毛細胞)的平均表達量,除以所有組織的整體平均表達量。比值 > 1 代表在目標組織富集。
**評分方式**
```
expression_score = 0.4 × 富集度百分位數 + 0.3 × τ 特異性 + 0.3 × 目標組織最大表達百分位數
```
使用百分位數正規化,避免絕對 TPM 值因定序深度差異而產生偏差。
---
### 4. 亞細胞定位
**生物學問題**:此蛋白質是否定位在纖毛、中心體或基體等與 Usher 病理機轉相關的亞細胞結構?
**資料來源**
- **HPA 亞細胞定位資料** — 免疫螢光顯微鏡
- **纖毛蛋白質體學** — 已發表的纖毛質譜資料集CiliaCarta 等)
- **中心體蛋白質體學** — 已發表的中心體質譜資料集
**科學依據**
已知的 Usher 蛋白MYO7A、USH1C、CDH23 等)都是纖毛或纖毛周邊蛋白。蛋白質定位在纖毛/基體/中心體,是參與纖毛病變途徑的**強力機制性證據**。
**評分方式**
| 定位區域 | 基礎分數 |
|---------|---------|
| 纖毛、中心體、基體、轉換區、靜纖毛 | 1.0 |
| 細胞骨架、微管、細胞連接 | 0.5 |
| 僅出現在蛋白質體學資料集 | 0.3 |
| 有定位資料但與纖毛無關 | 0.0 |
| 無定位資料 | NULL |
**證據權重**
- 實驗證據HPA Enhanced/Supported 可靠度、蛋白質體學):× 1.0
- 計算預測HPA Approved/Uncertain 可靠度):× 0.6
```
localization_score = 纖毛接近度基礎分數 × 證據類型權重
```
---
### 5. 動物模型表型
**生物學問題**在小鼠或斑馬魚中此基因的直系同源物ortholog被破壞時是否會產生感覺或纖毛相關的表型
**資料來源**
- **HCOP**HUGO Gene Nomenclature Committee Ortholog Predictions— 人-鼠/人-魚直系同源物對應
- **MGI**Mouse Genome Informatics— 小鼠基因剔除表型MP ontology
- **ZFIN**Zebrafish Information Network— 斑馬魚突變表型
- **IMPC**International Mouse Phenotyping Consortium— 系統性基因剔除篩選
**科學依據**
跨物種的感覺表型保守性提供了**功能驗證**。小鼠的聽力/平衡缺陷及斑馬魚的側線系統缺陷,都能重現人類 Usher 症候群的病理特徵。
**表型關鍵字篩選**
- 小鼠hearing、vision、retina、photoreceptor、cochlea、stereocilia、cilia、vestibular、balance
- 斑馬魚hearing、ear、otic、otolith、lateral line、hair cell、retina、vision、eye
**評分方式**
```
animal_model_score =
0.4 × 小鼠分數 × 小鼠信心度 +
0.3 × 斑馬魚分數 × 斑馬魚信心度 +
0.3 × IMPC 額外加分
信心度權重HIGH = 1.0, MEDIUM = 0.7, LOW = 0.4
表型數量縮放log₂(感覺表型數量 + 1) / log₂(最大數量 + 1)
```
使用對數縮放防止註釋數量過多的基因主導排名(報酬遞減效應)。
---
### 6. 文獻探勘
**生物學問題**:此基因在科學文獻中是否被提及與纖毛或感覺功能相關?文獻證據的品質如何?
**資料來源**NCBI PubMed透過 E-utilities API
**科學依據**
文獻證據反映了累積的科學知識,但存在**研究偏差**study bias——明星基因如 TP53、BRCA1擁有大量文獻但這不代表它們與纖毛相關。本層的正規化策略特別針對這個問題進行校正。
**證據品質分級**
| 等級 | 說明 | 權重 |
|------|------|------|
| direct_experimental | 基因剔除/突變 + 纖毛/感覺語境 | 1.0 |
| functional_mention | 纖毛/感覺語境 + ≥3 篇文獻 | 0.6 |
| hts_hit | 高通量篩選命中 + 纖毛/感覺語境 | 0.3 |
| incidental | 有文獻但無纖毛/感覺語境 | 0.1 |
| none | 無文獻 | 0.0 |
**研究偏差校正**
```
raw_score = (語境分數 × 品質權重) / log₂(PubMed 總文獻數 + 1)
literature_score = percentile_rank(raw_score)
```
除以 `log₂(總文獻數)` 的關鍵意義:一個基因若在 50 篇文獻中有 5 篇提及纖毛,會比在 100,000 篇文獻中有 5 篇提及纖毛得到更高的分數。這鼓勵發現被忽略但有線索的基因。
---
## 綜合評分與分層
### 加權綜合分數
六個證據層以可調整的權重合併為單一綜合分數:
| 證據層 | 預設權重 | 生物學意義 |
|--------|---------|-----------|
| gnomAD 約束性 | 20% | 基因的功能必要性 |
| 組織表達 | 20% | 在目標組織的特異性表達 |
| 功能註釋 | 15% | 已知的功能特徵 |
| 亞細胞定位 | 15% | 纖毛結構的接近程度 |
| 動物模型 | 15% | 跨物種的功能驗證 |
| 文獻探勘 | 15% | 文獻中的纖毛/感覺關聯 |
**NULL-aware 加權平均**
```
composite_score = Σ(scoreᵢ × weightᵢ) / Σ(weightᵢ)
─── 僅計算非 NULL 的層面 ───
```
例如:某基因只有 gnomAD0.8、表達0.6和定位0.9)三層有資料:
```
composite_score = (0.8×0.20 + 0.6×0.20 + 0.9×0.15) / (0.20 + 0.20 + 0.15)
= 0.415 / 0.55
= 0.755
```
### 信心分層
| 層級 | 條件 | 意義 |
|------|------|------|
| **HIGH** | 綜合分數 ≥ 0.7 且 ≥ 3 層有資料 | 高優先候選基因,建議進入實驗驗證 |
| **MEDIUM** | 綜合分數 ≥ 0.4 且 ≥ 2 層有資料 | 中等證據,值得進一步文獻調查 |
| **LOW** | 綜合分數 ≥ 0.2 | 證據薄弱,需要更多資料 |
| EXCLUDED | 低於以上條件 | 排除,不列入候選清單 |
### 品質旗標
| 旗標 | 條件 | 說明 |
|------|------|------|
| sufficient_evidence | ≥ 4 層有分數 | 資料涵蓋充分 |
| moderate_evidence | ≥ 2 層有分數 | 部分涵蓋 |
| sparse_evidence | ≥ 1 層有分數 | 資料稀疏 |
| no_evidence | 0 層有分數 | 完全無資料 |
---
## 驗證機制
管線內建正控制驗證,確認評分系統的有效性:
### 正控制基因集
1. **OMIM Usher 基因**10 個MYO7A、USH1C、CDH23、PCDH15、USH1G、CIB2、USH2A、ADGRV1、WHRN、CLRN1
2. **SYSCILIA SCGS v2 核心纖毛基因**28 個IFT88、IFT140、BBS1、CEP290、RPGR 等
### 驗證標準
- 已知基因的中位百分位排名應 ≥ 75%(前四分之一)
- 前 10% 候選基因應包含 > 70% 的已知基因Recall@10%
如果驗證未通過,表示權重設定或資料品質可能有問題,需要檢視並調整。
---
## 產出檔案說明
所有結果儲存在 `data/report/` 目錄:
| 檔案 | 說明 | 適用對象 |
|------|------|---------|
| `candidates.tsv` | 候選基因清單Tab 分隔) | 用 Excel 或文字編輯器開啟 |
| `candidates.parquet` | 同上,高效能二進位格式 | R`arrow::read_parquet`)或 Python`polars.read_parquet` |
| `score_distribution.png` | 分數分佈直方圖 | 快速了解整體分佈 |
| `evidence_coverage.png` | 各證據層覆蓋率圖 | 了解哪些層資料最完整 |
| `tier_distribution.png` | HIGH/MEDIUM/LOW 比例圓餅圖 | 快速概覽 |
| `reproducibility.md` | 完整的參數、資料版本與軟體環境記錄 | 論文方法學 / 重現分析 |
### candidates.tsv 欄位說明
**核心欄位**
- `gene_id` — Ensembl 基因 ID如 ENSG00000154229
- `gene_symbol` — HGNC 基因符號(如 MYO7A
- `composite_score` — 綜合分數01
- `confidence_tier` — 信心層級HIGH / MEDIUM / LOW
- `evidence_count` — 有非 NULL 分數的證據層數量06
**各層分數**01NULL 表示無資料):
- `gnomad_score` — 約束性分數
- `expression_score` — 表達特異性分數
- `annotation_score` — 功能註釋分數
- `localization_score` — 亞細胞定位分數
- `animal_model_score` — 動物模型表型分數
- `literature_score` — 文獻證據分數
**輔助欄位**
- `supporting_layers` — 有分數的層面清單(如 `gnomad,expression,localization`
- `evidence_gaps` — 缺少資料的層面清單
- `quality_flag` — 品質旗標
---
## 設定檔調整
管線設定檔位於 `config/default.yaml`
```yaml
# 資料版本
versions:
ensembl_release: 113
gnomad_version: v4.1
gtex_version: v8
hpa_version: "23.0"
# 評分權重(可依研究需求調整,總和必須為 1.0
scoring:
gnomad: 0.20
expression: 0.20
annotation: 0.15
localization: 0.15
animal_model: 0.15
literature: 0.15
# API 設定
api:
rate_limit_per_second: 5
max_retries: 5
timeout_seconds: 30
```
若要調整權重(例如更重視組織表達),修改 `scoring` 區塊後重新執行 `usher-pipeline score``usher-pipeline report` 即可。
---
## 已知限制
| 限制 | 說明 | 影響 |
|------|------|------|
| GTEx v8 無視網膜組織 | GTEx v8 不包含 "Eye - Retina" 組織 | 視網膜表達資料僅來自 HPA |
| HPA 基因符號比對落差 | HPA 使用 gene symbol管線以 gene ID 為主鍵 | 部分基因的 HPA 表達資料可能遺失 |
| gnomAD 轉錄本層級 ID | gnomAD 使用轉錄本 ID非基因層級 | 部分基因 JOIN 時可能產生 NaN |
| 文獻探勘速度 | NCBI API 限制,~8 基因/分鐘 | 完整執行需要較長時間;可使用 API key 加速 |
| 單一寫入限制 | DuckDB 不支援多行程同時寫入 | 請勿同時執行兩個管線步驟 |
| 內耳資料稀缺 | 人類耳蝸組織的大量轉錄體資料有限 | 以小腦(含纖毛)及 CellxGene 毛細胞資料作為替代 |
---
## 引用資料庫
本管線整合以下公開資料庫與工具:
- **gnomAD** v4.1 — Karczewski et al. (2020) *Nature* 581:434-443
- **Human Protein Atlas** v23 — Uhlén et al. (2015) *Science* 347:1260419
- **GTEx** v8 — GTEx Consortium (2020) *Science* 369:1318-1330
- **CellxGene** — Chan Zuckerberg Initiative single-cell atlas
- **Gene Ontology** — Gene Ontology Consortium (2021) *Nucleic Acids Res* 49:D325-D334
- **UniProt** — UniProt Consortium (2023) *Nucleic Acids Res* 51:D523-D531
- **MGI** — Mouse Genome Informatics, The Jackson Laboratory
- **ZFIN** — Zebrafish Information Network
- **IMPC** — International Mouse Phenotyping Consortium
- **SYSCILIA SCGS v2** — van Dam et al. (2021) *Mol Biol Cell* 32:br6
- **mygene.info** — Xin et al. (2016) *Genome Biol* 17:91
- **NCBI E-utilities** — PubMed literature search API
---
## 授權
MIT License