1073 lines
44 KiB
Python
1073 lines
44 KiB
Python
|
#!/usr/bin/env python3
|
|||
|
|
"""
|
|||
|
|
RNA-Seq QC阶段：对原始测序数据进行质量控制
|
|||
|
|
|
|||
|
|
输入：
|
|||
|
|
- 原始FASTQ文件（双端）
|
|||
|
|
- 样本配置文件（group_sample.yaml）
|
|||
|
|
|
|||
|
|
输出：
|
|||
|
|
- FastQC报告
|
|||
|
|
- Cutadapt处理后的clean reads
|
|||
|
|
- 质控报告和manifest
|
|||
|
|
"""
|
|||
|
|
|
|||
|
|
import argparse
|
|||
|
|
import json
|
|||
|
|
import subprocess
|
|||
|
|
import sys
|
|||
|
|
import time
|
|||
|
|
import yaml
|
|||
|
|
import shutil
|
|||
|
|
from concurrent.futures import ThreadPoolExecutor, as_completed
|
|||
|
|
from datetime import datetime
|
|||
|
|
from pathlib import Path
|
|||
|
|
from typing import Any, Dict, List, Optional, Sequence
|
|||
|
|
import hashlib
|
|||
|
|
import os
|
|||
|
|
|
|||
|
|
PROJECT_ROOT = Path(__file__).resolve().parents[2]
|
|||
|
|
if str(PROJECT_ROOT) not in sys.path:
|
|||
|
|
sys.path.insert(0, str(PROJECT_ROOT))
|
|||
|
|
|
|||
|
|
from gfa_core.rnaseq_context import apply_rnaseq_context, build_rnaseq_context, resolve_rnaseq_project_context
|
|||
|
|
from apps.cli.rnaseq_env import prepend_path, rnaseq_bin
|
|||
|
|
|
|||
|
|
|
|||
|
|
RNASEQ_BIN = str(rnaseq_bin())
|
|||
|
|
prepend_path(RNASEQ_BIN)
|
|||
|
|
|
|||
|
|
def validate_input_files(input_dir: Path, config: Dict[str, Any]) -> Dict[str, Dict[str, Path]]:
|
|||
|
|
"""
|
|||
|
|
严格验证输入文件：检查R1+R2配对、文件存在、文件大小>0
|
|||
|
|
|
|||
|
|
Returns:
|
|||
|
|
验证通过的样本字典 {sample_name: {"r1": Path, "r2": Path}}
|
|||
|
|
|
|||
|
|
Raises:
|
|||
|
|
SystemExit: 验证失败时退出并返回错误信息
|
|||
|
|
"""
|
|||
|
|
samples = config.get("samples", {})
|
|||
|
|
if not samples:
|
|||
|
|
print("❌ ERROR: No samples defined in configuration file", file=sys.stderr)
|
|||
|
|
sys.exit(1)
|
|||
|
|
|
|||
|
|
validated_samples = {}
|
|||
|
|
errors = []
|
|||
|
|
|
|||
|
|
for sample_name, paths in samples.items():
|
|||
|
|
# 解析路径
|
|||
|
|
if isinstance(paths, list):
|
|||
|
|
if len(paths) < 2:
|
|||
|
|
errors.append(f"Sample '{sample_name}': insufficient paths (need 2, got {len(paths)})")
|
|||
|
|
continue
|
|||
|
|
r1_path, r2_path = paths[0], paths[1]
|
|||
|
|
elif isinstance(paths, dict):
|
|||
|
|
r1_path = paths.get("r1") or paths.get("read1")
|
|||
|
|
r2_path = paths.get("r2") or paths.get("read2")
|
|||
|
|
if not r1_path or not r2_path:
|
|||
|
|
errors.append(f"Sample '{sample_name}': missing r1/r2 paths")
|
|||
|
|
continue
|
|||
|
|
else:
|
|||
|
|
errors.append(f"Sample '{sample_name}': unknown format {type(paths)}")
|
|||
|
|
continue
|
|||
|
|
|
|||
|
|
# 解析完整路径
|
|||
|
|
r1 = Path(r1_path) if Path(r1_path).is_absolute() else input_dir / r1_path
|
|||
|
|
r2 = Path(r2_path) if Path(r2_path).is_absolute() else input_dir / r2_path
|
|||
|
|
|
|||
|
|
# 验证文件存在
|
|||
|
|
if not r1.exists():
|
|||
|
|
errors.append(f"Sample '{sample_name}': R1 file not found: {r1}")
|
|||
|
|
continue
|
|||
|
|
if not r2.exists():
|
|||
|
|
errors.append(f"Sample '{sample_name}': R2 file not found: {r2}")
|
|||
|
|
continue
|
|||
|
|
|
|||
|
|
# 验证文件大小>0
|
|||
|
|
r1_size = r1.stat().st_size
|
|||
|
|
r2_size = r2.stat().st_size
|
|||
|
|
if r1_size == 0:
|
|||
|
|
errors.append(f"Sample '{sample_name}': R1 file is empty: {r1}")
|
|||
|
|
continue
|
|||
|
|
if r2_size == 0:
|
|||
|
|
errors.append(f"Sample '{sample_name}': R2 file is empty: {r2}")
|
|||
|
|
continue
|
|||
|
|
|
|||
|
|
# 验证R1和R2大小合理（差距不超过10倍）
|
|||
|
|
if r1_size > 0 and r2_size > 0:
|
|||
|
|
ratio = max(r1_size, r2_size) / min(r1_size, r2_size)
|
|||
|
|
if ratio > 10:
|
|||
|
|
errors.append(f"Sample '{sample_name}': R1/R2 file size mismatch (ratio: {ratio:.1f})")
|
|||
|
|
continue
|
|||
|
|
|
|||
|
|
validated_samples[sample_name] = {"r1": r1, "r2": r2, "r1_size": r1_size, "r2_size": r2_size}
|
|||
|
|
|
|||
|
|
# 报告错误并退出
|
|||
|
|
if errors:
|
|||
|
|
print(f"\n❌ ERROR: Input validation failed for {len(errors)} sample(s):\n", file=sys.stderr)
|
|||
|
|
for err in errors:
|
|||
|
|
print(f" - {err}", file=sys.stderr)
|
|||
|
|
print(f"\nTotal samples in config: {len(samples)}", file=sys.stderr)
|
|||
|
|
print(f"Valid samples: {len(validated_samples)}", file=sys.stderr)
|
|||
|
|
sys.exit(1)
|
|||
|
|
|
|||
|
|
print(f"✅ Input validation passed: {len(validated_samples)} sample(s) ready")
|
|||
|
|
for name, info in validated_samples.items():
|
|||
|
|
print(f" - {name}: R1={info['r1_size']/1024/1024:.1f}MB, R2={info['r2_size']/1024/1024:.1f}MB")
|
|||
|
|
|
|||
|
|
return validated_samples
|
|||
|
|
|
|||
|
|
|
|||
|
|
def compute_file_hash(filepath: Path) -> str:
|
|||
|
|
"""计算文件的SHA256哈希值"""
|
|||
|
|
sha256_hash = hashlib.sha256()
|
|||
|
|
try:
|
|||
|
|
with open(filepath, "rb") as f:
|
|||
|
|
for byte_block in iter(lambda: f.read(4096), b""):
|
|||
|
|
sha256_hash.update(byte_block)
|
|||
|
|
return sha256_hash.hexdigest()
|
|||
|
|
except Exception:
|
|||
|
|
return ""
|
|||
|
|
|
|||
|
|
|
|||
|
|
def fastp_version(fastp_bin: Optional[Path] = None) -> str:
|
|||
|
|
if not fastp_bin:
|
|||
|
|
return "not_configured"
|
|||
|
|
try:
|
|||
|
|
result = subprocess.run(
|
|||
|
|
[str(fastp_bin), "--version"],
|
|||
|
|
capture_output=True,
|
|||
|
|
text=True,
|
|||
|
|
check=False,
|
|||
|
|
timeout=5,
|
|||
|
|
)
|
|||
|
|
output = result.stdout.strip() or result.stderr.strip()
|
|||
|
|
return output.split("\n")[0] if output else "unknown"
|
|||
|
|
except Exception:
|
|||
|
|
return "not_found"
|
|||
|
|
|
|||
|
|
|
|||
|
|
def resolve_fastp_bin(fastp_bin_arg: str = "") -> Optional[Path]:
|
|||
|
|
candidates = []
|
|||
|
|
if fastp_bin_arg:
|
|||
|
|
candidates.append(Path(fastp_bin_arg))
|
|||
|
|
env_fastp = os.environ.get("FASTP_BIN", "").strip()
|
|||
|
|
if env_fastp:
|
|||
|
|
candidates.append(Path(env_fastp))
|
|||
|
|
path_fastp = shutil.which("fastp")
|
|||
|
|
if path_fastp:
|
|||
|
|
candidates.append(Path(path_fastp))
|
|||
|
|
for candidate in candidates:
|
|||
|
|
if candidate.exists() and candidate.is_file() and os.access(candidate, os.X_OK):
|
|||
|
|
return candidate
|
|||
|
|
return None
|
|||
|
|
|
|||
|
|
|
|||
|
|
def get_tool_versions(fastp_bin: Optional[Path] = None) -> Dict[str, str]:
|
|||
|
|
"""获取所需工具版本信息"""
|
|||
|
|
versions = {}
|
|||
|
|
|
|||
|
|
# fastqc
|
|||
|
|
try:
|
|||
|
|
result = subprocess.run(
|
|||
|
|
[f"{RNASEQ_BIN}/fastqc", "--version"],
|
|||
|
|
capture_output=True,
|
|||
|
|
text=True,
|
|||
|
|
check=False,
|
|||
|
|
timeout=5
|
|||
|
|
)
|
|||
|
|
output = result.stdout.strip() or result.stderr.strip()
|
|||
|
|
versions["fastqc"] = output.split('\n')[0] if output else "unknown"
|
|||
|
|
except Exception:
|
|||
|
|
versions["fastqc"] = "not_found"
|
|||
|
|
|
|||
|
|
# cutadapt
|
|||
|
|
try:
|
|||
|
|
result = subprocess.run(
|
|||
|
|
[f"{RNASEQ_BIN}/cutadapt", "--version"],
|
|||
|
|
capture_output=True,
|
|||
|
|
text=True,
|
|||
|
|
check=False,
|
|||
|
|
timeout=5
|
|||
|
|
)
|
|||
|
|
output = result.stdout.strip() or result.stderr.strip()
|
|||
|
|
versions["cutadapt"] = output.split('\n')[0] if output else "unknown"
|
|||
|
|
except Exception:
|
|||
|
|
versions["cutadapt"] = "not_found"
|
|||
|
|
|
|||
|
|
# MultiQC
|
|||
|
|
try:
|
|||
|
|
result = subprocess.run(
|
|||
|
|
[f"{RNASEQ_BIN}/multiqc", "--version"],
|
|||
|
|
capture_output=True,
|
|||
|
|
text=True,
|
|||
|
|
check=False,
|
|||
|
|
timeout=10
|
|||
|
|
)
|
|||
|
|
output = result.stdout.strip() or result.stderr.strip()
|
|||
|
|
versions["multiqc"] = output.split('\n')[0] if output else "unknown"
|
|||
|
|
except Exception:
|
|||
|
|
versions["multiqc"] = "not_found"
|
|||
|
|
|
|||
|
|
versions["fastp"] = fastp_version(fastp_bin)
|
|||
|
|
|
|||
|
|
return versions
|
|||
|
|
|
|||
|
|
|
|||
|
|
def load_samples_config(config_file: Path) -> Dict[str, Any]:
|
|||
|
|
"""加载样本配置文件（支持合并格式：samples + groups + compare）"""
|
|||
|
|
with open(config_file, 'r') as f:
|
|||
|
|
config = yaml.safe_load(f)
|
|||
|
|
|
|||
|
|
# 确保有 samples 字段
|
|||
|
|
if "samples" not in config:
|
|||
|
|
raise ValueError(f"Config file {config_file} must contain 'samples' field")
|
|||
|
|
|
|||
|
|
return config
|
|||
|
|
|
|||
|
|
|
|||
|
|
RNASEQ_BIN = str(rnaseq_bin())
|
|||
|
|
prepend_path(RNASEQ_BIN)
|
|||
|
|
|
|||
|
|
|
|||
|
|
def resolve_thread_plan(cpu: int, sample_parallel_jobs: int = 1, threads_per_sample: Optional[int] = None) -> Dict[str, int]:
|
|||
|
|
"""Resolve conservative RNA-seq sample/thread settings from CLI args."""
|
|||
|
|
total_cpu = max(1, int(cpu or 1))
|
|||
|
|
jobs = max(1, int(sample_parallel_jobs or 1))
|
|||
|
|
if threads_per_sample is None or int(threads_per_sample) < 1:
|
|||
|
|
threads = max(1, total_cpu // jobs)
|
|||
|
|
else:
|
|||
|
|
threads = max(1, int(threads_per_sample))
|
|||
|
|
if jobs * threads > total_cpu:
|
|||
|
|
threads = max(1, total_cpu // jobs)
|
|||
|
|
return {"cpu": total_cpu, "sample_parallel_jobs": jobs, "threads_per_sample": threads}
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
def fastqc_output_base(fastq_path: Path) -> str:
|
|||
|
|
"""Return FastQC output basename without _fastqc suffix for .fastq/.fq(.gz)."""
|
|||
|
|
name = fastq_path.name
|
|||
|
|
for suffix in (".fastq.gz", ".fq.gz", ".fastq", ".fq"):
|
|||
|
|
if name.lower().endswith(suffix):
|
|||
|
|
return name[:-len(suffix)]
|
|||
|
|
return fastq_path.stem
|
|||
|
|
|
|||
|
|
|
|||
|
|
def sample_qc_outputs_complete(sample_name: str, file_info: Dict[str, Path], analysis_dir: Path, qc_tool: str = "fastp") -> bool:
|
|||
|
|
r1 = file_info["r1"]
|
|||
|
|
r2 = file_info["r2"]
|
|||
|
|
if qc_tool == "fastp":
|
|||
|
|
required = [
|
|||
|
|
analysis_dir / f"{sample_name}_1.clean.fastq.gz",
|
|||
|
|
analysis_dir / f"{sample_name}_2.clean.fastq.gz",
|
|||
|
|
analysis_dir / f"{sample_name}_fastp.html",
|
|||
|
|
analysis_dir / f"{sample_name}_fastp.json",
|
|||
|
|
analysis_dir / f"{sample_name}_fastp.log",
|
|||
|
|
]
|
|||
|
|
else:
|
|||
|
|
required = [
|
|||
|
|
analysis_dir / f"{fastqc_output_base(r1)}_fastqc.html",
|
|||
|
|
analysis_dir / f"{fastqc_output_base(r1)}_fastqc.zip",
|
|||
|
|
analysis_dir / f"{fastqc_output_base(r2)}_fastqc.html",
|
|||
|
|
analysis_dir / f"{fastqc_output_base(r2)}_fastqc.zip",
|
|||
|
|
analysis_dir / f"{sample_name}_1.clean.fastq.gz",
|
|||
|
|
analysis_dir / f"{sample_name}_2.clean.fastq.gz",
|
|||
|
|
analysis_dir / f"{sample_name}_cutadapt.log",
|
|||
|
|
]
|
|||
|
|
return all(path.exists() and path.stat().st_size > 0 for path in required)
|
|||
|
|
|
|||
|
|
def run_fastqc(fastq_files: List[Path], output_dir: Path, cpu: int = 4) -> Dict[str, Any]:
|
|||
|
|
"""
|
|||
|
|
运行FastQC质控
|
|||
|
|
|
|||
|
|
Returns:
|
|||
|
|
Dict with success status and timing info
|
|||
|
|
"""
|
|||
|
|
results = {"success": True, "files": {}, "total_time_sec": 0}
|
|||
|
|
fastqc_path = f"{RNASEQ_BIN}/fastqc"
|
|||
|
|
|
|||
|
|
for fq in fastq_files:
|
|||
|
|
cmd = [fastqc_path, "-o", str(output_dir), "-t", str(cpu), str(fq)]
|
|||
|
|
file_size_mb = fq.stat().st_size / 1024 / 1024
|
|||
|
|
|
|||
|
|
print(f" [FastQC] Starting: {fq.name} ({file_size_mb:.1f} MB)...", flush=True)
|
|||
|
|
start_ts = time.time()
|
|||
|
|
|
|||
|
|
try:
|
|||
|
|
result = subprocess.run(cmd, check=True, capture_output=True, text=True)
|
|||
|
|
elapsed = time.time() - start_ts
|
|||
|
|
results["files"][fq.name] = {"success": True, "time_sec": elapsed}
|
|||
|
|
results["total_time_sec"] += elapsed
|
|||
|
|
print(f" [FastQC] Completed: {fq.name} in {elapsed:.1f}s", flush=True)
|
|||
|
|
except subprocess.CalledProcessError as e:
|
|||
|
|
elapsed = time.time() - start_ts
|
|||
|
|
results["files"][fq.name] = {"success": False, "time_sec": elapsed, "error": str(e)}
|
|||
|
|
results["success"] = False
|
|||
|
|
print(f" [FastQC] FAILED: {fq.name} after {elapsed:.1f}s: {e.stderr}", file=sys.stderr)
|
|||
|
|
raise RuntimeError(f"FastQC failed for {fq.name}: {e.stderr}")
|
|||
|
|
|
|||
|
|
return results
|
|||
|
|
|
|||
|
|
|
|||
|
|
def run_cutadapt(
|
|||
|
|
r1: Path,
|
|||
|
|
r2: Path,
|
|||
|
|
out_r1: Path,
|
|||
|
|
out_r2: Path,
|
|||
|
|
log_file: Path,
|
|||
|
|
adapter_r1: str = "AGATCGGAAGAGC",
|
|||
|
|
adapter_r2: str = "AGATCGGAAGAGC",
|
|||
|
|
quality_cutoff: int = 20,
|
|||
|
|
min_length: int = 50,
|
|||
|
|
cpu: int = 1
|
|||
|
|
) -> Dict[str, Any]:
|
|||
|
|
"""
|
|||
|
|
运行Cutadapt去除接头
|
|||
|
|
|
|||
|
|
Returns:
|
|||
|
|
Dict with success status, timing info, and stats
|
|||
|
|
"""
|
|||
|
|
cmd = [
|
|||
|
|
f"{RNASEQ_BIN}/cutadapt",
|
|||
|
|
"-a", adapter_r1,
|
|||
|
|
"-A", adapter_r2,
|
|||
|
|
"-o", str(out_r1),
|
|||
|
|
"-p", str(out_r2),
|
|||
|
|
"-q", str(quality_cutoff),
|
|||
|
|
"-m", str(min_length),
|
|||
|
|
"-j", str(max(1, int(cpu or 1))),
|
|||
|
|
str(r1),
|
|||
|
|
str(r2)
|
|||
|
|
]
|
|||
|
|
|
|||
|
|
r1_size_mb = r1.stat().st_size / 1024 / 1024
|
|||
|
|
r2_size_mb = r2.stat().st_size / 1024 / 1024
|
|||
|
|
|
|||
|
|
print(f" [Cutadapt] Starting: {r1.name} + {r2.name} ({r1_size_mb+r2_size_mb:.1f} MB total)...", flush=True)
|
|||
|
|
start_ts = time.time()
|
|||
|
|
|
|||
|
|
try:
|
|||
|
|
result = subprocess.run(cmd, check=True, capture_output=True, text=True)
|
|||
|
|
elapsed = time.time() - start_ts
|
|||
|
|
|
|||
|
|
with open(log_file, 'w') as f:
|
|||
|
|
f.write(result.stdout)
|
|||
|
|
|
|||
|
|
# 验证输出文件
|
|||
|
|
if not out_r1.exists() or out_r1.stat().st_size == 0:
|
|||
|
|
raise RuntimeError(f"Cutadapt output R1 is empty or missing: {out_r1}")
|
|||
|
|
if not out_r2.exists() or out_r2.stat().st_size == 0:
|
|||
|
|
raise RuntimeError(f"Cutadapt output R2 is empty or missing: {out_r2}")
|
|||
|
|
|
|||
|
|
out_size_mb = (out_r1.stat().st_size + out_r2.stat().st_size) / 1024 / 1024
|
|||
|
|
print(f" [Cutadapt] Completed in {elapsed:.1f}s (output: {out_size_mb:.1f} MB)", flush=True)
|
|||
|
|
|
|||
|
|
return {"success": True, "time_sec": elapsed, "log": result.stdout}
|
|||
|
|
|
|||
|
|
except subprocess.CalledProcessError as e:
|
|||
|
|
elapsed = time.time() - start_ts
|
|||
|
|
with open(log_file, 'w') as f:
|
|||
|
|
f.write(e.stdout + "\n" + e.stderr)
|
|||
|
|
print(f" [Cutadapt] FAILED after {elapsed:.1f}s: {e.stderr}", file=sys.stderr)
|
|||
|
|
raise RuntimeError(f"Cutadapt failed: {e.stderr}")
|
|||
|
|
|
|||
|
|
except Exception as e:
|
|||
|
|
elapsed = time.time() - start_ts
|
|||
|
|
print(f" [Cutadapt] FAILED after {elapsed:.1f}s: {str(e)}", file=sys.stderr)
|
|||
|
|
raise
|
|||
|
|
|
|||
|
|
|
|||
|
|
def run_fastp(
|
|||
|
|
fastp_bin: Path,
|
|||
|
|
r1: Path,
|
|||
|
|
r2: Path,
|
|||
|
|
out_r1: Path,
|
|||
|
|
out_r2: Path,
|
|||
|
|
html_report: Path,
|
|||
|
|
json_report: Path,
|
|||
|
|
log_file: Path,
|
|||
|
|
cpu: int = 1,
|
|||
|
|
detect_adapter_for_pe: bool = True,
|
|||
|
|
cut_front: bool = True,
|
|||
|
|
cut_tail: bool = True,
|
|||
|
|
cut_window_size: int = 4,
|
|||
|
|
cut_mean_quality: int = 20,
|
|||
|
|
qualified_quality_phred: int = 20,
|
|||
|
|
unqualified_percent_limit: int = 30,
|
|||
|
|
n_base_limit: int = 5,
|
|||
|
|
length_required: int = 30,
|
|||
|
|
) -> Dict[str, Any]:
|
|||
|
|
cmd = [
|
|||
|
|
str(fastp_bin),
|
|||
|
|
"-i", str(r1),
|
|||
|
|
"-I", str(r2),
|
|||
|
|
"-o", str(out_r1),
|
|||
|
|
"-O", str(out_r2),
|
|||
|
|
"--html", str(html_report),
|
|||
|
|
"--json", str(json_report),
|
|||
|
|
"--thread", str(max(1, int(cpu or 1))),
|
|||
|
|
"--cut_window_size", str(cut_window_size),
|
|||
|
|
"--cut_mean_quality", str(cut_mean_quality),
|
|||
|
|
"--qualified_quality_phred", str(qualified_quality_phred),
|
|||
|
|
"--unqualified_percent_limit", str(unqualified_percent_limit),
|
|||
|
|
"--n_base_limit", str(n_base_limit),
|
|||
|
|
"--length_required", str(length_required),
|
|||
|
|
]
|
|||
|
|
if detect_adapter_for_pe:
|
|||
|
|
cmd.append("--detect_adapter_for_pe")
|
|||
|
|
if cut_front:
|
|||
|
|
cmd.append("--cut_front")
|
|||
|
|
if cut_tail:
|
|||
|
|
cmd.append("--cut_tail")
|
|||
|
|
|
|||
|
|
print(f" [fastp] Starting: {r1.name} + {r2.name}...", flush=True)
|
|||
|
|
start_ts = time.time()
|
|||
|
|
result = subprocess.run(cmd, capture_output=True, text=True, check=False)
|
|||
|
|
elapsed = time.time() - start_ts
|
|||
|
|
log_file.write_text(
|
|||
|
|
"Command: " + " ".join(cmd) + "\n\n"
|
|||
|
|
+ "--- STDOUT ---\n"
|
|||
|
|
+ result.stdout
|
|||
|
|
+ "\n--- STDERR ---\n"
|
|||
|
|
+ result.stderr,
|
|||
|
|
encoding="utf-8",
|
|||
|
|
)
|
|||
|
|
if result.returncode != 0:
|
|||
|
|
print(f" [fastp] FAILED after {elapsed:.1f}s: {result.stderr}", file=sys.stderr)
|
|||
|
|
raise RuntimeError(f"fastp failed with exit code {result.returncode}: {result.stderr}")
|
|||
|
|
for path in (out_r1, out_r2, html_report, json_report):
|
|||
|
|
if not path.exists() or path.stat().st_size == 0:
|
|||
|
|
raise RuntimeError(f"fastp output is missing or empty: {path}")
|
|||
|
|
print(f" [fastp] Completed in {elapsed:.1f}s", flush=True)
|
|||
|
|
return {"success": True, "time_sec": elapsed, "cmd": cmd, "log": str(log_file)}
|
|||
|
|
|
|||
|
|
|
|||
|
|
def run_multiqc(search_paths: Sequence[Path], output_dir: Path, report_name: str) -> Optional[Path]:
|
|||
|
|
"""Run MultiQC over selected QC outputs and return the report path."""
|
|||
|
|
multiqc = Path(f"{RNASEQ_BIN}/multiqc")
|
|||
|
|
existing = [str(path) for path in search_paths if path.exists()]
|
|||
|
|
if not multiqc.exists() or not existing:
|
|||
|
|
return None
|
|||
|
|
cmd = [str(multiqc), *existing, "-o", str(output_dir), "-n", report_name, "-f"]
|
|||
|
|
try:
|
|||
|
|
subprocess.run(cmd, check=True, capture_output=True, text=True)
|
|||
|
|
report = output_dir / report_name
|
|||
|
|
return report if report.exists() else None
|
|||
|
|
except subprocess.CalledProcessError as e:
|
|||
|
|
print(f"MultiQC failed: {e.stderr}", file=sys.stderr)
|
|||
|
|
return None
|
|||
|
|
|
|||
|
|
|
|||
|
|
def get_sample_stats(sample_dir: Path, sample_name: str) -> Dict[str, Any]:
|
|||
|
|
"""获取样本统计信息"""
|
|||
|
|
import re
|
|||
|
|
|
|||
|
|
def first_count_after_colon(text: str) -> Optional[int]:
|
|||
|
|
value = text.split(":", 1)[1] if ":" in text else text
|
|||
|
|
match = re.search(r'([\d,]+)', value)
|
|||
|
|
if match:
|
|||
|
|
return int(match.group(1).replace(',', ''))
|
|||
|
|
return None
|
|||
|
|
|
|||
|
|
stats = {
|
|||
|
|
"raw_reads": 0,
|
|||
|
|
"clean_reads": 0,
|
|||
|
|
"adapter_removed": 0,
|
|||
|
|
"quality_filtered": 0
|
|||
|
|
}
|
|||
|
|
fastp_json = sample_dir / f"{sample_name}_fastp.json"
|
|||
|
|
if fastp_json.exists():
|
|||
|
|
try:
|
|||
|
|
payload = json.loads(fastp_json.read_text(encoding="utf-8"))
|
|||
|
|
before = payload.get("summary", {}).get("before_filtering", {})
|
|||
|
|
after = payload.get("summary", {}).get("after_filtering", {})
|
|||
|
|
stats["raw_reads"] = int(before.get("total_reads", 0))
|
|||
|
|
stats["clean_reads"] = int(after.get("total_reads", 0))
|
|||
|
|
stats["fastp_q20_rate_before"] = before.get("q20_rate")
|
|||
|
|
stats["fastp_q30_rate_before"] = before.get("q30_rate")
|
|||
|
|
stats["fastp_q20_rate_after"] = after.get("q20_rate")
|
|||
|
|
stats["fastp_q30_rate_after"] = after.get("q30_rate")
|
|||
|
|
stats["read1_mean_length"] = after.get("read1_mean_length")
|
|||
|
|
stats["read2_mean_length"] = after.get("read2_mean_length")
|
|||
|
|
stats["gc_content_after"] = after.get("gc_content")
|
|||
|
|
except Exception as e:
|
|||
|
|
print(f"Warning: Could not read fastp JSON for {sample_name}: {e}")
|
|||
|
|
|
|||
|
|
# 读取FastQC数据（从zip文件中提取）
|
|||
|
|
fastqc_zip = sample_dir / f"{sample_name}_1_fastqc.zip"
|
|||
|
|
if fastqc_zip.exists():
|
|||
|
|
try:
|
|||
|
|
import zipfile
|
|||
|
|
with zipfile.ZipFile(fastqc_zip, 'r') as zf:
|
|||
|
|
# 找到 fastqc_data.txt 文件
|
|||
|
|
data_file = None
|
|||
|
|
for name in zf.namelist():
|
|||
|
|
if name.endswith('fastqc_data.txt'):
|
|||
|
|
data_file = name
|
|||
|
|
break
|
|||
|
|
if data_file:
|
|||
|
|
with zf.open(data_file) as f:
|
|||
|
|
for line in f:
|
|||
|
|
line = line.decode('utf-8')
|
|||
|
|
if line.startswith("Total Sequences"):
|
|||
|
|
stats["raw_reads"] = int(line.split("\t")[1])
|
|||
|
|
elif line.startswith("Sequence length"):
|
|||
|
|
stats["read_length"] = line.split("\t")[1].strip()
|
|||
|
|
except Exception as e:
|
|||
|
|
print(f"Warning: Could not read FastQC data for {sample_name}: {e}")
|
|||
|
|
|
|||
|
|
# 读取Cutadapt日志
|
|||
|
|
log_file = sample_dir / f"{sample_name}_cutadapt.log"
|
|||
|
|
if log_file.exists():
|
|||
|
|
with open(log_file, 'r') as f:
|
|||
|
|
for line in f:
|
|||
|
|
if "Total reads processed:" in line or "Total read pairs processed:" in line:
|
|||
|
|
count = first_count_after_colon(line)
|
|||
|
|
if count is not None:
|
|||
|
|
stats["raw_reads"] = count
|
|||
|
|
elif "Reads with adapters:" in line or "Read 1 with adapter:" in line or "Read 2 with adapter:" in line:
|
|||
|
|
count = first_count_after_colon(line)
|
|||
|
|
if count is not None:
|
|||
|
|
stats["adapter_removed"] += count
|
|||
|
|
elif "Pairs that were too short:" in line or "Reads that were too short:" in line:
|
|||
|
|
count = first_count_after_colon(line)
|
|||
|
|
if count is not None:
|
|||
|
|
stats["quality_filtered"] += count
|
|||
|
|
elif "Reads written (passing filters):" in line or "Pairs written (passing filters):" in line:
|
|||
|
|
count = first_count_after_colon(line)
|
|||
|
|
if count is not None:
|
|||
|
|
stats["clean_reads"] = count
|
|||
|
|
|
|||
|
|
return stats
|
|||
|
|
|
|||
|
|
|
|||
|
|
def generate_report(
|
|||
|
|
output_dir: Path,
|
|||
|
|
sample_stats: Dict[str, Dict[str, Any]],
|
|||
|
|
versions: Dict[str, str],
|
|||
|
|
start_time: datetime,
|
|||
|
|
end_time: datetime,
|
|||
|
|
prefix: str = "",
|
|||
|
|
qc_tool: str = "fastp",
|
|||
|
|
) -> Path:
|
|||
|
|
"""生成质控报告"""
|
|||
|
|
report_file = output_dir / f"{prefix}.rnaseq_qc_report.md" if prefix else output_dir / "rnaseq_qc_report.md"
|
|||
|
|
|
|||
|
|
elapsed = (end_time - start_time).total_seconds()
|
|||
|
|
|
|||
|
|
content = f"""# RNA-Seq Quality Control Report
|
|||
|
|
|
|||
|
|
## Summary
|
|||
|
|
|
|||
|
|
- **Analysis Date**: {end_time.strftime("%Y-%m-%d %H:%M:%S")}
|
|||
|
|
- **Total Samples**: {len(sample_stats)}
|
|||
|
|
- **Total Runtime**: {elapsed:.1f} seconds
|
|||
|
|
- **QC Tool**: {qc_tool}
|
|||
|
|
|
|||
|
|
## Tool Versions
|
|||
|
|
|
|||
|
|
"""
|
|||
|
|
for tool, version in versions.items():
|
|||
|
|
content += f"- **{tool}**: {version}\n"
|
|||
|
|
|
|||
|
|
content += "\n## Sample Statistics\n\n"
|
|||
|
|
content += "| Sample | Raw Reads | Clean Reads | Adapter Removed |\n"
|
|||
|
|
content += "|--------|-----------|-------------|-----------------|\n"
|
|||
|
|
|
|||
|
|
for sample, stats in sample_stats.items():
|
|||
|
|
content += f"| {sample} | {stats.get('raw_reads', 'N/A')} | {stats.get('clean_reads', 'N/A')} | {stats.get('adapter_removed', 'N/A')} |\n"
|
|||
|
|
|
|||
|
|
with open(report_file, 'w') as f:
|
|||
|
|
f.write(content)
|
|||
|
|
|
|||
|
|
return report_file
|
|||
|
|
|
|||
|
|
|
|||
|
|
def generate_manifest(
|
|||
|
|
output_dir: Path,
|
|||
|
|
sample_stats: Dict[str, Dict[str, Any]],
|
|||
|
|
versions: Dict[str, str],
|
|||
|
|
start_time: datetime,
|
|||
|
|
end_time: datetime,
|
|||
|
|
cmd: List[str],
|
|||
|
|
input_dir: Path,
|
|||
|
|
workspace_dir: Path,
|
|||
|
|
job_id: str,
|
|||
|
|
run_id: str,
|
|||
|
|
prefix: str = "",
|
|||
|
|
analysis_dir: Path = None, # analysis 目录路径
|
|||
|
|
cpu: int = 4,
|
|||
|
|
sample_parallel_jobs: int = 1,
|
|||
|
|
threads_per_sample: int = 4,
|
|||
|
|
sample_config: Optional[Path] = None,
|
|||
|
|
qc_tool: str = "fastp",
|
|||
|
|
fastp_bin: str = "",
|
|||
|
|
fastp_params: Optional[Dict[str, Any]] = None,
|
|||
|
|
) -> Path:
|
|||
|
|
"""生成manifest文件"""
|
|||
|
|
manifest_file = output_dir / f"{prefix}.rnaseq_qc.manifest.json" if prefix else output_dir / "rnaseq_qc.manifest.json"
|
|||
|
|
|
|||
|
|
total_raw = sum(s.get('raw_reads', 0) for s in sample_stats.values())
|
|||
|
|
total_clean = sum(s.get('clean_reads', 0) for s in sample_stats.values())
|
|||
|
|
|
|||
|
|
manifest = {
|
|||
|
|
"manifest_version": "1.0",
|
|||
|
|
"stage": "rnaseq_qc",
|
|||
|
|
"status": "success",
|
|||
|
|
"prefix": prefix,
|
|||
|
|
"job_id": job_id,
|
|||
|
|
"run_id": run_id,
|
|||
|
|
"timestamp": end_time.isoformat(),
|
|||
|
|
"start_time": start_time.isoformat(),
|
|||
|
|
"end_time": end_time.isoformat(),
|
|||
|
|
"elapsed_sec": (end_time - start_time).total_seconds(),
|
|||
|
|
"cmd": cmd,
|
|||
|
|
"workspace_root": str(output_dir.parent.parent),
|
|||
|
|
"workspace": str(workspace_dir),
|
|||
|
|
"analysis_dir": str(analysis_dir) if analysis_dir else str(workspace_dir / "analysis"),
|
|||
|
|
"input_dir": str(input_dir),
|
|||
|
|
"metrics": {
|
|||
|
|
"total_samples": len(sample_stats),
|
|||
|
|
"total_raw_reads": total_raw,
|
|||
|
|
"total_clean_reads": total_clean,
|
|||
|
|
"samples": sample_stats
|
|||
|
|
},
|
|||
|
|
"parameters": {
|
|||
|
|
"cpu": cpu,
|
|||
|
|
"sample_parallel_jobs": sample_parallel_jobs,
|
|||
|
|
"threads_per_sample": threads_per_sample,
|
|||
|
|
"qc_tool": qc_tool,
|
|||
|
|
"fastqc_params": f"-t {threads_per_sample}",
|
|||
|
|
"cutadapt_params": f"-q 20 -m 50 -j {threads_per_sample}" if qc_tool == "cutadapt" else "",
|
|||
|
|
"fastp_bin": fastp_bin,
|
|||
|
|
"fastp_params": fastp_params or {},
|
|||
|
|
},
|
|||
|
|
"artifacts": {},
|
|||
|
|
"tools": versions
|
|||
|
|
}
|
|||
|
|
context = build_rnaseq_context(
|
|||
|
|
prefix=prefix,
|
|||
|
|
sample_config=sample_config,
|
|||
|
|
foreground_mode="global",
|
|||
|
|
stage="rnaseq_qc",
|
|||
|
|
)
|
|||
|
|
apply_rnaseq_context(manifest, context)
|
|||
|
|
|
|||
|
|
# 添加artifacts（文件生成在 final/ 目录，使用 prefix 命名）
|
|||
|
|
# 参考 synteny 契约，包含 full_path、size、sha256 等字段
|
|||
|
|
report_name = f"{prefix}.rnaseq_qc_report.md" if prefix else "rnaseq_qc_report.md"
|
|||
|
|
manifest_name = f"{prefix}.rnaseq_qc.manifest.json" if prefix else "rnaseq_qc.manifest.json"
|
|||
|
|
|
|||
|
|
def add_artifact(name: str, filepath: Path, artifact_type: str, **extra):
|
|||
|
|
"""添加 artifact 到 manifest"""
|
|||
|
|
if filepath.exists():
|
|||
|
|
manifest["artifacts"][name] = {
|
|||
|
|
"artifact": name,
|
|||
|
|
"path": f"final/{filepath.name}",
|
|||
|
|
"full_path": str(filepath.absolute()),
|
|||
|
|
"type": artifact_type,
|
|||
|
|
"size": filepath.stat().st_size,
|
|||
|
|
"sha256": compute_file_hash(filepath),
|
|||
|
|
**extra
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
# analysis_dir 用于引用质控报告文件（保留在 analysis 目录）
|
|||
|
|
analysis_path = analysis_dir if analysis_dir else workspace_dir / "analysis"
|
|||
|
|
|
|||
|
|
# final_dir 是 output_dir 的父目录下的 final 目录
|
|||
|
|
final_dir = output_dir.parent / "final"
|
|||
|
|
|
|||
|
|
artifact_prefix = prefix if prefix else "rnaseq"
|
|||
|
|
for sample in sample_stats.keys():
|
|||
|
|
# FastQC HTML 报告（复制到 final，前端展示）
|
|||
|
|
html_r1 = final_dir / f"{artifact_prefix}.rnaseq_qc.{sample}_1_fastqc.html"
|
|||
|
|
html_r2 = final_dir / f"{artifact_prefix}.rnaseq_qc.{sample}_2_fastqc.html"
|
|||
|
|
add_artifact(f"qc_report_{sample}_r1", html_r1, "html", format="html")
|
|||
|
|
add_artifact(f"qc_report_{sample}_r2", html_r2, "html", format="html")
|
|||
|
|
|
|||
|
|
# FastQC 数据（zip，保留在 analysis 目录）
|
|||
|
|
zip_r1 = analysis_path / f"{sample}_1_fastqc.zip"
|
|||
|
|
zip_r2 = analysis_path / f"{sample}_2_fastqc.zip"
|
|||
|
|
if zip_r1.exists():
|
|||
|
|
manifest["artifacts"][f"qc_data_{sample}_r1"] = {
|
|||
|
|
"artifact": f"qc_data_{sample}_r1",
|
|||
|
|
"path": f"rnaseq_qc/analysis/{zip_r1.name}",
|
|||
|
|
"full_path": str(zip_r1.absolute()),
|
|||
|
|
"type": "zip",
|
|||
|
|
"size": zip_r1.stat().st_size,
|
|||
|
|
"sha256": compute_file_hash(zip_r1),
|
|||
|
|
"format": "zip"
|
|||
|
|
}
|
|||
|
|
if zip_r2.exists():
|
|||
|
|
manifest["artifacts"][f"qc_data_{sample}_r2"] = {
|
|||
|
|
"artifact": f"qc_data_{sample}_r2",
|
|||
|
|
"path": f"rnaseq_qc/analysis/{zip_r2.name}",
|
|||
|
|
"full_path": str(zip_r2.absolute()),
|
|||
|
|
"type": "zip",
|
|||
|
|
"size": zip_r2.stat().st_size,
|
|||
|
|
"sha256": compute_file_hash(zip_r2),
|
|||
|
|
"format": "zip"
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
# Clean reads（FASTQ）- 大文件，复制到 final 后使用标准前缀命名
|
|||
|
|
clean_r1 = output_dir / f"{artifact_prefix}.rnaseq_qc.{sample}_1.clean.fastq.gz"
|
|||
|
|
clean_r2 = output_dir / f"{artifact_prefix}.rnaseq_qc.{sample}_2.clean.fastq.gz"
|
|||
|
|
add_artifact(f"clean_reads_r1_{sample}", clean_r1, "fastq", format="fastq.gz")
|
|||
|
|
add_artifact(f"clean_reads_r2_{sample}", clean_r2, "fastq", format="fastq.gz")
|
|||
|
|
|
|||
|
|
# Cutadapt 日志
|
|||
|
|
log_file = output_dir / f"{artifact_prefix}.rnaseq_qc.{sample}_cutadapt.log"
|
|||
|
|
add_artifact(f"cutadapt_log_{sample}", log_file, "log", format="txt")
|
|||
|
|
|
|||
|
|
fastp_html = output_dir / f"{artifact_prefix}.rnaseq_qc.{sample}_fastp.html"
|
|||
|
|
fastp_json = output_dir / f"{artifact_prefix}.rnaseq_qc.{sample}_fastp.json"
|
|||
|
|
fastp_log = output_dir / f"{artifact_prefix}.rnaseq_qc.{sample}_fastp.log"
|
|||
|
|
add_artifact(f"fastp_report_{sample}", fastp_html, "html", format="html")
|
|||
|
|
add_artifact(f"fastp_json_{sample}", fastp_json, "json", format="json")
|
|||
|
|
add_artifact(f"fastp_log_{sample}", fastp_log, "log", format="txt")
|
|||
|
|
|
|||
|
|
# 报告、MultiQC 和 manifest
|
|||
|
|
report_file = output_dir / report_name
|
|||
|
|
add_artifact("report", report_file, "markdown", format="md")
|
|||
|
|
multiqc_name = f"{prefix}.rnaseq_qc_multiqc_report.html" if prefix else "rnaseq_qc_multiqc_report.html"
|
|||
|
|
multiqc_file = output_dir / multiqc_name
|
|||
|
|
add_artifact("multiqc_report", multiqc_file, "html", format="html")
|
|||
|
|
multiqc_data = output_dir / multiqc_name.replace(".html", "_data")
|
|||
|
|
if multiqc_data.exists():
|
|||
|
|
manifest["artifacts"]["multiqc_data"] = {
|
|||
|
|
"artifact": "multiqc_data",
|
|||
|
|
"path": f"final/{multiqc_data.name}",
|
|||
|
|
"full_path": str(multiqc_data.absolute()),
|
|||
|
|
"type": "directory",
|
|||
|
|
"format": "multiqc_data"
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
# manifest 自身（在生成后添加）
|
|||
|
|
manifest["artifacts"]["manifest"] = {
|
|||
|
|
"artifact": "manifest",
|
|||
|
|
"path": f"final/{manifest_name}",
|
|||
|
|
"type": "json",
|
|||
|
|
"format": "json"
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
with open(manifest_file, 'w') as f:
|
|||
|
|
json.dump(manifest, f, indent=2)
|
|||
|
|
|
|||
|
|
return manifest_file
|
|||
|
|
|
|||
|
|
|
|||
|
|
def main():
|
|||
|
|
parser = argparse.ArgumentParser(description="RNA-Seq Quality Control")
|
|||
|
|
parser.add_argument("-i", "--input", required=True, help="Input directory with FASTQ files")
|
|||
|
|
parser.add_argument("-c", "--config", required=True, help="Sample configuration file (YAML)")
|
|||
|
|
parser.add_argument("-o", "--output", required=True, help="Output directory")
|
|||
|
|
parser.add_argument("-j", "--job-id", default="unknown", help="Job ID")
|
|||
|
|
parser.add_argument("-r", "--run-id", default="unknown", help="Run ID")
|
|||
|
|
parser.add_argument("-t", "--cpu", type=int, default=4, help="Total CPU threads for this stage")
|
|||
|
|
parser.add_argument("--sample-parallel-jobs", type=int, default=1, help="Number of samples to process concurrently")
|
|||
|
|
parser.add_argument("--threads-per-sample", type=int, help="Threads used by each per-sample tool invocation")
|
|||
|
|
parser.add_argument("--resume-existing", action="store_true", help="Skip samples whose FastQC, clean reads, and cutadapt outputs already exist in the output directory")
|
|||
|
|
parser.add_argument("-p", "--prefix", default="", help="Output file prefix")
|
|||
|
|
parser.add_argument("--adapter-r1", default="AGATCGGAAGAGC", help="Adapter sequence for read1")
|
|||
|
|
parser.add_argument("--adapter-r2", default="AGATCGGAAGAGC", help="Adapter sequence for read2")
|
|||
|
|
parser.add_argument("--qc-tool", choices=["cutadapt", "fastp"], default="fastp", help="QC/cleaning backend. Default uses fastp; cutadapt remains available by explicit selection.")
|
|||
|
|
parser.add_argument("--fastp-bin", default="tools/fastp", help="Path to fastp executable. If omitted, FASTP_BIN or PATH is used.")
|
|||
|
|
parser.add_argument("--detect-adapter-for-pe", action=argparse.BooleanOptionalAction, default=True, help="Enable fastp paired-end adapter auto-detection.")
|
|||
|
|
parser.add_argument("--cut-front", action=argparse.BooleanOptionalAction, default=True, help="Enable fastp front quality cutting.")
|
|||
|
|
parser.add_argument("--cut-tail", action=argparse.BooleanOptionalAction, default=True, help="Enable fastp tail quality cutting.")
|
|||
|
|
parser.add_argument("--cut-window-size", type=int, default=4, help="fastp quality cutting window size.")
|
|||
|
|
parser.add_argument("--cut-mean-quality", type=int, default=20, help="fastp mean quality threshold for cutting.")
|
|||
|
|
parser.add_argument("--qualified-quality-phred", type=int, default=20, help="fastp qualified quality phred threshold.")
|
|||
|
|
parser.add_argument("--unqualified-percent-limit", type=int, default=30, help="fastp maximum unqualified base percent.")
|
|||
|
|
parser.add_argument("--n-base-limit", type=int, default=5, help="fastp maximum N bases per read.")
|
|||
|
|
parser.add_argument("--length-required", type=int, default=30, help="fastp minimum output read length.")
|
|||
|
|
|
|||
|
|
args = parser.parse_args()
|
|||
|
|
thread_plan = resolve_thread_plan(args.cpu, args.sample_parallel_jobs, args.threads_per_sample)
|
|||
|
|
threads_per_sample = thread_plan["threads_per_sample"]
|
|||
|
|
|
|||
|
|
start_time = datetime.now()
|
|||
|
|
input_dir = Path(args.input)
|
|||
|
|
output_dir = Path(args.output) # rnaseq_qc 工作目录
|
|||
|
|
config_file = Path(args.config)
|
|||
|
|
context_for_prefix = resolve_rnaseq_project_context(
|
|||
|
|
input_dir=input_dir,
|
|||
|
|
requested_prefix=args.prefix,
|
|||
|
|
sample_config=config_file,
|
|||
|
|
)
|
|||
|
|
resolved_prefix = context_for_prefix.get("output_prefix")
|
|||
|
|
if resolved_prefix and resolved_prefix != args.prefix:
|
|||
|
|
print(
|
|||
|
|
f"[INFO] Using RNA-seq output prefix {resolved_prefix} "
|
|||
|
|
f"instead of requested {args.prefix or '<empty>'} "
|
|||
|
|
f"(source={context_for_prefix.get('source', '')})"
|
|||
|
|
)
|
|||
|
|
args.prefix = resolved_prefix
|
|||
|
|
|
|||
|
|
# 创建工作目录
|
|||
|
|
output_dir.mkdir(parents=True, exist_ok=True)
|
|||
|
|
|
|||
|
|
# 最终结果目录: workspace/<job_id>/final/
|
|||
|
|
# 根据契约，最终产物必须位于 workspace/<job_id>/final/
|
|||
|
|
workspace_dir = output_dir.parent # job_id 目录
|
|||
|
|
final_dir = workspace_dir / "final"
|
|||
|
|
final_dir.mkdir(exist_ok=True)
|
|||
|
|
|
|||
|
|
# 中间文件目录: rnaseq_qc/analysis/
|
|||
|
|
analysis_dir = output_dir / "analysis"
|
|||
|
|
analysis_dir.mkdir(exist_ok=True)
|
|||
|
|
|
|||
|
|
# 加载样本配置并严格验证输入文件
|
|||
|
|
config = load_samples_config(config_file)
|
|||
|
|
|
|||
|
|
# 严格验证输入文件（失败时直接退出）
|
|||
|
|
validated_samples = validate_input_files(input_dir, config)
|
|||
|
|
|
|||
|
|
print(f"\nProcessing {len(validated_samples)} samples...")
|
|||
|
|
|
|||
|
|
fastp_bin = resolve_fastp_bin(args.fastp_bin) if args.qc_tool == "fastp" else None
|
|||
|
|
if args.qc_tool == "fastp" and not fastp_bin:
|
|||
|
|
print("❌ ERROR: --qc-tool fastp requires an executable --fastp-bin, FASTP_BIN, or fastp on PATH", file=sys.stderr)
|
|||
|
|
sys.exit(1)
|
|||
|
|
|
|||
|
|
fastp_params = {
|
|||
|
|
"detect_adapter_for_pe": bool(args.detect_adapter_for_pe),
|
|||
|
|
"cut_front": bool(args.cut_front),
|
|||
|
|
"cut_tail": bool(args.cut_tail),
|
|||
|
|
"cut_window_size": args.cut_window_size,
|
|||
|
|
"cut_mean_quality": args.cut_mean_quality,
|
|||
|
|
"qualified_quality_phred": args.qualified_quality_phred,
|
|||
|
|
"unqualified_percent_limit": args.unqualified_percent_limit,
|
|||
|
|
"n_base_limit": args.n_base_limit,
|
|||
|
|
"length_required": args.length_required,
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
# 获取工具版本
|
|||
|
|
versions = get_tool_versions(fastp_bin)
|
|||
|
|
print(f"Tool versions: {versions}")
|
|||
|
|
|
|||
|
|
sample_stats = {}
|
|||
|
|
total_samples = len(validated_samples)
|
|||
|
|
|
|||
|
|
def qc_one_sample(sample_name: str, file_info: Dict[str, Path]) -> tuple[str, Optional[Dict[str, Any]], Optional[str]]:
|
|||
|
|
r1 = file_info["r1"]
|
|||
|
|
r2 = file_info["r2"]
|
|||
|
|
print(f"\n[PROGRESS] Processing: {sample_name}", flush=True)
|
|||
|
|
print(f" Input: {r1.name} + {r2.name}", flush=True)
|
|||
|
|
if args.resume_existing and sample_qc_outputs_complete(sample_name, file_info, analysis_dir, args.qc_tool):
|
|||
|
|
stats = get_sample_stats(analysis_dir, sample_name)
|
|||
|
|
stats["resumed_from_existing"] = True
|
|||
|
|
stats["fastqc_time_sec"] = 0
|
|||
|
|
stats["cutadapt_time_sec"] = 0
|
|||
|
|
stats["fastp_time_sec"] = 0
|
|||
|
|
print(f" [RESUME] Skipping completed sample {sample_name}; existing QC outputs are complete", flush=True)
|
|||
|
|
return sample_name, stats, None
|
|||
|
|
|
|||
|
|
out_r1 = analysis_dir / f"{sample_name}_1.clean.fastq.gz"
|
|||
|
|
out_r2 = analysis_dir / f"{sample_name}_2.clean.fastq.gz"
|
|||
|
|
if args.qc_tool == "fastp":
|
|||
|
|
try:
|
|||
|
|
fastp_result = run_fastp(
|
|||
|
|
fastp_bin=fastp_bin,
|
|||
|
|
r1=r1,
|
|||
|
|
r2=r2,
|
|||
|
|
out_r1=out_r1,
|
|||
|
|
out_r2=out_r2,
|
|||
|
|
html_report=analysis_dir / f"{sample_name}_fastp.html",
|
|||
|
|
json_report=analysis_dir / f"{sample_name}_fastp.json",
|
|||
|
|
log_file=analysis_dir / f"{sample_name}_fastp.log",
|
|||
|
|
cpu=threads_per_sample,
|
|||
|
|
**fastp_params,
|
|||
|
|
)
|
|||
|
|
except RuntimeError as e:
|
|||
|
|
return sample_name, None, f"fastp failed for {sample_name}: {e}"
|
|||
|
|
stats = get_sample_stats(analysis_dir, sample_name)
|
|||
|
|
stats["fastqc_time_sec"] = 0
|
|||
|
|
stats["cutadapt_time_sec"] = 0
|
|||
|
|
stats["fastp_time_sec"] = fastp_result.get("time_sec", 0)
|
|||
|
|
else:
|
|||
|
|
try:
|
|||
|
|
fastqc_result = run_fastqc([r1, r2], analysis_dir, threads_per_sample)
|
|||
|
|
except RuntimeError as e:
|
|||
|
|
return sample_name, None, f"FastQC failed for {sample_name}: {e}"
|
|||
|
|
log_file = analysis_dir / f"{sample_name}_cutadapt.log"
|
|||
|
|
try:
|
|||
|
|
cutadapt_result = run_cutadapt(
|
|||
|
|
r1, r2, out_r1, out_r2, log_file,
|
|||
|
|
args.adapter_r1, args.adapter_r2,
|
|||
|
|
cpu=threads_per_sample,
|
|||
|
|
)
|
|||
|
|
except RuntimeError as e:
|
|||
|
|
return sample_name, None, f"Cutadapt failed for {sample_name}: {e}"
|
|||
|
|
stats = get_sample_stats(analysis_dir, sample_name)
|
|||
|
|
stats["fastqc_time_sec"] = fastqc_result.get("total_time_sec", 0)
|
|||
|
|
stats["cutadapt_time_sec"] = cutadapt_result.get("time_sec", 0)
|
|||
|
|
stats["fastp_time_sec"] = 0
|
|||
|
|
|
|||
|
|
print(f" ✓ Sample {sample_name} completed", flush=True)
|
|||
|
|
return sample_name, stats, None
|
|||
|
|
|
|||
|
|
print(
|
|||
|
|
f"Thread plan: sample_parallel_jobs={thread_plan['sample_parallel_jobs']}, "
|
|||
|
|
f"threads_per_sample={threads_per_sample}, total_cpu={thread_plan['cpu']}"
|
|||
|
|
)
|
|||
|
|
with ThreadPoolExecutor(max_workers=max(1, int(thread_plan["sample_parallel_jobs"]))) as executor:
|
|||
|
|
future_map = {
|
|||
|
|
executor.submit(qc_one_sample, sample_name, file_info): sample_name
|
|||
|
|
for sample_name, file_info in validated_samples.items()
|
|||
|
|
}
|
|||
|
|
for future in as_completed(future_map):
|
|||
|
|
sample_name = future_map[future]
|
|||
|
|
try:
|
|||
|
|
completed_name, stats, error = future.result()
|
|||
|
|
except Exception as exc:
|
|||
|
|
print(f"\n❌ ERROR: QC failed for {sample_name}: {exc}", file=sys.stderr)
|
|||
|
|
sys.exit(1)
|
|||
|
|
if error:
|
|||
|
|
print(f"\n❌ ERROR: {error}", file=sys.stderr)
|
|||
|
|
sys.exit(1)
|
|||
|
|
if stats is not None:
|
|||
|
|
sample_stats[completed_name] = stats
|
|||
|
|
|
|||
|
|
# 复制结果到 final 目录，使用标准命名格式
|
|||
|
|
safe_prefix = args.prefix if args.prefix else "rnaseq"
|
|||
|
|
|
|||
|
|
print(f"\n[PROGRESS] Copying results to final directory...")
|
|||
|
|
copied_files = 0
|
|||
|
|
expected_files_per_sample = 5 # cutadapt: 2 clean + 2 FastQC HTML + 1 log; fastp: 2 clean + HTML + JSON + log
|
|||
|
|
|
|||
|
|
for sample_name, file_info in validated_samples.items():
|
|||
|
|
sample_copied = 0
|
|||
|
|
r1_name = file_info["r1"].stem.replace('.fastq', '').replace('.fq', '') # 原始R1文件名（不含扩展名）
|
|||
|
|
r2_name = file_info["r2"].stem.replace('.fastq', '').replace('.fq', '') # 原始R2文件名（不含扩展名）
|
|||
|
|
|
|||
|
|
# 复制 clean reads 到 final（下游分析需要）
|
|||
|
|
for suffix in ["_1.clean.fastq.gz", "_2.clean.fastq.gz"]:
|
|||
|
|
src = analysis_dir / f"{sample_name}{suffix}"
|
|||
|
|
dst = final_dir / f"{safe_prefix}.rnaseq_qc.{sample_name}{suffix}"
|
|||
|
|
if src.exists() and src.stat().st_size > 0:
|
|||
|
|
shutil.copy2(src, dst)
|
|||
|
|
copied_files += 1
|
|||
|
|
sample_copied += 1
|
|||
|
|
else:
|
|||
|
|
print(f" ⚠ Missing/empty: {src.name}", file=sys.stderr)
|
|||
|
|
|
|||
|
|
if args.qc_tool == "fastp":
|
|||
|
|
for src_name, dst_suffix in [
|
|||
|
|
(f"{sample_name}_fastp.html", f"{sample_name}_fastp.html"),
|
|||
|
|
(f"{sample_name}_fastp.json", f"{sample_name}_fastp.json"),
|
|||
|
|
(f"{sample_name}_fastp.log", f"{sample_name}_fastp.log"),
|
|||
|
|
]:
|
|||
|
|
src = analysis_dir / src_name
|
|||
|
|
dst = final_dir / f"{safe_prefix}.rnaseq_qc.{dst_suffix}"
|
|||
|
|
if src.exists() and src.stat().st_size > 0:
|
|||
|
|
shutil.copy2(src, dst)
|
|||
|
|
copied_files += 1
|
|||
|
|
sample_copied += 1
|
|||
|
|
else:
|
|||
|
|
print(f" ⚠ Missing/empty fastp output: {src_name}", file=sys.stderr)
|
|||
|
|
else:
|
|||
|
|
# 复制 FastQC HTML 报告到 final（前端展示需要）
|
|||
|
|
# FastQC使用原始FASTQ文件名生成报告
|
|||
|
|
fastqc_patterns = [
|
|||
|
|
(f"{r1_name}_fastqc.html", f"{sample_name}_1_fastqc.html"),
|
|||
|
|
(f"{r2_name}_fastqc.html", f"{sample_name}_2_fastqc.html")
|
|||
|
|
]
|
|||
|
|
for src_pattern, dst_suffix in fastqc_patterns:
|
|||
|
|
src = analysis_dir / src_pattern
|
|||
|
|
dst = final_dir / f"{safe_prefix}.rnaseq_qc.{dst_suffix}"
|
|||
|
|
if src.exists() and src.stat().st_size > 0:
|
|||
|
|
shutil.copy2(src, dst)
|
|||
|
|
copied_files += 1
|
|||
|
|
sample_copied += 1
|
|||
|
|
else:
|
|||
|
|
print(f" ⚠ Missing/empty FastQC HTML: {src_pattern}", file=sys.stderr)
|
|||
|
|
|
|||
|
|
# 复制 cutadapt 日志到 final（小文件，包含统计信息）
|
|||
|
|
src = analysis_dir / f"{sample_name}_cutadapt.log"
|
|||
|
|
dst = final_dir / f"{safe_prefix}.rnaseq_qc.{sample_name}_cutadapt.log"
|
|||
|
|
if src.exists() and src.stat().st_size > 0:
|
|||
|
|
shutil.copy2(src, dst)
|
|||
|
|
copied_files += 1
|
|||
|
|
sample_copied += 1
|
|||
|
|
else:
|
|||
|
|
print(f" ⚠ Missing/empty: {src.name}", file=sys.stderr)
|
|||
|
|
|
|||
|
|
print(f" {sample_name}: {sample_copied}/{expected_files_per_sample} files copied")
|
|||
|
|
|
|||
|
|
# 验证输出完整性
|
|||
|
|
total_expected = len(validated_samples) * expected_files_per_sample
|
|||
|
|
if copied_files < total_expected:
|
|||
|
|
print(f"\n❌ ERROR: Output validation failed", file=sys.stderr)
|
|||
|
|
print(f" Expected: {total_expected} files", file=sys.stderr)
|
|||
|
|
print(f" Copied: {copied_files} files", file=sys.stderr)
|
|||
|
|
sys.exit(1)
|
|||
|
|
|
|||
|
|
print(f"\n✅ All {copied_files} output files validated and copied")
|
|||
|
|
|
|||
|
|
# 生成 MultiQC 报告
|
|||
|
|
multiqc_name = f"{args.prefix}.rnaseq_qc_multiqc_report.html" if args.prefix else "rnaseq_qc_multiqc_report.html"
|
|||
|
|
multiqc_report = run_multiqc([analysis_dir, final_dir], final_dir, multiqc_name)
|
|||
|
|
if multiqc_report:
|
|||
|
|
print(f"MultiQC report: {multiqc_report}")
|
|||
|
|
|
|||
|
|
# 生成报告
|
|||
|
|
end_time = datetime.now()
|
|||
|
|
generate_report(final_dir, sample_stats, versions, start_time, end_time, args.prefix, qc_tool=args.qc_tool)
|
|||
|
|
|
|||
|
|
# 生成manifest
|
|||
|
|
cmd = ["python", __file__, "-i", str(input_dir), "-c", str(config_file),
|
|||
|
|
"-o", str(output_dir), "-j", args.job_id, "-r", args.run_id, "-p", args.prefix,
|
|||
|
|
"-t", str(thread_plan["cpu"]),
|
|||
|
|
"--sample-parallel-jobs", str(thread_plan["sample_parallel_jobs"]),
|
|||
|
|
"--threads-per-sample", str(thread_plan["threads_per_sample"]),
|
|||
|
|
"--qc-tool", args.qc_tool]
|
|||
|
|
if args.qc_tool == "fastp":
|
|||
|
|
cmd.extend(["--fastp-bin", str(fastp_bin)])
|
|||
|
|
for flag, value in {
|
|||
|
|
"--cut-window-size": args.cut_window_size,
|
|||
|
|
"--cut-mean-quality": args.cut_mean_quality,
|
|||
|
|
"--qualified-quality-phred": args.qualified_quality_phred,
|
|||
|
|
"--unqualified-percent-limit": args.unqualified_percent_limit,
|
|||
|
|
"--n-base-limit": args.n_base_limit,
|
|||
|
|
"--length-required": args.length_required,
|
|||
|
|
}.items():
|
|||
|
|
cmd.extend([flag, str(value)])
|
|||
|
|
if args.detect_adapter_for_pe:
|
|||
|
|
cmd.append("--detect-adapter-for-pe")
|
|||
|
|
else:
|
|||
|
|
cmd.append("--no-detect-adapter-for-pe")
|
|||
|
|
if args.cut_front:
|
|||
|
|
cmd.append("--cut-front")
|
|||
|
|
else:
|
|||
|
|
cmd.append("--no-cut-front")
|
|||
|
|
if args.cut_tail:
|
|||
|
|
cmd.append("--cut-tail")
|
|||
|
|
else:
|
|||
|
|
cmd.append("--no-cut-tail")
|
|||
|
|
generate_manifest(
|
|||
|
|
final_dir, sample_stats, versions, start_time, end_time,
|
|||
|
|
cmd, input_dir, output_dir, args.job_id, args.run_id, args.prefix,
|
|||
|
|
analysis_dir=analysis_dir, # 传入 analysis 目录
|
|||
|
|
cpu=args.cpu,
|
|||
|
|
sample_parallel_jobs=thread_plan["sample_parallel_jobs"],
|
|||
|
|
threads_per_sample=threads_per_sample,
|
|||
|
|
sample_config=config_file,
|
|||
|
|
qc_tool=args.qc_tool,
|
|||
|
|
fastp_bin=str(fastp_bin) if fastp_bin else "",
|
|||
|
|
fastp_params=fastp_params if args.qc_tool == "fastp" else {},
|
|||
|
|
)
|
|||
|
|
|
|||
|
|
total_elapsed = (end_time - start_time).total_seconds()
|
|||
|
|
print(f"\n{'='*60}")
|
|||
|
|
print(f"✅ RNA-seq QC completed successfully")
|
|||
|
|
print(f" Samples processed: {len(validated_samples)}")
|
|||
|
|
print(f" Total runtime: {total_elapsed:.1f}s ({total_elapsed/60:.1f} min)")
|
|||
|
|
print(f" Results: {output_dir}")
|
|||
|
|
print(f"{'='*60}")
|
|||
|
|
|
|||
|
|
|
|||
|
|
if __name__ == "__main__":
|
|||
|
|
main()
|