\n Generiert am: {pd.Timestamp.now().strftime('%Y-%m-%d %H:%M:%S')}\n
\n \n \n \"\"\"\n \n with open(output_path, 'w', encoding='utf-8') as f:\n f.write(html_content)\n \n print(f\"\\n✅ HTML-Report erstellt: {output_path}\")\n\ncreate_eda_html_report()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-30T06:17:52.383322Z","iopub.execute_input":"2025-11-30T06:17:52.383630Z","iopub.status.idle":"2025-11-30T06:17:52.472524Z","shell.execute_reply.started":"2025-11-30T06:17:52.383607Z","shell.execute_reply":"2025-11-30T06:17:52.471381Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"## Optimized Feature Engineering with Dimensionality Reduction\n\nclass OptimizedFeatureEngineering(FeatureEngineering):\n \"\"\"\n Optimized feature engineering with dimensionality reduction\n \"\"\"\n \n def __init__(self, kmer_sizes: List[int] = [3, 4], \n use_vj_features: bool = True,\n use_physicochemical: bool = True,\n use_diversity: bool = True,\n max_kmer_features: int = 1000,\n use_pca: bool = False,\n pca_components: int = 100):\n \"\"\"\n Initialize optimized feature engineering\n \n Args:\n kmer_sizes: List of k-mer sizes (reduced to [3,4] for speed)\n use_vj_features: Whether to include V-J gene features\n use_physicochemical: Whether to include physicochemical properties\n use_diversity: Whether to include diversity metrics\n max_kmer_features: Maximum number of k-mer features to keep\n use_pca: Whether to use PCA for dimensionality reduction\n pca_components: Number of PCA components\n \"\"\"\n super().__init__(kmer_sizes, use_vj_features, use_physicochemical, use_diversity)\n self.max_kmer_features = max_kmer_features\n self.use_pca = use_pca\n self.pca_components = pca_components\n self.pca = None\n self.top_kmers = {}\n self.top_vj = {}\n \n def select_top_features(self, features_df: pd.DataFrame, labels: pd.Series) -> pd.DataFrame:\n \"\"\"\n Select top features based on variance and correlation with labels\n \n Args:\n features_df: Features DataFrame\n labels: Target labels\n \n Returns:\n Reduced features DataFrame\n \"\"\"\n print(\"Selecting top features...\")\n \n # Separate feature types\n kmer_cols = [col for col in features_df.columns if col.startswith('kmer_')]\n vj_cols = [col for col in features_df.columns if col.startswith('vj_')]\n other_cols = [col for col in features_df.columns if not (col.startswith('kmer_') or col.startswith('vj_'))]\n \n selected_cols = other_cols.copy() # Always keep physicochemical and diversity features\n \n # Select top k-mer features based on variance\n if kmer_cols:\n kmer_df = features_df[kmer_cols]\n variances = kmer_df.var()\n top_kmer_cols = variances.nlargest(min(self.max_kmer_features, len(kmer_cols))).index.tolist()\n selected_cols.extend(top_kmer_cols)\n \n # Store for transform\n self.top_kmers = {col: True for col in top_kmer_cols}\n \n # Select top V-J features based on variance\n if vj_cols:\n vj_df = features_df[vj_cols]\n variances = vj_df.var()\n top_vj_cols = variances.nlargest(min(500, len(vj_cols))).index.tolist()\n selected_cols.extend(top_vj_cols)\n \n # Store for transform\n self.top_vj = {col: True for col in top_vj_cols}\n \n reduced_df = features_df[selected_cols]\n \n print(f\"Reduced features from {len(features_df.columns)} to {len(selected_cols)}\")\n \n return reduced_df\n \n def fit_transform(self, data_dir: str) -> Tuple[pd.DataFrame, pd.Series]:\n \"\"\"\n Fit feature engineering pipeline and transform training data\n \n Args:\n data_dir: Path to training data directory\n \n Returns:\n Tuple of (features_df, labels)\n \"\"\"\n print(\"Extracting features from training data (optimized)...\")\n \n metadata = pd.read_csv(os.path.join(data_dir, 'metadata.csv'))\n \n all_features = []\n all_labels = []\n repertoire_ids = []\n \n # Use parallel processing for feature extraction\n from joblib import Parallel, delayed\n \n def process_repertoire(row):\n \"\"\"Process single repertoire\"\"\"\n file_path = os.path.join(data_dir, row['filename'])\n repertoire_df = pd.read_csv(file_path, sep='\\t')\n features = self.extract_repertoire_features(repertoire_df)\n return features, row['label_positive'], row['repertoire_id']\n \n # Parallel processing\n results = Parallel(n_jobs=self.n_jobs if hasattr(self, 'n_jobs') else 4)(\n delayed(process_repertoire)(row) \n for _, row in tqdm(metadata.iterrows(), total=len(metadata), desc=\"Processing repertoires\")\n )\n \n # Unpack results\n for features, label, rep_id in results:\n all_features.append(features)\n all_labels.append(label)\n repertoire_ids.append(rep_id)\n \n # Convert to DataFrame\n features_df = pd.DataFrame(all_features)\n features_df.index = repertoire_ids\n features_df = features_df.fillna(0)\n \n labels = pd.Series(all_labels, index=repertoire_ids)\n \n # Feature selection\n features_df = self.select_top_features(features_df, labels)\n \n # Store feature names\n self.feature_names = features_df.columns.tolist()\n \n # Apply PCA if requested\n if self.use_pca:\n print(f\"Applying PCA (n_components={self.pca_components})...\")\n self.pca = PCA(n_components=min(self.pca_components, len(features_df.columns)))\n features_pca = self.pca.fit_transform(features_df)\n \n # Create new DataFrame with PCA components\n pca_cols = [f'pca_{i}' for i in range(features_pca.shape[1])]\n features_df = pd.DataFrame(features_pca, columns=pca_cols, index=features_df.index)\n self.feature_names = pca_cols\n \n print(f\"Explained variance ratio: {self.pca.explained_variance_ratio_.sum():.3f}\")\n else:\n # Fit scaler\n self.scaler = StandardScaler()\n features_scaled = self.scaler.fit_transform(features_df)\n features_df = pd.DataFrame(features_scaled, columns=features_df.columns, index=features_df.index)\n \n print(f\"Final feature count: {len(self.feature_names)}\")\n \n return features_df, labels\n \n def transform(self, data_dir: str) -> pd.DataFrame:\n \"\"\"\n Transform test data using fitted pipeline (optimized)\n \n Args:\n data_dir: Path to test data directory\n \n Returns:\n Features DataFrame\n \"\"\"\n print(f\"Extracting features from test data: {data_dir}\")\n \n metadata_path = os.path.join(data_dir, 'metadata.csv')\n \n all_features = []\n repertoire_ids = []\n \n # Use parallel processing\n from joblib import Parallel, delayed\n \n def process_repertoire(file_path, rep_id):\n \"\"\"Process single repertoire\"\"\"\n repertoire_df = pd.read_csv(file_path, sep='\\t')\n features = self.extract_repertoire_features(repertoire_df)\n return features, rep_id\n \n if os.path.exists(metadata_path):\n metadata = pd.read_csv(metadata_path)\n \n results = Parallel(n_jobs=self.n_jobs if hasattr(self, 'n_jobs') else 4)(\n delayed(process_repertoire)(os.path.join(data_dir, row['filename']), row['repertoire_id'])\n for _, row in tqdm(metadata.iterrows(), total=len(metadata), desc=\"Processing repertoires\")\n )\n else:\n tsv_files = glob.glob(os.path.join(data_dir, '*.tsv'))\n \n results = Parallel(n_jobs=self.n_jobs if hasattr(self, 'n_jobs') else 4)(\n delayed(process_repertoire)(file_path, os.path.basename(file_path).replace('.tsv', ''))\n for file_path in tqdm(tsv_files, desc=\"Processing repertoires\")\n )\n \n # Unpack results\n for features, rep_id in results:\n all_features.append(features)\n repertoire_ids.append(rep_id)\n \n # Convert to DataFrame\n features_df = pd.DataFrame(all_features)\n features_df.index = repertoire_ids\n \n # Select only top features\n for col in features_df.columns:\n if col.startswith('kmer_') and col not in self.top_kmers:\n features_df = features_df.drop(columns=[col])\n elif col.startswith('vj_') and col not in self.top_vj:\n features_df = features_df.drop(columns=[col])\n \n # Align with training features\n for col in self.feature_names:\n if col not in features_df.columns and not col.startswith('pca_'):\n features_df[col] = 0\n \n if not self.use_pca:\n features_df = features_df[self.feature_names]\n \n features_df = features_df.fillna(0)\n \n # Apply transformation\n if self.use_pca:\n # First scale, then PCA\n temp_scaler = StandardScaler()\n features_scaled = temp_scaler.fit_transform(features_df)\n features_pca = self.pca.transform(features_scaled)\n features_df = pd.DataFrame(features_pca, columns=self.feature_names, index=features_df.index)\n else:\n features_scaled = self.scaler.transform(features_df)\n features_df = pd.DataFrame(features_scaled, columns=features_df.columns, index=features_df.index)\n \n return features_df\n\n## Optimized Model Trainer with Early Stopping and Reduced CV\n\nclass OptimizedModelTrainer(ModelTrainer):\n \"\"\"\n Optimized model trainer with faster training\n \"\"\"\n \n def __init__(self, model_type: str = 'xgboost', n_jobs: int = 4, device: str = 'cpu',\n fast_mode: bool = True):\n \"\"\"\n Initialize optimized model trainer\n \n Args:\n model_type: Type of model\n n_jobs: Number of parallel jobs\n device: Device to use\n fast_mode: Whether to use faster hyperparameters\n \"\"\"\n super().__init__(model_type, n_jobs, device)\n self.fast_mode = fast_mode\n \n def create_model(self):\n \"\"\"Create model with optimized hyperparameters for speed\"\"\"\n \n if self.model_type == 'xgboost':\n if self.fast_mode:\n # Faster XGBoost settings\n self.model = xgb.XGBClassifier(\n n_estimators=100, # Reduced from 200\n max_depth=4, # Reduced from 6\n learning_rate=0.1, # Increased from 0.05\n subsample=0.8,\n colsample_bytree=0.8,\n random_state=42,\n n_jobs=self.n_jobs,\n tree_method='hist', # Faster than exact\n max_bin=256, # Reduced from default\n eval_metric='logloss'\n )\n else:\n # Standard settings\n self.model = xgb.XGBClassifier(\n n_estimators=200,\n max_depth=6,\n learning_rate=0.05,\n subsample=0.8,\n colsample_bytree=0.8,\n random_state=42,\n n_jobs=self.n_jobs,\n tree_method='hist',\n eval_metric='logloss'\n )\n \n elif self.model_type == 'lightgbm':\n if self.fast_mode:\n # Faster LightGBM settings\n self.model = lgb.LGBMClassifier(\n n_estimators=100,\n max_depth=4,\n learning_rate=0.1,\n subsample=0.8,\n colsample_bytree=0.8,\n random_state=42,\n n_jobs=self.n_jobs,\n num_leaves=31, # Reduced\n verbose=-1\n )\n else:\n self.model = lgb.LGBMClassifier(\n n_estimators=200,\n max_depth=6,\n learning_rate=0.05,\n subsample=0.8,\n colsample_bytree=0.8,\n random_state=42,\n n_jobs=self.n_jobs,\n verbose=-1\n )\n \n elif self.model_type == 'rf':\n if self.fast_mode:\n # Faster Random Forest\n self.model = RandomForestClassifier(\n n_estimators=50, # Reduced from 200\n max_depth=8, # Reduced from 10\n min_samples_split=10, # Increased from 5\n min_samples_leaf=4, # Increased from 2\n random_state=42,\n n_jobs=self.n_jobs,\n max_features='sqrt' # Faster than 'auto'\n )\n else:\n self.model = RandomForestClassifier(\n n_estimators=200,\n max_depth=10,\n min_samples_split=5,\n min_samples_leaf=2,\n random_state=42,\n n_jobs=self.n_jobs\n )\n \n elif self.model_type == 'lr':\n # Logistic Regression is already fast\n self.model = LogisticRegression(\n C=1.0,\n max_iter=500, # Reduced from 1000\n random_state=42,\n n_jobs=self.n_jobs,\n solver='saga' # Faster for large datasets\n )\n \n else:\n raise ValueError(f\"Unknown model type: {self.model_type}\")\n \n return self.model\n \n def cross_validate(self, X: pd.DataFrame, y: pd.Series, cv: int = 3):\n \"\"\"\n Perform cross-validation with reduced folds for speed\n \n Args:\n X: Features\n y: Labels\n cv: Number of folds (reduced to 3 for speed)\n \n Returns:\n Cross-validation scores\n \"\"\"\n print(f\"Performing {cv}-fold cross-validation...\")\n \n if self.model is None:\n self.create_model()\n \n skf = StratifiedKFold(n_splits=cv, shuffle=True, random_state=42)\n \n scores = cross_val_score(\n self.model, X, y,\n cv=skf,\n scoring='roc_auc',\n n_jobs=1, # Set to 1 to avoid nested parallelism\n verbose=0\n )\n \n print(f\"Cross-validation AUC: {scores.mean():.4f} (+/- {scores.std():.4f})\")\n \n return scores\n\n## Fast Important Sequences Identifier\n\nclass FastImportantSequencesIdentifier(ImportantSequencesIdentifier):\n \"\"\"\n Faster version of important sequences identifier\n \"\"\"\n \n def identify_frequency_based(self, data_dir: str, top_k: int = 50000, \n sample_size: int = None) -> pd.DataFrame:\n \"\"\"\n Identify important sequences with optional sampling for speed\n \n Args:\n data_dir: Path to training data directory\n top_k: Number of top sequences to return\n sample_size: Maximum number of sequences to process per repertoire (None = all)\n \n Returns:\n DataFrame with important sequences\n \"\"\"\n print(\"Identifying important sequences (frequency-based, optimized)...\")\n \n metadata = pd.read_csv(os.path.join(data_dir, 'metadata.csv'))\n \n # Collect sequences by label\n pos_sequences = Counter()\n neg_sequences = Counter()\n \n for idx, row in tqdm(metadata.iterrows(), total=len(metadata), desc=\"Processing repertoires\"):\n file_path = os.path.join(data_dir, row['filename'])\n repertoire_df = pd.read_csv(file_path, sep='\\t')\n \n sequences = repertoire_df['junction_aa'].dropna().tolist()\n \n # Sample if requested\n if sample_size and len(sequences) > sample_size:\n sequences = np.random.choice(sequences, sample_size, replace=False).tolist()\n \n if row['label_positive'] == 1:\n pos_sequences.update(sequences)\n else:\n neg_sequences.update(sequences)\n \n # Calculate enrichment scores (vectorized)\n all_sequences = set(pos_sequences.keys()) | set(neg_sequences.keys())\n \n # Pre-calculate totals\n pos_total = sum(pos_sequences.values())\n neg_total = sum(neg_sequences.values())\n \n sequence_scores = []\n \n # Batch processing\n for seq in tqdm(all_sequences, desc=\"Calculating enrichment scores\"):\n pos_count = pos_sequences[seq]\n neg_count = neg_sequences[seq]\n total = pos_count + neg_count\n \n if total >= 5: # Minimum count threshold\n pos_freq = pos_count / pos_total if pos_total > 0 else 0\n neg_freq = neg_count / neg_total if neg_total > 0 else 0\n \n if neg_freq > 0:\n log_odds = np.log2((pos_freq + 1e-10) / (neg_freq + 1e-10))\n else:\n log_odds = 10\n \n importance = abs(log_odds) * np.log(total + 1)\n \n sequence_scores.append({\n 'junction_aa': seq,\n 'pos_count': pos_count,\n 'neg_count': neg_count,\n 'log_odds_ratio': log_odds,\n 'importance_score': importance\n })\n \n # Sort and take top k\n scores_df = pd.DataFrame(sequence_scores).sort_values('importance_score', ascending=False)\n top_sequences = scores_df.head(top_k)\n \n # Add placeholder V and J genes\n top_sequences['v_call'] = 'TRBV20-1'\n top_sequences['j_call'] = 'TRBJ2-7'\n \n print(f\"Identified {len(top_sequences)} important sequences\")\n \n return top_sequences\n\n## Updated ImmuneStatePredictor with Optimizations\n\nclass OptimizedImmuneStatePredictor(ImmuneStatePredictor):\n \"\"\"\n Optimized immune state predictor with faster training\n \"\"\"\n \n def __init__(self, n_jobs: int = 1, device: str = 'cpu', \n model_type: str = 'xgboost',\n kmer_sizes: List[int] = [3, 4], # Reduced from [3,4,5]\n importance_method: str = 'frequency_based',\n fast_mode: bool = True,\n max_kmer_features: int = 1000,\n use_pca: bool = False,\n cv_folds: int = 3, # Reduced from 5\n **kwargs):\n \"\"\"\n Initialize optimized predictor\n \n Args:\n n_jobs: Number of CPU cores\n device: Device to use\n model_type: Type of model to use\n kmer_sizes: K-mer sizes for feature extraction\n importance_method: Method for identifying important sequences\n fast_mode: Whether to use faster settings\n max_kmer_features: Maximum number of k-mer features\n use_pca: Whether to use PCA\n cv_folds: Number of cross-validation folds\n \"\"\"\n total_cores = os.cpu_count()\n if n_jobs == -1:\n self.n_jobs = total_cores\n else:\n self.n_jobs = min(n_jobs, total_cores)\n \n self.device = device\n self.model_type = model_type\n self.kmer_sizes = kmer_sizes\n self.importance_method = importance_method\n self.fast_mode = fast_mode\n self.cv_folds = cv_folds\n \n # Initialize optimized components\n self.feature_eng = OptimizedFeatureEngineering(\n kmer_sizes=kmer_sizes,\n use_vj_features=True,\n use_physicochemical=True,\n use_diversity=True,\n max_kmer_features=max_kmer_features,\n use_pca=use_pca\n )\n self.feature_eng.n_jobs = self.n_jobs\n \n self.model_trainer = OptimizedModelTrainer(\n model_type=model_type,\n n_jobs=self.n_jobs,\n device=self.device,\n fast_mode=fast_mode\n )\n \n self.seq_identifier = FastImportantSequencesIdentifier(\n method=importance_method\n )\n \n self.important_sequences_ = None\n self.model = None\n \n def fit(self, train_dir_path: str):\n \"\"\"\n Train the model on training data (optimized)\n \n Args:\n train_dir_path: Path to training data directory\n \n Returns:\n self\n \"\"\"\n print(f\"\\n{'='*80}\")\n print(f\"TRAINING IMMUNE STATE PREDICTOR (OPTIMIZED)\")\n print(f\"{'='*80}\\n\")\n \n # Step 1: Feature extraction\n print(\"Step 1: Feature Extraction (Optimized)\")\n start_time = time.time()\n X_train, y_train = self.feature_eng.fit_transform(train_dir_path)\n print(f\"Feature extraction time: {time.time() - start_time:.2f}s\")\n \n # Step 2: Model training\n print(\"\\nStep 2: Model Training (Optimized)\")\n start_time = time.time()\n self.model_trainer.train(X_train, y_train)\n self.model = self.model_trainer.model\n print(f\"Training time: {time.time() - start_time:.2f}s\")\n \n # Step 3: Cross-validation (optional, can be skipped for speed)\n if self.cv_folds > 0:\n print(f\"\\nStep 3: Cross-Validation ({self.cv_folds} folds)\")\n start_time = time.time()\n cv_scores = self.model_trainer.cross_validate(X_train, y_train, cv=self.cv_folds)\n print(f\"Cross-validation time: {time.time() - start_time:.2f}s\")\n \n # Step 4: Identify important sequences\n print(\"\\nStep 4: Identifying Important Sequences (Optimized)\")\n start_time = time.time()\n self.important_sequences_ = self.seq_identifier.identify_frequency_based(\n train_dir_path,\n top_k=50000,\n sample_size=5000 if self.fast_mode else None # Sample for speed\n )\n print(f\"Sequence identification time: {time.time() - start_time:.2f}s\")\n \n print(f\"\\n{'='*80}\")\n print(\"TRAINING COMPLETE!\")\n print(f\"{'='*80}\\n\")\n \n return self\n\n## Usage Example with Timing\n\nimport time\n\nprint(\"\\n\" + \"=\"*80)\nprint(\"STARTING OPTIMIZED IMMUNE STATE PREDICTION PIPELINE\")\nprint(\"=\"*80 + \"\\n\")\n\n# Process each dataset pair\nfor train_dir, test_dirs in train_test_dataset_pairs:\n print(f\"\\n{'='*80}\")\n print(f\"Processing: {os.path.basename(train_dir)}\")\n print(f\"{'='*80}\\n\")\n \n total_start = time.time()\n \n # Initialize optimized predictor\n predictor = OptimizedImmuneStatePredictor(\n n_jobs=-1, # Use all available cores\n device='cpu',\n model_type='xgboost', # Or 'lightgbm' for even faster training\n kmer_sizes=[3, 4], # Reduced from [3,4,5]\n importance_method='frequency_based',\n fast_mode=True, # Enable fast mode\n max_kmer_features=1000, # Limit k-mer features\n use_pca=False, # Set to True for even more speed (but may reduce accuracy)\n cv_folds=3 # Reduced from 5\n )\n \n # Train\n predictor.fit(train_dir)\n \n # Save model\n model_save_dir = os.path.join(results_dir, 'models', os.path.basename(train_dir))\n predictor.save(model_save_dir)\n \n # Predict on test sets\n for test_dir in test_dirs:\n predictions = predictor.predict_proba(test_dir)\n \n # Save predictions\n pred_path = os.path.join(\n results_dir,\n f\"{os.path.basename(train_dir)}_test_predictions.tsv\"\n )\n save_tsv(predictions, pred_path)\n \n # Save important sequences\n important_seqs = predictor.identify_associated_sequences(\n dataset_name=os.path.basename(train_dir),\n top_k=50000\n )\n \n seq_path = os.path.join(\n results_dir,\n f\"{os.path.basename(train_dir)}_important_sequences.tsv\"\n )\n save_tsv(important_seqs, seq_path)\n \n total_time = time.time() - total_start\n print(f\"\\nTotal processing time for {os.path.basename(train_dir)}: {total_time:.2f}s ({total_time/60:.2f} min)\")\n\n# Concatenate all outputs\nconcatenate_output_files(out_dir=results_dir)\n\nprint(\"\\n\" + \"=\"*80)\nprint(\"OPTIMIZED PIPELINE COMPLETE!\")\nprint(\"=\"*80)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-06T10:07:07.636941Z","iopub.execute_input":"2025-11-06T10:07:07.637321Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true,"_kg_hide-input":true},"outputs":[],"execution_count":null}]}