# 多物种基因家族分析工作流设计

## 核心概念

### 物种角色定义

| 角色 | 说明 | 示例 |
|:---|:---|:---|
| **目标物种** (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)
```

这样设计是否满足你的需求？需要我立即实施哪个方案？