code: def match_subfeats(f1, f2, dbx1, dbx2, featuretype=None, slop=False):
      
      f1c, f2c = list(dbx1.children(f1, featuretype=featuretype, order_by='start')), \
              list(dbx2.children(f2, featuretype=featuretype, order_by='start'))
  
      lf1c, lf2c = len(f1c), len(f2c)
      if match_nchildren(f1c, f2c):
          if lf1c > 0 and lf2c > 0:
              exclN = set()
              if featuretype.endswith('UTR') or featuretype == 'exon':
                  N = []
                  if featuretype.startswith('five_prime'):
                      N = [1] if f1.strand == "+" else [lf1c]
                  elif featuretype.startswith('three_prime'):
                      N = [lf1c] if f1.strand == "+" else [1]
                  else:   # infer UTR from exon collection
                      N = [1] if 1 == lf1c else [1, lf1c]
  
                  for n in N:
                      if match_Nth_child(f1c, f2c, N=n, slop=slop):
                          exclN.add(n-1)
                      else:
                          return False
  
              for i, (cf1, cf2) in enumerate(zip(f1c, f2c)):
                  if i in exclN: continue
                  if not match_span(cf1, cf2):
                      return False
      else:
          if (lf1c, lf2c) in [(0, 1), (1, 0)] and slop \
                  and featuretype.endswith('UTR'):
              return True
  
          return False
  
      return True
docs: Given 2 gffutils features located in 2 separate gffutils databases,
      iterate through all subfeatures of a certain type and check whether
      they are identical or not
  
      The `slop` parameter allows for variation in the terminal UTR region