NHS_WES_generate_DEC_IGV.py

        else:
            if (key in DAD_DICT['x-linked dominant']) and (DAD_STAT == "UNAFFECTED"):
                DAD_GT = DAD_DICT['x-linked dominant'][key][0]
                if DAD_GT == 'HOM':                                             # i.e., x-linked dominant variant in unnafected father
                    print "***[X-linked model (x-linked dominant)]*** Excluded CHILD var %s CHILD_GT = %s, DAD_GT = %s, DAD_STAT = %s" % (key,CHILD_GT,DAD_GT,DAD_STAT)
                    continue

        # if a non-normalized INDEL in child G2P - must adjust (should not happen really, we split, normalized and left-aligned the family VCF before sending it to VEP+G2P)
        chr,start,end,ref,alt = key.split(":")
        if len(ref) > 1 and len(alt) > 1:                           # an INDEL - not normalized
            if len(ref) < len(alt):                                 # an INS
                orig_start = start
                orig_ref = ref
                orig_alt = alt
                start = orig_start
                ref = '-'
                alt = orig_alt[len(orig_ref):]
                print "    WARNING: original INS = %s:%s:%s:%s:%s --> replaced with INS = %s:%s:%s:%s" % (chr,orig_start,end,orig_ref,orig_alt,chr,start,ref,alt)
            else:                                                   # a DEL
                print "ERROR: At the momemnt, cannot deal with this non-normalized deletion"
                print line
                raise SystemExit

        new_key = '%s:%s:%s:%s' % (chr,start,ref,alt)

        # record the data for CHILD G2P variants (for OBS=x-linked dominant)
        if new_key not in G2P_DICT:
            G2P_DICT[new_key] = 0
        else:
            # print "ERROR: duplicate G2P variant new_key = %s" % (new_key)
            # raise SystemExit
            # this will happen if a gene is e.g. hemizygous,x-linked dominant - there will be two separate lines in the output for each req
            pass

        # and record the required data (CHILD_TRANS,CHILD_GENE,CHILD_GT) in G2P_DATA
        if new_key not in G2P_DATA:
            G2P_DATA[new_key] = (CHILD_TRANS,CHILD_GENE,CHILD_GT)
        else:
            # print "ERROR: duplicate G2P variant new_key = %s" % (new_key)
            # raise SystemExit
            # this will happen if a gene is e.g. hemizygous,x-linked dominant - there will be two separate lines in the output for each req
            pass


    ########################################################################
    ### process x-linked over-dominance  genes - no filtering to be done ###
    ########################################################################

    for key in CHILD_DICT['x-linked over-dominance']:       # this the second key: chr:start:end:ref:alt; value: (ZYG,gene,trans)

        CHILD_GT = CHILD_DICT['x-linked over-dominance'][key][0]
        CHILD_GENE = CHILD_DICT['x-linked over-dominance'][key][1]
        CHILD_TRANS = CHILD_DICT['x-linked over-dominance'][key][2]

        # if a non-normalized INDEL in child G2P - must adjust (should not happen really, we split, normalized and left-aligned the family VCF before sending it to VEP+G2P)
        chr,start,end,ref,alt = key.split(":")
        if len(ref) > 1 and len(alt) > 1:                           # an INDEL - not normalized
            if len(ref) < len(alt):                                 # an INS
                orig_start = start
                orig_ref = ref
                orig_alt = alt
                start = orig_start
                ref = '-'
                alt = orig_alt[len(orig_ref):]
                print "    WARNING: original INS = %s:%s:%s:%s:%s --> replaced with INS = %s:%s:%s:%s" % (chr,orig_start,end,orig_ref,orig_alt,chr,start,ref,alt)
            else:                                                   # a DEL
                print "ERROR: At the momemnt, cannot deal with this non-normalized deletion"
                print line
                raise SystemExit

        new_key = '%s:%s:%s:%s' % (chr,start,ref,alt)

        # record the data for CHILD G2P variants (for OBS=x-linked over-dominance)
        if new_key not in G2P_DICT:
            G2P_DICT[new_key] = 0
        else:
            # print "ERROR: duplicate G2P variant new_key = %s" % (new_key)
            # raise SystemExit
            # this will happen if a gene is e.g. hemizygous,x-linked dominant - there will be two separate lines in the output for each req
            pass

        # and record the required data (CHILD_TRANS,CHILD_GENE,CHILD_GT) in G2P_DATA
        if new_key not in G2P_DATA:
            G2P_DATA[new_key] = (CHILD_TRANS,CHILD_GENE,CHILD_GT)
        else:
            # print "ERROR: duplicate G2P variant new_key = %s" % (new_key)
            # raise SystemExit
            # this will happen if a gene is e.g. hemizygous,x-linked dominant - there will be two separate lines in the output for each req
            pass


    NUM_UNIQ_G2P_VARS = len(G2P_DICT)
    print "Found %s unique G2P variants in CHILD (%s) after considering MONOALLELIC and X-LINKED genes" % (NUM_UNIQ_G2P_VARS,CHILD_ID)
    sys.stdout.flush()

    print ""





    ##############################################################################################################
    ####    BIALLELIC filtering                                                                               ####
    ####    under the biallelic model (OBS == biallelic) - consider ALL variants per gene                     ####
    ####    must all be HET in CHILD, GT in parent does not matter                                            ####
    ####    all of them must *clearly* come from only one of the parents (maternally/paternally + biparental) ####
    ####    and this parent must be unaffected                                                                ####
    ####    if all these: then exclude all child variants in this gene                                        ####
    ##############################################################################################################


    print ""
    print "===   BIALLELIC filtering   ==="


    GENE_KEY_GT = defaultdict(dict)		# for child - 1st level key: gene_name; 2nd level key: chr:start:end:ref:alt; value: (GT,trans)

    # process all variants in biallelic genes in child
    for key in CHILD_DICT['biallelic']:		# this the second key: chr:start:end:ref:alt; value: (ZYG,gene,trans)
        b_GT = CHILD_DICT['biallelic'][key][0]
        b_gene = CHILD_DICT['biallelic'][key][1]  
        b_trans = CHILD_DICT['biallelic'][key][2]
        GENE_KEY_GT[b_gene][key] = (b_GT,b_trans)

    # iterate over genes in GENE_KEY_GT
    for g in GENE_KEY_GT: 			# this is the biallelic gene name 
        all_HET = True

        # iterate over variants in this gene 
        for kx in GENE_KEY_GT[g]:		# this the second key: chr:start:end:ref:alt
