バイオインフォマティクス実践ガイド

def preprocess_sequence(raw_sequence):
    # 品質管理
    quality_filtered = quality_control(raw_sequence)
    
    # アダプター除去
    trimmed = adapter_trimming(quality_filtered)
    
    # 重複除去
    deduplicated = remove_duplicates(trimmed)
    
    return deduplicated

import logging
import os
from typing import List, Dict, Optional, Tuple
from dataclasses import dataclass
from Bio import SeqIO
import pysam
import numpy as np

# ロギング設定
logging.basicConfig(
    level=logging.INFO,
    format='%(asctime)s - %(name)s - %(levelname)s - %(message)s'
)
logger = logging.getLogger(__name__)

@dataclass
class QualityMetrics:
    """品質管理メトリクス"""
    total_reads: int
    passed_reads: int
    mean_quality: float
    q20_percentage: float
    q30_percentage: float
    gc_content: float

class SequenceAnalysisPipeline:
    """エラーハンドリングを含む配列解析パイプライン"""
    
    def __init__(self, min_quality: int = 20, min_length: int = 50):
        self.min_quality = min_quality
        self.min_length = min_length
        
    def quality_control(self, fastq_file: str) -> Tuple[str, QualityMetrics]:
        """
        品質管理とフィルタリング
        
        Args:
            fastq_file: 入力FASTQファイル
            
        Returns:
            tuple: (フィルタ済みファイル名, 品質メトリクス)
        """
        try:
            output_file = fastq_file.replace('.fastq', '_filtered.fastq')
            
            total_reads = 0
            passed_reads = 0
            quality_scores = []
            gc_counts = []
            
            with open(output_file, 'w') as out_handle:
                for record in SeqIO.parse(fastq_file, 'fastq'):
                    total_reads += 1
                    
                    # 品質スコアの計算
                    quals = record.letter_annotations['phred_quality']
                    mean_qual = np.mean(quals)
                    
                    # 長さチェック
                    if len(record.seq) < self.min_length:
                        logger.debug(f"Read {record.id} too short: {len(record.seq)}")
                        continue
                    
                    # 品質チェック
                    if mean_qual < self.min_quality:
                        logger.debug(f"Read {record.id} low quality: {mean_qual:.2f}")
                        continue
                    
                    # GC含量計算
                    gc_count = record.seq.count('G') + record.seq.count('C')
                    gc_content = gc_count / len(record.seq)
                    
                    # フィルタを通過
                    passed_reads += 1
                    quality_scores.extend(quals)
                    gc_counts.append(gc_content)
                    
                    # 出力
                    SeqIO.write(record, out_handle, 'fastq')
            
            # メトリクス計算
            if passed_reads == 0:
                raise ValueError("No reads passed quality filters")
            
            quality_array = np.array(quality_scores)
            metrics = QualityMetrics(
                total_reads=total_reads,
                passed_reads=passed_reads,
                mean_quality=np.mean(quality_array),
                q20_percentage=np.sum(quality_array >= 20) / len(quality_array) * 100,
                q30_percentage=np.sum(quality_array >= 30) / len(quality_array) * 100,
                gc_content=np.mean(gc_counts)
            )
            
            logger.info(f"Quality control completed: {passed_reads}/{total_reads} reads passed")
            return output_file, metrics
            
        except FileNotFoundError:
            logger.error(f"Input file not found: {fastq_file}")
            raise
        except Exception as e:
            logger.error(f"Quality control failed: {e}")
            raise

class SomaticVariantCaller:
    """腫瘍-正常ペアからの体細胞変異検出"""
    
    def __init__(self, tumor_bam, normal_bam, reference):
        self.tumor_bam = tumor_bam
        self.normal_bam = normal_bam
        self.reference = reference
    
    def call_variants(self, min_vaf=0.05, min_depth=20):
        """
        体細胞変異の検出
        
        Args:
            min_vaf: 最小変異アレル頻度
            min_depth: 最小リード深度
        """
        # MuTect2スタイルの変異検出
        variants = []
        
        # 腫瘍特異的変異の統計的検定
        # Fisher's exact testによる有意性評価
        
        return variants

class LongReadAssembler:
    """ロングリードアセンブリ"""
    
    def overlap_layout_consensus(self, reads):
        # 1. オーバーラップグラフ構築
        overlap_graph = self.build_overlap_graph(reads)
        
        # 2. レイアウト（グラフ簡約化）
        layout = self.simplify_graph(overlap_graph)
        
        # 3. コンセンサス配列生成
        consensus = self.generate_consensus(layout)
        
        return consensus