annotate pafcount.py @ 22:2ddd41a0c2d5 draft

planemo upload for repository https://github.com/usegalaxy-eu/temporary-tools/tree/master/jbrowse2 commit be4f98b07024b59ff5e1ae0d8b467eecb15c7521-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]))