###            if GENE_KEY_GT[g][kx] == 'HOM':	# there is a HOM variant in the child - NO filtering
            if GENE_KEY_GT[g][kx][0] == 'HOM':     # there is a HOM variant in the child - NO filtering
                all_HET = False
                break

        if all_HET:				# for this gene
        # all variants in this gene in the CHILD are HET, check if all come from a single unaffected parent
        # if yes, filter out and write a message to the log file
        # if not, to be added to G2P_DICT and G2P_DATA for further processing

            all_from_one_parent = True

            # iterate again over the variants in this gene
            VAR_SOURCE_LIST = {}		# key: chr:start:end:ref:alt in child; value: (NONE) or (MOM or DAD or BOTH and the parent is UNAFFECTED)

            for ky in GENE_KEY_GT[g]:		# this the second key: chr:start:end:ref:alt

                this_var_status = 'NONE'

                if ((ky in MOM_DICT['biallelic']) or (ky in MOM_DICT['monoallelic'])) and (MOM_STAT == "UNAFFECTED"):		
                    this_var_status = 'MOM'
                if ((ky in DAD_DICT['biallelic']) or (ky in DAD_DICT['monoallelic'])) and (DAD_STAT == "UNAFFECTED"):
                    if this_var_status == 'NONE': 
                        this_var_status = 'DAD'
                    elif this_var_status == 'MOM':
                        this_var_status = 'BOTH'

                VAR_SOURCE_LIST[ky] = this_var_status

            # have collected the parent source for all variants in this gene
            tot_num_vars = len(VAR_SOURCE_LIST)
            num_mom = 0
            num_dad = 0
            num_none = 0
            for kt,v in VAR_SOURCE_LIST.iteritems():
                if v == 'NONE':
                    num_none += 1
                elif v == 'MOM':
                    num_mom += 1
                elif v == 'DAD':
                    num_dad += 1
                elif v == 'BOTH':
                    num_mom += 1 
                    num_dad += 1
                else: 
                    print "ERROR: cannot understand the source parent = %s" % (v)
                    raise SystemExit 

            if num_none > 0:
                all_from_one_parent = False 
            elif num_mom < tot_num_vars and num_dad < tot_num_vars:
                all_from_one_parent = False  

            # if all variants in the child in this gene are found in single unaffected parent - filter out 
            if all_from_one_parent:
                for kz in GENE_KEY_GT[g]:
                    print "***[Biallelic model]*** Excluded CHILD HET var %s in gene = %s, found in = %s, PARENT_STAT = UNAFFECTED" % (kz,g,VAR_SOURCE_LIST[kz])
                continue

        # end processing all HET variants in the proband - if all from single unaffected parent they have been excluded, message to the log written
        # and gone to evaluating the next biallelic gene in the child

        # if here
        # - either not all CHILD variants in this gene are not HET, or
        # - not all of them can be traced to a single unaffected parent
        # --> add to be processed

        # here we are at gene level, must iterate over all variants in this gene 
        # iterate over variants in this gene
        for kkk in GENE_KEY_GT[g]:                # this the second key: chr:start:end:ref:alt
            
            CHILD_GT = CHILD_DICT['biallelic'][kkk][0]
            CHILD_GENE = CHILD_DICT['biallelic'][kkk][1]
            CHILD_TRANS = CHILD_DICT['biallelic'][kkk][2]

            # if a non-normalized INDEL in child G2P - must adjust (should not happen really, we split, normalized and left-aligned the family VCF before sending it to VEP+G2P)
            chr,start,end,ref,alt = kkk.split(":")
            if len(ref) > 1 and len(alt) > 1:                           # an INDEL - not normalized
                if len(ref) < len(alt):                                 # an INS
                    orig_start = start
                    orig_ref = ref
                    orig_alt = alt
                    start = orig_start
                    ref = '-'
                    alt = orig_alt[len(orig_ref):]
                    print "    WARNING: original INS = %s:%s:%s:%s:%s --> replaced with INS = %s:%s:%s:%s" % (chr,orig_start,end,orig_ref,orig_alt,chr,start,ref,alt)
                else:                                                   # a DEL
                    print "ERROR: At the momemnt, cannot deal with this non-normalized deletion"
                    print line
                    raise SystemExit

            new_key = '%s:%s:%s:%s' % (chr,start,ref,alt)

            # record the data for CHILD G2P variants (for OBS=biallelic)
            if new_key not in G2P_DICT:
                G2P_DICT[new_key] = 0
            else:
                # print "ERROR: duplicate G2P variant new_key = %s" % (new_key)
                # raise SystemExit
                # this will happen if a gene is e.g. hemizygous,x-linked dominant - there will be two separate lines in the output for each req
                pass

            # and record the required data (CHILD_TRANS,CHILD_GENE,CHILD_GT) in G2P_DATA
            if new_key not in G2P_DATA:
                G2P_DATA[new_key] = (CHILD_TRANS,CHILD_GENE,CHILD_GT)
            else:
                # print "ERROR: duplicate G2P variant new_key = %s" % (new_key)
                # raise SystemExit
                # this will happen if a gene is e.g. hemizygous,x-linked dominant - there will be two separate lines in the output for each req
                pass

    NUM_UNIQ_G2P_VARS = len(G2P_DICT)
    print "Found %s unique G2P variants in CHILD (%s) after considering MONOALLELIC, X-LINKED and BIALLELIC genes" % (NUM_UNIQ_G2P_VARS,CHILD_ID)
    sys.stdout.flush()
    print ""
    print ""













