75 lines
2.1 KiB
Python
75 lines
2.1 KiB
Python
"""Get summary statistics from GFF file."""
|
|
from collections import Counter
|
|
from pathlib import Path
|
|
import pickle
|
|
|
|
import gffutils
|
|
from .util import wf_parser # noqa: ABS101
|
|
|
|
|
|
def argparser():
|
|
"""Argument parser for entrypoint."""
|
|
parser = wf_parser("summ_gff")
|
|
parser.add_argument(
|
|
"gff",
|
|
help="Report output file",
|
|
type=Path)
|
|
parser.add_argument(
|
|
"sample_id",
|
|
help="Output TSV file path")
|
|
parser.add_argument(
|
|
"out",
|
|
default="gff_summary.tsv",
|
|
help="Output TSV file path",
|
|
type=Path)
|
|
|
|
return parser
|
|
|
|
|
|
def main(args):
|
|
"""Entry point."""
|
|
db = gffutils.create_db(
|
|
str(args.gff), dbfn=':memory:', force=True, keep_order=True,
|
|
merge_strategy='merge', sort_attribute_values=True
|
|
)
|
|
|
|
num_transcripts = db.count_features_of_type('transcript')
|
|
num_genes = db.count_features_of_type('gene')
|
|
|
|
transcript_lens = []
|
|
exons_per_transcript = Counter()
|
|
isoforms_per_gene = Counter()
|
|
|
|
for gene in db.features_of_type('gene'):
|
|
|
|
n_isos = len(list(db.children(gene, featuretype='transcript')))
|
|
isoforms_per_gene[n_isos] += 1
|
|
|
|
for transcript in db.children(
|
|
gene, featuretype='transcript', order_by='start'):
|
|
tr_len = 0
|
|
exons = list(db.children(transcript, featuretype='exon'))
|
|
if len(exons) == 0:
|
|
continue
|
|
exons_per_transcript[len(exons)] += 1
|
|
for ex in exons:
|
|
tr_len += abs(ex.end - ex.start)
|
|
|
|
transcript_lens.append(tr_len)
|
|
|
|
results = {
|
|
'sample_id': args.sample_id,
|
|
'summaries': {
|
|
'Total genes': [num_genes],
|
|
'Total transcripts': [num_transcripts],
|
|
'Max trans. len': max(transcript_lens),
|
|
'Min trans. len': min(transcript_lens)
|
|
},
|
|
'transcript_lengths': transcript_lens,
|
|
'exons_per_transcript': exons_per_transcript,
|
|
'isoforms_per_gene': isoforms_per_gene
|
|
}
|
|
|
|
with open(args.out, 'wb') as fh:
|
|
pickle.dump(results, fh)
|