diff --git a/rna_tools/rna_pdb_tools.py b/rna_tools/rna_pdb_tools.py index a0994a0a..eeb38e0b 100755 --- a/rna_tools/rna_pdb_tools.py +++ b/rna_tools/rna_pdb_tools.py @@ -133,6 +133,8 @@ def get_parser(): parser.add_argument('--get-seq', help='get seq', action='store_true') + parser.add_argument('--rgyration', help='get seq', action='store_true') + parser.add_argument('--color-seq', help='color seq, works with --get-seq', action='store_true') parser.add_argument('--ignore-files', help='files to be ingored, .e.g, \'solution\'') @@ -1194,6 +1196,94 @@ def get_parser(): else: print(txt) + if args.rgyration: + if list != type(args.file): + args.file = [args.file] + # quick fix - make a list on the spot + if list != type(args.file): + args.file = [args.file] + ################################## + analyzed = [] + for f in args.file: + ##################################### + if args.uniq: + subname = eval('f' + args.uniq) + if subname in analyzed: + continue + else: + analyzed.append(subname) + + s = RNAStructure(f) + if not s.is_pdb(): + print('Error: Not a PDB file %s' % f) + sys.exit(1) + s.decap_gtp() + s.std_resn() + s.remove_hydrogen() + s.remove_ion() + s.remove_water() + if args.renum_atoms: + s.renum_atoms() + s.fix_O_in_UC() + s.fix_op_atoms() + + output = '' + + # with # is easier to grep this out + if args.fasta: + # s.fn vs s.name + output += s.get_seq(compact=args.compact, chainfirst=args.chain_first, fasta=args.fasta, addfn=s.name, color=args.color_seq) + '\n' + elif args.oneline: + output += s.get_seq(compact=args.compact, chainfirst=args.chain_first, color=args.color_seq).strip() + ' # '+ os.path.basename(f.replace('.pdb', '')) + '\n' + else: + N = len(s.get_seq(compact=args.compact, chainfirst=args.chain_first, color=args.color_seq)) + output += '' + os.path.basename(f.replace('.pdb', '')).ljust(60) + ' seq len: ' + str(N) + ' ' + try: + sys.stdout.write(output) + sys.stdout.flush() + except IOError: + pass + + import math + def calc_rg(): + """modified after https://github.com/sarisabban/Rg""" + coord = list() + mass = list() + for line in s.lines: + try: + line = line.split() + x = float(line[6]) + y = float(line[7]) + z = float(line[8]) + coord.append([x, y, z]) + atom = line[2][0] + if atom == 'C': + mass.append(12.0107) + elif atom == 'O': + mass.append(15.9994) + elif atom == 'N': + mass.append(14.0067) + elif atom == 'S': + mass.append(32.065) + elif atom == 'P': + mass.append(30.97) + except: + pass + + xm = [(m*i, m*j, m*k) for (i, j, k), m in zip(coord, mass)] + tmass = sum(mass) + rr = sum(mi*i + mj*j + mk*k for (i, j, k), (mi, mj, mk) in zip(coord, xm)) + mm = sum((sum(i) / tmass)**2 for i in zip(*xm)) + rg = math.sqrt(rr / tmass-mm) + return(round(rg, 3)) + + def radius_of_gyration(N): + return round(5.5 * (N ** (1/3)), 2) + + rg = calc_rg() + rgn = radius_of_gyration(N) + print(f'radius of gyration: {rg}; expected for this seqlen ({N}): {rgn}') + if args.renum_nmr: if list != type(args.file): args.file = [args.file]