PlantAI-dev/docs/multi_species_workflow_design.md

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

核心概念

物种角色定义

角色	说明	示例
目标物种 (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. 自动区分目标物种和参考物种

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. 多参考物种支持

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. 智能物种识别

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: 配置文件（灵活）

# 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

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

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)

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