86 lines
2.9 KiB
Python
Executable File
86 lines
2.9 KiB
Python
Executable File
from __future__ import annotations
|
|
|
|
import importlib.util
|
|
import sys
|
|
from pathlib import Path
|
|
|
|
|
|
SCRIPT = Path(__file__).resolve().parents[2] / "scripts" / "check_rnaseq_strandness.py"
|
|
|
|
|
|
def load_module():
|
|
assert SCRIPT.exists(), f"Missing script: {SCRIPT}"
|
|
spec = importlib.util.spec_from_file_location("check_rnaseq_strandness", SCRIPT)
|
|
assert spec and spec.loader
|
|
module = importlib.util.module_from_spec(spec)
|
|
sys.modules[spec.name] = module
|
|
spec.loader.exec_module(module)
|
|
return module
|
|
|
|
|
|
def test_discovers_paired_fastqs(tmp_path):
|
|
module = load_module()
|
|
for name in ["S1_1.fastq.gz", "S1_2.fastq.gz", "S2_R1.fq.gz", "S2_R2.fq.gz"]:
|
|
(tmp_path / name).write_text("")
|
|
|
|
pairs = module.discover_fastq_pairs(tmp_path)
|
|
|
|
assert [p.sample for p in pairs] == ["S1", "S2"]
|
|
assert pairs[0].r1.name == "S1_1.fastq.gz"
|
|
assert pairs[0].r2.name == "S1_2.fastq.gz"
|
|
assert pairs[1].r1.name == "S2_R1.fq.gz"
|
|
assert pairs[1].r2.name == "S2_R2.fq.gz"
|
|
|
|
|
|
def test_gtf_conversion_writes_bed12(tmp_path):
|
|
module = load_module()
|
|
gtf = tmp_path / "genes.gtf"
|
|
bed = tmp_path / "genes.bed12"
|
|
gtf.write_text(
|
|
"\n".join(
|
|
[
|
|
'chr1\ttest\ttranscript\t101\t300\t.\t+\t.\tgene_id "g1"; transcript_id "t1";',
|
|
'chr1\ttest\texon\t101\t150\t.\t+\t.\tgene_id "g1"; transcript_id "t1";',
|
|
'chr1\ttest\texon\t251\t300\t.\t+\t.\tgene_id "g1"; transcript_id "t1";',
|
|
]
|
|
)
|
|
+ "\n"
|
|
)
|
|
|
|
count = module.convert_annotation_to_bed12(gtf, bed)
|
|
|
|
assert count == 1
|
|
fields = bed.read_text().strip().split("\t")
|
|
assert fields[:6] == ["chr1", "100", "300", "t1", "0", "+"]
|
|
assert fields[9:] == ["2", "50,50", "0,150"]
|
|
|
|
|
|
def test_internal_orientation_classifies_reverse_from_sam_lines(tmp_path):
|
|
module = load_module()
|
|
bed = tmp_path / "genes.bed12"
|
|
bed.write_text("chr1\t100\t500\tt1\t0\t+\t100\t500\t0\t1\t400\t0\n")
|
|
sam_lines = [
|
|
"r1\t83\tchr1\t120\t60\t50M\t=\t220\t150\t*\t*",
|
|
"r1\t163\tchr1\t220\t60\t50M\t=\t120\t-150\t*\t*",
|
|
"r2\t83\tchr1\t140\t60\t50M\t=\t240\t150\t*\t*",
|
|
"r2\t163\tchr1\t240\t60\t50M\t=\t140\t-150\t*\t*",
|
|
]
|
|
|
|
result = module.infer_orientation_from_sam_lines(sam_lines, bed, threshold=0.8)
|
|
|
|
assert result["classification"] == "reverse"
|
|
assert result["hisat2_rna_strandness"] == "RF"
|
|
assert result["featurecounts_s"] == "2"
|
|
assert result["informative_alignments"] == 4
|
|
|
|
|
|
def test_simulation_dataset_contains_required_files(tmp_path):
|
|
module = load_module()
|
|
dataset = module.create_simulation_dataset(tmp_path, library="reverse", read_count=4)
|
|
|
|
assert dataset["genome"].exists()
|
|
assert dataset["annotation"].exists()
|
|
assert dataset["fastq_dir"].exists()
|
|
assert (dataset["fastq_dir"] / "SIM_reverse_1.fastq.gz").exists()
|
|
assert (dataset["fastq_dir"] / "SIM_reverse_2.fastq.gz").exists()
|