让程序提速 120 倍

1
def call_peaks_unified_redux(...):
2
    ...
3
    total = len(probes)
4
    for i in range(total):
5
        chrom, start, end, score = probes[i]
6
        if not score:
7
            continue
8
    for chrom, start, end, score in filter(lambda x: x[3], probes):
9
    ...

1
def load_gff(fn: str) -> list[PROBE]:
2
    global args, RAW_READS_NUM
3
    line_num = 0
4
    total_coverage = 0
5
    parsed_result: list[PROBE] = list()
6
    sys.stderr.write(f"Reading input file: {fn} ...\n")
7
    with open(fn, "r") as f:
8
        for line in f:
已折叠 14 行
9
            line_num += 1
10
            if line_num % 10000 == 0:
11
                sys.stderr.write(f"Read {line_num} lines\r")
12
            ll = line.strip().split("\t")
13
            # skip empty lines
14
            if len(ll) < 4:
15
                continue
16
            if len(ll) == 4:
17
                # bedgraph
18
                chrom, start, end, score = ll
19
            else:
20
                # GFF
21
                chrom = ll[0]
22
                start, end, score = ll[3:6]
23
            # increment raw reads number
24
            RAW_READS_NUM += 1
25
            # skip empty reads
26
            if not args.no_discard_zeros and (score == "NA" or not float(score)):
27
                continue
28
            # record read
29
            parsed_result.append(
30
                (
31
                    CHROM(chrom),
32
                    START(POS(int(start))),
33
                    END(POS(int(end))),
34
                    SCORE(float(score) if score != "NA" else 0),
35
                )
36
            )
37
            # record total coverage
38
            total_coverage += int(end) - int(start)
39
    sys.stderr.write(f"Read {line_num} lines\n")
40
    sys.stderr.write("Sorting ...\n")
41
    parsed_result = sorted(parsed_result, key=lambda x: (x[0], x[1]))
42
    sys.stderr.write(f"Total coverage was {total_coverage} bp\n")
43
    return parsed_result

1
def find_quant(probes: list[PROBE]):
2
    global args
3
    total_coverage = 0
4
    frags: list[SCORE] = list()
5
    for (chrom, start, end, score) in filter(lambda x: x[3], probes):
6
        total_coverage += end - start
7
        frags.append(score)
8
    frags = [x[3] for x in probes]
9
    frags = sorted(frags)
10
    sys.stderr.write(f"Total coverage was {total_coverage} bp\n")
已折叠 8 行
11
    quants = [
12
        (q * args.step, int(q * args.step * len(frags)) - 1)
13
        for q in range(math.ceil(args.min_quant / args.step), math.ceil(1 / args.step))
14
    ]
15
    for cut_off, score_idx in quants:
16
        sys.stdout.write(f"\tQuantile {cut_off:0.2f}: {frags[score_idx]:0.2f}\n")
17
    peakmins = [THRESH(frags[score_idx]) for (_, score_idx) in quants]
18
    return peakmins

1
def call_peaks_unified_redux(...):
2
    ...
3
    for pm in peakmins:
4
        ...
5
        for chrom, start, end, score in filter(lambda x: x[3], probes):
6
        for chrom, start, end, score in probes:
7
        ...

1
def find_randomised_peaks(probes: list[PROBE], peakmins: list[THRESH]):
已折叠 5 行
2
    global args, RAW_READS_NUM
3
    peak_count = None
4
    sys.stdout.write("Duplicating ...\n")
5
    pbs = probes.copy()
6
    sys.stdout.write("Calling peaks on input file ...\n")
7
    for iter_num in range(args.n):
8
        sys.stdout.write(f"Iteration {iter_num+1}: [shuffling ...]      \r")
9
        # This is a naive approximation to randomly sample a fraction
10
        # as the full sequence doesn't contain empty reads
11
        # (but no worse than the original approach anyway)
12
        if args.frac:
13
            pbs = pbs[: args.frac]
14
            num_to_sample = sum(
15
                map(
16
                    # This makes sure that it works for both
17
                    # data including and excluding empty reads
18
                    lambda x: x <= int(len(probes) * args.frac),
19
                    random.sample(range(RAW_READS_NUM), int(RAW_READS_NUM * args.frac)),
20
                )
21
            )
已折叠 7 行
22
            pbs = random.sample(probes, num_to_sample)
23
        # The built-in shuffle uses the same algorithm (Fisher-Yates)
24
        # as the original Perl program
25
        random.shuffle(pbs)
26
        _, peak_count = call_peaks_unified_redux(
27
            iter_num, pbs, peakmins, None, peak_count
28
        )
29
    return peak_count