ONT_ARTIC_Summaries/simple_stats.py at main · LooseLab/ONT_ARTIC_Summaries · GitHub

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
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
import os, sys
import numpy as np
import pandas as pd
import argparse

if __name__ == "__main__":
    p = argparse.ArgumentParser()
    p.add_argument(
        "-o",
        "--output",
        help="File name to write out the summary stats.",
        metavar="FILE",
        required=True,
    )
    p.add_argument(
        "input",
        help="One or many input files. These must be ONT sequencing summary files as called by MinKNOW or Guppy with bardoing options set. Reads demultiplexed with PoreChop will not work. The files can be optionally compressed.",
        nargs="+",
        metavar="FILE",
    )
    p.add_argument(
        "-c",
        "--centre",
        help="Sequencing Center Generating this Data.",
        metavar="STRING",
        required=True,
    )
    args = p.parse_args()
    # Flowcell_id	run_id	experiment_id	sample_id	pore_count	run_time	number_of_barcodes	barcode_id	pass_filtering	read_count	yield	mean_length	median_length	std_length	min_length	max_length

    dtypes = {
        "filename_fastq": "str",
        "mux": "int64",
        "channel": "int64",
        "start_time": "float64",
        "duration": "float64",
        "run_id": "category",
        "experiment_id": "category",
        "sample_id": "category",
        "passes_filtering": "category",
        "barcode_arrangement": "category",
        "sequence_length_template": "int64",
    }

    dtypes2 = {
        "filename": "str",
        "mux": "int64",
        "channel": "int64",
        "start_time": "float64",
        "duration": "float64",
        "run_id": "category",
        "passes_filtering": "category",
        "barcode_arrangement": "category",
        "sequence_length_template": "int64",
    }

    for file_to_read in args.input:

        try:
            df2 = pd.read_csv(
                file_to_read,
                sep="\t",
                usecols=dtypes.keys(),
            )
            df2 = df2.astype(dtypes)
        except ValueError as e:
            print("{!r} could not be read.".format(file_to_read), file=sys.stderr)
            print("{!r}.".format(e), file=sys.stderr)
            print("Trying alternate approach.", file=sys.stderr)
            try:
                df2 = pd.read_csv(
                    file_to_read,
                    sep="\t",
                    usecols=dtypes2.keys(),
                )
                df2 = df2.astype(dtypes2)
                df2["filename_fastq"]=df2["filename"]
                df2["experiment_id"]="NaN"
                df2["sample_id"]="NaN"

            except ValueError as p:
                print("{!r} could not be read.".format(file_to_read), file=sys.stderr)
                print("{!r}.".format(p), file=sys.stderr)
                continue

        df2["flowcell_id"] = df2["filename_fastq"].str[0:8]

        pore_count = {}
        run_time = {}
        barcode_count = {}
        for run_id, gb in df2.groupby(["run_id"]):
            pore_count[run_id] = len(
                gb.groupby(["mux", "channel"])
                .size()
                .reset_index()
                .rename(columns={0: "count"})
            )
            run_time[run_id] = np.max(gb["start_time"] + gb["duration"])
            barcode_count[run_id] = len(gb["barcode_arrangement"].unique())

        df_final = df2.groupby(
            [
                "run_id",
                "experiment_id",
                "sample_id",
                "flowcell_id",
                "passes_filtering",
                "barcode_arrangement",
            ]
        ).agg(
            {
                "sequence_length_template": [
                    "sum",
                    "count",
                    "min",
                    "max",
                    "mean",
                    "median",
                    "std",
                ]
            }
        )

        df_final.columns = df_final.columns.droplevel(0)
        df_final = df_final.reset_index()

        df_final["pore_count"] = df_final["run_id"].map(pore_count)
        df_final["barcode_count"] = df_final["run_id"].map(barcode_count)
        df_final["run_time"] = df_final["run_id"].map(run_time)
        df_final["sequencing_centre"] = args.centre

        df_final = df_final.rename(
            columns={
                "sum": "yield",
                "count": "read_count",
                "min": "min_length",
                "max": "max_length",
                "mean": "mean_length",
                "median": "median_length",
                "std": "std_length",
            }
        )

        output_order = [
            "sequencing_centre",
            "run_id",
            "experiment_id",
            "sample_id",
            "flowcell_id",
            "run_time",
            "pore_count",
            "barcode_count",
            "passes_filtering",
            "barcode_arrangement",
            "yield",
            "read_count",
            "min_length",
            "max_length",
            "mean_length",
            "median_length",
            "std_length",
        ]
        if not os.path.isfile(args.output):
            df_final[output_order].to_csv(
                args.output, sep="\t", header=True, index=False
            )
        else:
            df_final[output_order].to_csv(
                args.output, sep="\t", header=False, index=False, mode="a"
            )