319 lines
8.8 KiB
Markdown
Executable File
319 lines
8.8 KiB
Markdown
Executable File
# 多物种基因家族分析工作流设计
|
||
|
||
## 核心概念
|
||
|
||
### 物种角色定义
|
||
|
||
| 角色 | 说明 | 示例 |
|
||
|:---|:---|:---|
|
||
| **目标物种** (Target) | 要分析基因家族的物种 | 苦荞 (Ft) |
|
||
| **参考物种** (Reference) | 用于辅助鉴定的已知物种 | 拟南芥 (Ath)、水稻 (Os)、番茄 (Sl) |
|
||
|
||
### 文件命名规范
|
||
|
||
```
|
||
# 目标物种文件（必需）
|
||
{prefix}.pep.fasta # 目标物种蛋白序列
|
||
{prefix}.gff3 # 目标物种注释文件
|
||
{prefix}.genome.fasta # 目标物种基因组（可选）
|
||
|
||
# 参考物种文件（可选，可多个）
|
||
{ref_species}_ref.pep.fasta # 参考物种蛋白（如 Ath_ref.pep.fasta）
|
||
{ref_species}_ref.gff3 # 参考物种注释（如 Ath_ref.gff3）
|
||
|
||
# 原始数据库文件（自动识别）
|
||
Athaliana_167_TAIR10.protein_primaryTranscriptOnly.fa.gz # 拟南芥JGI格式
|
||
Osativa_323_v7.protein_primaryTranscriptOnly.fa.gz # 水稻JGI格式
|
||
```
|
||
|
||
## 工作流程
|
||
|
||
### 阶段1: 数据准备 (genome_prepare)
|
||
|
||
**输入**: 原始数据库文件（JGI/Ensembl/NCBI格式）
|
||
**处理**:
|
||
1. 识别目标物种文件（匹配prefix）
|
||
2. 识别参考物种文件（其他文件）
|
||
3. 分别处理每个物种
|
||
4. 生成标准化输出
|
||
|
||
**输出**:
|
||
```
|
||
input/{job_id}/
|
||
├── {prefix}.pep.fasta # 目标物种蛋白
|
||
├── {prefix}.gff3 # 目标物种GFF
|
||
├── {prefix}.genome.fasta # 目标物种基因组
|
||
├── Ath_ref.pep.fasta # 参考物种1蛋白
|
||
├── Ath_ref.gff3 # 参考物种1 GFF
|
||
├── Os_ref.pep.fasta # 参考物种2蛋白
|
||
└── ...
|
||
```
|
||
|
||
### 阶段2: 基因家族鉴定 (identify)
|
||
|
||
**输入**:
|
||
- 目标物种蛋白: `{prefix}.pep.fasta`
|
||
- 参考物种蛋白: `*_ref.pep.fasta`（多个）
|
||
- HMM文件或Pfam ID
|
||
|
||
**处理**:
|
||
1. HMM搜索目标物种蛋白组
|
||
2. 用参考物种蛋白构建BLAST数据库
|
||
3. 双向BLAST验证
|
||
4. 整合结果
|
||
|
||
**输出**:
|
||
```
|
||
final/{prefix}.final_members.pep.fa # 目标物种基因家族成员
|
||
final/{prefix}.identify.manifest.json # 鉴定结果清单
|
||
```
|
||
|
||
### 阶段3: 进化分析 (phylogeny)
|
||
|
||
**输入**:
|
||
- 目标物种基因家族成员
|
||
- 参考物种同源基因
|
||
|
||
**处理**:
|
||
1. 多序列比对（目标+参考物种）
|
||
2. 构建系统发育树
|
||
3. 亚家族分类
|
||
|
||
**输出**:
|
||
```
|
||
final/{prefix}.tree.nwk # 系统发育树
|
||
final/{prefix}.geneID.subfamily # 亚家族分类
|
||
```
|
||
|
||
## 关键改进点
|
||
|
||
### 1. 自动区分目标物种和参考物种
|
||
|
||
```python
|
||
class SpeciesClassifier:
|
||
"""自动分类目标物种和参考物种"""
|
||
|
||
def classify(self, files: List[Path], prefix: str) -> Tuple[SpeciesFiles, List[SpeciesFiles]]:
|
||
"""
|
||
返回: (目标物种, 参考物种列表)
|
||
"""
|
||
target = None
|
||
references = []
|
||
|
||
for file in files:
|
||
# 检查是否是目标物种
|
||
if self.is_target_species(file, prefix):
|
||
target = self.parse_target_file(file)
|
||
else:
|
||
# 参考物种
|
||
ref = self.parse_reference_file(file)
|
||
references.append(ref)
|
||
|
||
return target, references
|
||
|
||
def is_target_species(self, file: Path, prefix: str) -> bool:
|
||
"""判断文件是否属于目标物种"""
|
||
# 策略1: 文件名匹配prefix
|
||
if file.name.startswith(prefix.split('.')[0]):
|
||
return True
|
||
|
||
# 策略2: 用户明确标记
|
||
if '.target.' in file.name or '_target.' in file.name:
|
||
return True
|
||
|
||
# 策略3: 只有一个主要物种文件
|
||
return False
|
||
```
|
||
|
||
### 2. 多参考物种支持
|
||
|
||
```python
|
||
class MultiReferenceBlast:
|
||
"""多参考物种BLAST验证"""
|
||
|
||
def build_combined_db(self, reference_files: List[Path]) -> Path:
|
||
"""合并多个参考物种构建BLAST数据库"""
|
||
combined_fasta = self.work_dir / "combined_ref.pep.fa"
|
||
|
||
with open(combined_fasta, 'w') as out:
|
||
for ref_file in reference_files:
|
||
species_code = self.extract_species_code(ref_file)
|
||
with open(ref_file) as f:
|
||
for line in f:
|
||
if line.startswith('>'):
|
||
# 添加物种标识到序列ID
|
||
out.write(f">{species_code}|{line[1:]}")
|
||
else:
|
||
out.write(line)
|
||
|
||
return combined_fasta
|
||
|
||
def validate_hits(self, query_hits: List[BlastHit]) -> Dict[str, Any]:
|
||
"""验证hit来自哪个参考物种"""
|
||
results = {}
|
||
for hit in query_hits:
|
||
species = hit.target_id.split('|')[0]
|
||
if species not in results:
|
||
results[species] = []
|
||
results[species].append(hit)
|
||
|
||
return results
|
||
```
|
||
|
||
### 3. 智能物种识别
|
||
|
||
```python
|
||
class SpeciesRecognizer:
|
||
"""智能识别文件属于哪个物种"""
|
||
|
||
SPECIES_SIGNATURES = {
|
||
'Ath': {
|
||
'patterns': ['Athaliana', 'TAIR', 'AT1G', 'AT2G'],
|
||
'name': 'Arabidopsis thaliana'
|
||
},
|
||
'Os': {
|
||
'patterns': ['Osativa', 'Os01g', 'Os02g'],
|
||
'name': 'Oryza sativa'
|
||
},
|
||
'Ft': {
|
||
'patterns': ['Fagopyrum', 'FtPinG', 'tataricum'],
|
||
'name': 'Fagopyrum tataricum'
|
||
},
|
||
'Sl': {
|
||
'patterns': ['Solanum', 'Soly', 'Solyc'],
|
||
'name': 'Solanum lycopersicum'
|
||
},
|
||
'Zm': {
|
||
'patterns': ['Zea', 'Zm00001'],
|
||
'name': 'Zea mays'
|
||
}
|
||
}
|
||
|
||
def recognize_from_content(self, fasta_file: Path) -> str:
|
||
"""从FASTA文件内容识别物种"""
|
||
with open(fasta_file) as f:
|
||
# 读取前100个序列ID
|
||
ids = []
|
||
for line in f:
|
||
if line.startswith('>'):
|
||
ids.append(line[1:].split()[0])
|
||
if len(ids) >= 100:
|
||
break
|
||
|
||
# 匹配物种特征
|
||
for species_code, info in self.SPECIES_SIGNATURES.items():
|
||
match_count = sum(1 for id in ids if any(p in id for p in info['patterns']))
|
||
if match_count > len(ids) * 0.5: # 超过50%匹配
|
||
return species_code
|
||
|
||
return "UNKNOWN"
|
||
```
|
||
|
||
## 实现建议
|
||
|
||
### 方案A: 文件名约定（简单）
|
||
|
||
用户按约定命名文件：
|
||
```
|
||
# 目标物种
|
||
Ft.pep.fasta # 必需
|
||
Ft.gff3 # 可选
|
||
|
||
# 参考物种
|
||
Ath_ref.pep.fasta # 拟南芥参考
|
||
Os_ref.pep.fasta # 水稻参考
|
||
Sl_ref.pep.fasta # 番茄参考
|
||
```
|
||
|
||
### 方案B: 配置文件（灵活）
|
||
|
||
```yaml
|
||
# species_config.yaml
|
||
target:
|
||
code: "Ft"
|
||
name: "Fagopyrum tataricum"
|
||
files:
|
||
pep: "Ft.pep.fasta"
|
||
gff: "Ft.gff3"
|
||
|
||
references:
|
||
- code: "Ath"
|
||
name: "Arabidopsis thaliana"
|
||
files:
|
||
pep: "Athaliana_167_TAIR10.protein_primaryTranscriptOnly.fa.gz"
|
||
gff: "Athaliana_167_TAIR10.gene_exons.gff3.gz"
|
||
|
||
- code: "Os"
|
||
name: "Oryza sativa"
|
||
files:
|
||
pep: "Osativa_323_v7.protein_primaryTranscriptOnly.fa.gz"
|
||
|
||
- code: "Sl"
|
||
name: "Solanum lycopersicum"
|
||
files:
|
||
pep: "Slycopersicum_...fa.gz"
|
||
```
|
||
|
||
### 方案C: 自动识别（智能）
|
||
|
||
系统自动识别：
|
||
1. 扫描所有文件
|
||
2. 从内容识别物种
|
||
3. 用户确认或修正
|
||
|
||
## 推荐实现
|
||
|
||
**短期**（本周）：方案A - 文件名约定
|
||
- 简单，用户容易理解
|
||
- 改动最小
|
||
|
||
**中期**（本月）：方案B - 配置文件
|
||
- 更灵活，支持复杂场景
|
||
- 可扩展性好
|
||
|
||
**长期**（未来）：方案C - 自动识别+LLM
|
||
- 完全自动化
|
||
- 智能推荐参考物种
|
||
|
||
## 立即实施的修改
|
||
|
||
### 1. 修改 `genome_prepare`
|
||
|
||
```python
|
||
def main():
|
||
parser.add_argument("--target-prefix", required=True, help="目标物种前缀")
|
||
parser.add_argument("--reference-species", nargs='+', help="参考物种代码列表")
|
||
|
||
# 处理逻辑
|
||
target_files = find_target_files(input_dir, args.target_prefix)
|
||
reference_files = find_reference_files(input_dir, exclude=target_files)
|
||
|
||
# 分别处理
|
||
process_target(target_files)
|
||
for ref in reference_files:
|
||
process_reference(ref)
|
||
```
|
||
|
||
### 2. 修改 `identify`
|
||
|
||
```python
|
||
def run_identify():
|
||
# 加载目标物种
|
||
target_pep = load_target_protein(prefix)
|
||
|
||
# 加载参考物种
|
||
reference_peps = load_reference_proteins(workspace)
|
||
|
||
# HMM搜索目标物种
|
||
hmm_results = run_hmmsearch(target_pep)
|
||
|
||
# BLAST验证（使用所有参考物种）
|
||
combined_db = build_combined_blast_db(reference_peps)
|
||
blast_results = run_blastp(target_pep, combined_db)
|
||
|
||
# 整合结果
|
||
final_members = integrate_results(hmm_results, blast_results)
|
||
```
|
||
|
||
这样设计是否满足你的需求？需要我立即实施哪个方案？
|