-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathmain.py
More file actions
78 lines (65 loc) · 3.11 KB
/
Copy pathmain.py
File metadata and controls
78 lines (65 loc) · 3.11 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
import argparse
import os
import time
import pandas as pd
from circnote.read_data.input_validation import validate_junction, validate_gtf
from circnote.read_data.junction_parser import parse_junction
from circnote.read_data.gtf_parser import load_exons
from circnote.detection.bsj_detector import find_bsj_candidates
def main():
parser = argparse.ArgumentParser(description="circnote: circRNA BSJ detection and annotation")
parser.add_argument("--junction", required=True, help="Path to STAR Chimeric.out.junction file")
parser.add_argument("--gtf", required=True, help="Path to GTF annotation file")
parser.add_argument("--chromosome", nargs="*", default=None, help="Chromosome(s) to process, e.g. chr1 chr14 (default: all)")
parser.add_argument("--min_reads", type=int, default=2, help="Minimum read support (default: 2)")
parser.add_argument("--max_offset", type=int, default=0, help="Maximum offset allowed between junction and exon boundary (default: 0)")
parser.add_argument("--output", required=True, help="Output TSV file path")
args = parser.parse_args()
validate_junction(args.junction)
validate_gtf(args.gtf)
chromosomes = set(args.chromosome) if args.chromosome else None
t0 = time.time()
junction_df = parse_junction(
args.junction,
chromosomes=chromosomes,
min_reads=args.min_reads
)
print(f"Junction parsing: {time.time() - t0:.2f}s")
# caching frequently used data for optimization
cache_dir = "cache"
os.makedirs(cache_dir, exist_ok=True)
cache_filename = os.path.basename(args.gtf) + ".exons_all.pkl"
cache_path = os.path.join(cache_dir, cache_filename)
t1 = time.time()
if os.path.exists(cache_path) and os.path.getmtime(cache_path) > os.path.getmtime(args.gtf):
print(f"Loading cached exon annotations from {cache_path}")
exon_df = pd.read_pickle(cache_path)
print(f"Cache load: {time.time() - t1:.2f}s")
else:
print("Parsing GTF exon annotations...")
exon_df = load_exons(
args.gtf,
chromosomes=None
)
exon_df.to_pickle(cache_path)
print(f"GTF parse + cache write: {time.time() - t1:.2f}s")
if chromosomes is not None:
exon_df = exon_df[exon_df["chromosome"].isin(chromosomes)].copy()
t2 = time.time()
summary_df, transcript_df = find_bsj_candidates(
junction_df,
exon_df,
max_offset=args.max_offset
)
print(f"BSJ detection: {time.time() - t2:.2f}s")
base, ext = (args.output.rsplit(".", 1) if "." in args.output else (args.output, ""))
transcript_output = f"{base}.transcripts.{ext}" if ext else f"{base}.transcripts"
summary_df.to_csv(args.output, sep="\t", index=False)
transcript_df.to_csv(transcript_output, sep="\t", index=False)
print(f"Done. {len(summary_df)} BSJ candidates written to {args.output}")
print(f" {len(transcript_df)} transcript-level rows written to {transcript_output}")
if __name__ == "__main__":
start = time.time()
main()
end = time.time()
print(f"\nTotal runtime: {end - start:.2f} seconds")