view pafcount.py @ 24:fb6cc7bc24df draft

planemo upload for repository https://github.com/usegalaxy-eu/temporary-tools/tree/master/jbrowse2 commit cde4b8a574ded34a0ff8df3ecafc1a057787dcfb-dirty

import sys

"""

Col     Type    Description
1   string  Query sequence name
2   int     Query sequence length
3   int     Query start (0-based; BED-like; closed)
4   int     Query end (0-based; BED-like; open)
5   char    Relative strand: "+" or "-"
6   string  Target sequence name
7   int     Target sequence length
8   int     Target start on original strand (0-based)
9   int     Target end on original strand (0-based)
10  int     Number of residue matches
11  int     Alignment block length
12  int     Mapping quality (0-255; 255 for missing)
"""

qcis = {}
tcis = {}
qtrans = {}
ttrans = {}
pafname = sys.argv[1]
pf = open(pafname, "r").readlines()
for row in pf:
    qn, ql, qs, qe, qrs, tn, tl, ts, te, nm, abl, mq = row.strip().split("\t")[:12]

    if qn == tn:  # cis
        print("cis", qn, tn)
        tcis.setdefault(tn, 0)
        tcis[tn] = tcis[tn] + 1
        qcis.setdefault(qn, 0)
        qcis[qn] = qcis[qn] + 1
    else:  # trans
        print("trans", qn, tn)
        k = "%s_%s" % (qn, tn)
        ttrans.setdefault(k, 0)
        ttrans[k] = ttrans[k] + 1
        qtrans.setdefault(k, 0)
        qtrans[k] = qtrans[k] + 1
# print('qcis', qcis,'\nqtrans', qtrans,'\ntcis', tcis,'\ntt', ttrans)
# print('\nqtrans', qtrans,'\nttrans', ttrans)
chroms = list(qtrans.keys())
print("chroms=", chroms)
# print('chrom\tqcis\ttcis\tqtrans\tttrans')
print("chrom\tqtrans\tttrans")
for cr in chroms:
    # print('%s\t%d\t%d\t%d\t%d' % (cr, qcis[cr], tcis[cr], qtrans[cr], ttrans[cr]))
    print("%s\t%d\t%d" % (cr, qtrans[cr], ttrans[cr]))