def read_ped(in_file):

    global CHILD_ID
    global CHILD_SEX
    global DEC_CHILD_SEX
    global MOM_ID
    global MOM_STAT 
    global DAD_ID
    global DAD_STAT

    CHILD_ID = 0
    CHILD_SEX = 0
    MOM_ID = 0
    MOM_STAT = 0
    DAD_ID = 0
    DAD_STAT = 0

    in_han = open(in_file,'r')
    for line in in_han:
        data = [x.strip() for x in line.strip().split('\t')]
        if data[2] != '0' and data[3] != '0':			# this is the child in the trio
            if CHILD_ID == 0:
                CHILD_ID = data[1]
            else:						# seen another child
                print "ERROR: already have seen a child (possibly a quad) - cannot handle at the moment"  
                raise SystemExit

            if DAD_ID == 0: 
                DAD_ID = data[2] 
            else:
                if data[2] != DAD_ID:
                    print "ERROR: DAD_ID mismatch - from child line dad_id = %s, from dad line dad_id = %s" % (data[2],DAD_ID)
                    raise SystemExit   
            if MOM_ID == 0:
                MOM_ID = data[3]
            else:
                if data[3] != MOM_ID:
                    print "ERROR: MOM_ID mismatch - from child line mom_id = %s, from mom line mom_id = %s" % (data[3],MOM_ID)
                    raise SystemExit

            CHILD_SEX = int(data[4])
            if CHILD_SEX == 1:		# boy
                DEC_CHILD_SEX = '46XY'
            elif CHILD_SEX == 2:	# girl
                DEC_CHILD_SEX = '46XX'
            else:
                print "ERROR: proband sex unknown"
                print line
                raise SystemExit

            if int(data[5]) != 2:
                print "ERROR: child not affected"
                print line
                raise SystemExit  


        elif int(data[2]) == 0 and int(data[3]) == 0:		# this is a parent record
            if int(data[4]) == 1:				# this is the dad
                if int(data[5]) == 1:
                    DAD_STAT = "UNAFFECTED"
                elif int(data[5]) == 2:
                    DAD_STAT = "AFFECTED" 
                else:  
                    print "ERROR: cannot establish the dad's status"
                    print line 
                    raise SystemExit

                if DAD_ID == 0:
                    DAD_ID = data[1]
                else:
                    if data[1] != DAD_ID:
                        print "ERROR: DAD_ID mismatch - from dad line dad_id = %s, from child line dad_id = %s" % (data[1],DAD_ID)
                        raise SystemExit
                
            if int(data[4]) == 2:                               # this is the mom
                if int(data[5]) == 1:
                    MOM_STAT = "UNAFFECTED"
                elif int(data[5]) == 2:
                    MOM_STAT = "AFFECTED"
                else:
                    print "ERROR: cannot establish mom's status"
                    print line
                    raise SystemExit

                if MOM_ID == 0:
                    MOM_ID = data[1]
                else:
                    if data[1] != MOM_ID:
                        print "ERROR: MOM_ID mismatch - from mom line mom_id = %s, from child line mom_id = %s" % (data[1],MOM_ID)
                        raise SystemExit
        else:
            print "ERROR: problematic PED line"
            print line
            raise SystemExit 






