diff galaxy-tools/biobank/tools/build_miniped.py @ 3:43be74e62bfe draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/galaxy-tools/biobank/tools/build_miniped.py	Thu Sep 22 08:57:04 2016 -0400
@@ -0,0 +1,155 @@
+# BEGIN_COPYRIGHT
+# END_COPYRIGHT
+
+"""
+A rough example of basic pedigree info generation.
+"""
+
+import csv, argparse, sys, os
+
+from bl.vl.kb import KnowledgeBase as KB
+from bl.vl.kb.drivers.omero.ehr import EHR
+import bl.vl.individual.pedigree as ped
+import bl.vl.utils.ome_utils as vlu
+from bl.vl.utils import LOG_LEVELS, get_logger
+
+DIAGNOSIS_ARCH = 'openEHR-EHR-EVALUATION.problem-diagnosis.v1'
+DIAGNOSIS_FIELD = 'at0002.1'
+T1D_ICD10 = 'icd10-cm:E10'
+MS_ICD10 = 'icd10-cm:G35'
+NEFRO_ICD10 = 'icd10-cm:E23.2'
+
+PLINK_MISSING = -9
+PLINK_UNAFFECTED = 1
+PLINK_AFFECTED = 2
+
+FIELDS = ["fam_label", "ind_label", "fat_label", "mot_label", "gender", "t1d_status", "ms_status", "nefro_status"]
+
+
+def make_parser():
+  parser = argparse.ArgumentParser(description='build the first columns of a ped file from VL')
+  parser.add_argument('--logfile', type=str, help='log file (default=stderr)')
+  parser.add_argument('--loglevel', type=str, choices = LOG_LEVELS,
+                      help='logging level', default='INFO')
+  parser.add_argument('-H', '--host', type=str, help='omero hostname')
+  parser.add_argument('-U', '--user', type=str, help='omero user')
+  parser.add_argument('-P', '--passwd', type=str, help='omero password')
+  parser.add_argument('-S', '--study', type=str, required=True,
+                      help='a list of comma separated studies used to retrieve individuals that will be written to ped file')
+  parser.add_argument('--ofile', type=str, help='output file path',
+                      required=True)
+  return parser
+
+def build_families(individuals, logger):
+  # Individuals with only one parent will be considered like founders
+  # for i in individuals:
+  #   if ((i.mother is None) or (i.father is None)):
+  #     i.mother = None
+  #     i.father = None
+  logger.info("individuals: %d" % len(individuals))
+  #logger.info("individuals: with 0 or 2 parents: %d" % len(not_one_parent))
+  logger.info("analyzing pedigree")
+  founders, non_founders, dangling, couples, children = ped.analyze(
+    individuals
+    )
+  logger.info("splitting into families")
+  return ped.split_disjoint(individuals, children)
+
+
+def main(argv):
+  parser = make_parser()
+  args = parser.parse_args(argv)
+
+  logger = get_logger('build_miniped', level=args.loglevel,
+                      filename=args.logfile)
+  
+  try:
+    host = args.host or vlu.ome_host()
+    user = args.user or vlu.ome_user()
+    passwd = args.passwd or vlu.ome_passwd()
+  except ValueError, ve:
+    logger.critical(ve)
+    sys.exit(ve)
+
+  kb = KB(driver='omero')(host, user, passwd)
+  logger.debug('Loading all individuals from omero')
+  all_inds = kb.get_objects(kb.Individual)  # store all inds to cache
+  logger.debug('%d individuals loaded' % len(all_inds))
+  studies = [kb.get_study(s) for s in args.study.split(',')]
+  # Removing None values
+  studies = set(studies)
+  try:
+    studies.remove(None)
+  except KeyError:
+    pass
+  studies = list(studies)
+  if len(studies) == 0:
+    logger.error('No matches found for labels %s, stopping program' % args.study)
+    sys.exit(2)
+  enrolled_map = {}
+  for study in studies:
+    logger.info('Loading enrolled individuals for study %s' % study.label)
+    enrolled = kb.get_enrolled(study)
+    logger.debug('%d individuals loaded' % len(enrolled))
+    for en in enrolled:
+      if en.individual.id not in enrolled_map:
+        enrolled_map[en.individual.id] = ('%s:%s' % (en.study.label, en.studyCode),
+                                          en.individual)
+      else:
+        logger.debug('Individual %s already mapped' % en.individual.id)
+    logger.debug('Loading EHR records')
+    ehr_records = kb.get_ehr_records()
+    logger.debug('%s EHR records loaded' % len(ehr_records))
+    ehr_records_map = {}
+    for r in ehr_records:
+      ehr_records_map.setdefault(r['i_id'], []).append(r)
+    affection_map = {}
+    for ind_id, ehr_recs in ehr_records_map.iteritems():
+      affection_map[ind_id] = dict(t1d=PLINK_UNAFFECTED, ms=PLINK_UNAFFECTED,
+                                   nefro=PLINK_UNAFFECTED)
+      ehr = EHR(ehr_recs)
+      if ehr.matches(DIAGNOSIS_ARCH, DIAGNOSIS_FIELD, T1D_ICD10):
+        affection_map[ind_id]['t1d'] = PLINK_AFFECTED
+      if ehr.matches(DIAGNOSIS_ARCH, DIAGNOSIS_FIELD, MS_ICD10):
+        affection_map[ind_id]['ms'] = PLINK_AFFECTED
+      if ehr.matches(DIAGNOSIS_ARCH, DIAGNOSIS_FIELD, NEFRO_ICD10):
+        affection_map[ind_id]['nefro'] = PLINK_AFFECTED
+
+  immuno_inds = [i for (ind_id, (st_code, i)) in enrolled_map.iteritems()]
+  families = build_families(immuno_inds, logger)
+  logger.info("found %d families" % len(families))
+  
+  def resolve_label(i):
+    try:
+      return enrolled_map[i.id][0]
+    except KeyError:
+      return i.id
+
+  def resolve_pheno(i):
+    try:
+      immuno_affection = affection_map[i.id]
+    except KeyError:
+      return PLINK_MISSING, PLINK_MISSING, PLINK_MISSING
+    return immuno_affection["t1d"], immuno_affection["ms"], immuno_affection["nefro"]
+
+  kb.Gender.map_enums_values(kb)
+  gender_map = lambda x: 2 if x == kb.Gender.FEMALE else 1
+
+  logger.info("writing miniped")
+  with open(args.ofile, "w") as f:
+    writer = csv.DictWriter(f, FIELDS, delimiter="\t", lineterminator="\n")
+    for k, fam in enumerate(families):
+      fam_label = "FAM_%d" % (k+1)
+      for i in fam:
+        r = {}
+        r["fam_label"] = fam_label
+        r["ind_label"] = resolve_label(i)
+        r["fat_label"] = 0 if (i.father is None or i.father not in fam) else resolve_label(i.father)
+        r["mot_label"] = 0 if (i.mother is None or i.mother not in fam) else resolve_label(i.mother)
+        r["gender"] = gender_map(i.gender)
+        r["t1d_status"], r["ms_status"], r["nefro_status"] = resolve_pheno(i)
+        writer.writerow(r)
+
+
+if __name__ == "__main__":
+  main(sys.argv[1:])