forked from admin/PlantAI-rag-new
chore(bioflow): 删除副本中 13 个 dev 历史备份文件(.bak/.orig)——副本全树零硬编码
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
This commit is contained in:
Claude Fable 5
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
@@ -1,617 +0,0 @@
|
||||
#!/home/xia-lab/miniconda3/envs/test4py311/bin/python
|
||||
"""
|
||||
RNA-Seq QC阶段：对原始测序数据进行质量控制
|
||||
|
||||
输入：
|
||||
- 原始FASTQ文件（双端）
|
||||
- 样本配置文件（group_sample.yaml）
|
||||
|
||||
输出：
|
||||
- FastQC报告
|
||||
- Cutadapt处理后的clean reads
|
||||
- 质控报告和manifest
|
||||
"""
|
||||
|
||||
import argparse
|
||||
import json
|
||||
import subprocess
|
||||
import sys
|
||||
import time
|
||||
import yaml
|
||||
from datetime import datetime
|
||||
from pathlib import Path
|
||||
from typing import Any, Dict, List, Optional
|
||||
import hashlib
|
||||
import os
|
||||
|
||||
|
||||
# 添加conda环境bin目录到PATHCONDA_BIN = "/home/xia-lab/miniconda3/envs/test4py311/bin"if CONDA_BIN not in os.environ.get("PATH", ""): os.environ["PATH"] = CONDA_BIN + ":" + os.environ.get("PATH", "")
|
||||
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 get_tool_versions() -> 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"
|
||||
|
||||
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 = "/home/xia-lab/miniconda3/envs/rnaseq/bin"
|
||||
|
||||
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
|
||||
) -> 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),
|
||||
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 get_sample_stats(sample_dir: Path, sample_name: str) -> Dict[str, Any]:
|
||||
"""获取样本统计信息"""
|
||||
import re
|
||||
stats = {
|
||||
"raw_reads": 0,
|
||||
"clean_reads": 0,
|
||||
"adapter_removed": 0,
|
||||
"quality_filtered": 0
|
||||
}
|
||||
|
||||
# 读取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:
|
||||
match = re.search(r'([\d,]+)', line)
|
||||
if match:
|
||||
stats["raw_reads"] = int(match.group(1).replace(',', ''))
|
||||
elif "Reads with adapters:" in line:
|
||||
match = re.search(r'([\d,]+)', line)
|
||||
if match:
|
||||
stats["adapter_removed"] = int(match.group(1).replace(',', ''))
|
||||
elif "Reads written (passing filters):" in line:
|
||||
match = re.search(r'([\d,]+)', line)
|
||||
if match:
|
||||
stats["clean_reads"] = int(match.group(1).replace(',', ''))
|
||||
|
||||
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 = ""
|
||||
) -> 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
|
||||
|
||||
## 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 目录路径
|
||||
) -> 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",
|
||||
"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": 4,
|
||||
"fastqc_params": "-t 4",
|
||||
"cutadapt_params": "-q 20 -m 50"
|
||||
},
|
||||
"artifacts": {},
|
||||
"tools": versions
|
||||
}
|
||||
|
||||
# 添加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"
|
||||
|
||||
for sample in sample_stats.keys():
|
||||
# FastQC HTML 报告（复制到 final，前端展示）
|
||||
html_r1 = final_dir / f"{prefix}.rnaseq_qc.{sample}_1_fastqc.html"
|
||||
html_r2 = final_dir / f"{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）- 大文件
|
||||
clean_r1 = output_dir / f"{sample}_1.clean.fastq.gz"
|
||||
clean_r2 = output_dir / f"{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"{sample}_cutadapt.log"
|
||||
add_artifact(f"cutadapt_log_{sample}", log_file, "log", format="txt")
|
||||
|
||||
# 报告和 manifest
|
||||
report_file = output_dir / report_name
|
||||
add_artifact("report", report_file, "markdown", format="md")
|
||||
|
||||
# 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="Number of CPU threads")
|
||||
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")
|
||||
|
||||
args = parser.parse_args()
|
||||
|
||||
start_time = datetime.now()
|
||||
input_dir = Path(args.input)
|
||||
output_dir = Path(args.output) # rnaseq_qc 工作目录
|
||||
config_file = Path(args.config)
|
||||
|
||||
# 创建工作目录
|
||||
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...")
|
||||
|
||||
# 获取工具版本
|
||||
versions = get_tool_versions()
|
||||
print(f"Tool versions: {versions}")
|
||||
|
||||
sample_stats = {}
|
||||
total_samples = len(validated_samples)
|
||||
current_sample = 0
|
||||
|
||||
# 处理每个样本
|
||||
for sample_name, file_info in validated_samples.items():
|
||||
current_sample += 1
|
||||
r1 = file_info["r1"]
|
||||
r2 = file_info["r2"]
|
||||
|
||||
print(f"\n[PROGRESS {current_sample}/{total_samples}] Processing: {sample_name}", flush=True)
|
||||
print(f" Input: {r1.name} + {r2.name}", flush=True)
|
||||
|
||||
# 运行FastQC - 输出到 analysis 目录
|
||||
try:
|
||||
fastqc_result = run_fastqc([r1, r2], analysis_dir, args.cpu)
|
||||
except RuntimeError as e:
|
||||
print(f"\n❌ ERROR: FastQC failed for {sample_name}: {e}", file=sys.stderr)
|
||||
sys.exit(1)
|
||||
|
||||
# 运行Cutadapt - 输出到 analysis 目录
|
||||
out_r1 = analysis_dir / f"{sample_name}_1.clean.fastq.gz"
|
||||
out_r2 = analysis_dir / f"{sample_name}_2.clean.fastq.gz"
|
||||
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)
|
||||
except RuntimeError as e:
|
||||
print(f"\n❌ ERROR: Cutadapt failed for {sample_name}: {e}", file=sys.stderr)
|
||||
sys.exit(1)
|
||||
|
||||
# 获取统计信息（从 analysis 目录读取）
|
||||
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)
|
||||
sample_stats[sample_name] = stats
|
||||
|
||||
print(f" ✓ Sample {sample_name} completed", flush=True)
|
||||
|
||||
# 复制结果到 final 目录，使用标准命名格式
|
||||
import shutil
|
||||
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 # 2 clean FASTQ + 2 FastQC HTML + 1 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)
|
||||
|
||||
# 复制 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
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
Reference in New Issue
Block a user