def read_map_file(in_file):
    in_han = open(in_file,'r')
    for line in in_han:
        data = [x.strip() for x in line.strip().split('\t')]
        dec_id = data[0]
        int_id = data[1]
        if dec_id not in MAP_DICT:
            MAP_DICT[dec_id] = int_id
        else:
            print "ERROR: duplicate DECIPHER/family ID = %s" % (dec_id)
            raise SystemExit
    in_han.close()  




def read_trans_map(in_file):
    in_han = open(in_file,'r')
    for line in in_han:
        data = [x.strip() for x in line.strip().split('\t')]
        old_trans_id = data[0]
        new_trans_id = data[1]
        if old_trans_id not in TRANS_DICT:
            TRANS_DICT[old_trans_id] = new_trans_id
        else:
            print "ERROR: duplicate old transcript ID = %s" % (old_trans_id)
            raise SystemExit
    in_han.close()






if __name__ == '__main__':
    if len(sys.argv) == 12:
        go(sys.argv[1],sys.argv[2],sys.argv[3],sys.argv[4],sys.argv[5],sys.argv[6],sys.argv[7],sys.argv[8],sys.argv[9],sys.argv[10],sys.argv[11])
    else:
        print "Suggested use: time python /home/u035/u035/shared/scripts/NHS_WES_generate_DEC_IGV.py \
        dec_map_file,trans_map_file,ped_file,in_g2p_file,in_vase_file,fam_igv_dir,vcf_dir,plate_id,fam_id,dec_dir,fam_bam_dir"
        raise SystemExit