<filename>exomeCRG/SNP_applyFamilyModel_2.py
#!/usr/bin/python

import argparse
import csv
import sys
import tools

def main(args):
   
    #### argument reading
    try:
        snpFH = open(sys.argv[1])
    except IOError:
        sys.exit("ERROR: could not open SNP file")
    except IndexError:
        usage()
        sys.exit("ERROR: not enough arguments given")
    
    try:
        famFH = open(sys.argv[2])
    except IOError:
        sys.exit("ERROR: could not open Fam file")
    except IndexError:
        usage()
        sys.exit("ERROR: not enough arguments given")
    
    try:
        job = sys.argv[3]
    except IOError:
        usage()
        sys.exit("ERROR: not enough arguments given")
    
    if job == 'compound':
        try:
            cutoff = float(sys.argv[4])
        except IndexError:
            print 'setting MAF cutoff to default = 1'
            cutoff = float(1)
            
    
    try:
        outname = '%s.familyfiltered' % sys.argv[1]
        outFH = open(outname, 'w')
    except IOError:
        sys.exit("ERROR: could not open output file")
    
    #### read family file
    family = dict()
    affection = dict()
    
    for line in famFH:
        if line.startswith('family'):
            continue
        
        splitline = line.split()
        
        family[splitline[1]] = splitline[0]
        affection[splitline[1]] = splitline[2]
    
    #### read snp file and judge family info and SNP distribution
    csvfile = csv.reader(snpFH)
    
    index_SNPhighqual = 0
    index_SNPlowqual = 0
    index_Gene = 0
    index_HomHet = 0
    index_Refhom = 0
    index_Refhet = 0
    index_EVSMAF = 0
    index_1000GMAF = 0
    
    genesHit = dict()
    
    if job == 'denovo':
            print 'highlighting de novo variants'
    elif job == 'recfamily':
        print 'highlighting variants fitting to a recessive disease model'
    elif job == 'domfamily':
        print 'highlighting variants fitting to a dominant disease model'
    elif job == 'homozygous':
        print 'looking for variants that are het in parents and hom in offspring'
    elif job == 'compound':
        print 'looking for variants affecting the same gene, applying a cutoff of %f' % cutoff
    else:
        sys.exit('job id not understood, sorry. Select one of: denovo, domfamily, recfamily, homozygous, compound')

    
    for row in csvfile:
        # find positions of some headers
        if row[0] == 'Chr':
            header = row[0:9]
            header.append(job)
            #header.append('unfit')
            header.extend(row[9:])
            header.append('\n')
            outFH.write(','.join(header))
            
            index_SNPhighqual = row.index('SNPcall:highQual')
            index_SNPlowqual = row.index('SNPcall:lowQual')
            index_Gene = row.index('Gene')
            index_HomHet = row.index('SampleDetails(SNP[0|1|2]Ref[0|1|2])')
            index_Refhom = row.index('Refcall:hom')
            index_Refhet = row.index('Refcall:het')
            try:
                index_EVSMAF = row.index('EVS_MAF_TAC')# indels do not harbor Allele Frequencies
                index_1000GMAF = row.index('MAF_1000G')
            except:
                index_EVSMAF = 'NA'
                index_1000GMAF = 'NA'
            
            continue
        
        SNPhighqual = row[index_SNPhighqual]
        SNPlowqual = row[index_SNPlowqual]
        Refhom = row[index_Refhom]
        Refhet = row[index_Refhet]
        
        HomHet = row[index_HomHet]
        Gene = row[index_Gene]
        
        if index_EVSMAF != 'NA' and index_1000GMAF != 'NA':
            EVSMAF = float(row[index_EVSMAF])
            MAF1000G = float(row[index_1000GMAF])
            
            MAF = max(EVSMAF, MAF1000G) # use the maximum MAF of EVS or 1000Genomes
            
        else:
            MAF = 0 # indels do not harbor Allele Frequencies that's why index_MAF is 'NA'
        
        
        
        
        detailed_family = dict()
        
        if job == 'denovo':
            #print 'highlighting de novo variants'
            detailed_family = tools.denovo(SNPhighqual, SNPlowqual, Refhom, Refhet, affection, family)
        elif job == 'recfamily':
            #print 'highlighting variants fitting to a recessive disease model'
            #sys.exit('currently not implemented')
            detailed_family = tools.recfamily(SNPhighqual, SNPlowqual, Refhom, Refhet, HomHet, affection, family)
        elif job == 'domfamily':
            #print 'highlighting variants fitting to a dominant disease model'
            detailed_family = tools.domfamily(SNPhighqual, SNPlowqual, Refhom, Refhet, HomHet, affection, family)
        elif job == 'homozygous':
            #print 'looking for variants that are het in parents and hom in offspring'
            detailed_family = tools.homozygous(HomHet, affection, family)
        elif job == 'compound':
            #print 'looking for variants affecting the same gene'
            # this has to be done in another part of the script...
            genesHit = tools.compoundGenelist(HomHet, Gene, affection, family, genesHit, MAF, cutoff)
        else:
            sys.exit('job id not understood, sorry. Select one of: denovo, domfamily, recfamily, homozygous, compound')
        
        if job != 'compound':
            #write output to file
            fit = list()
            unfit = list()
            
            for element in detailed_family: # element is the family ID given in pedigree information
                if detailed_family[element] == 1:
                    fit.append(element)
                elif detailed_family[element] == 0:
                    unfit.append(element)
                else:
                    sys.stderr.write('error in assigning pattern status for %s %s \n' % (row[0], row[1]))
            
            fit = tools.fillempty(fit)
            unfit = tools.fillempty(unfit)
            
            printrow = row[0:9]
            printrow.append(' '.join(map(str, fit)))
            printrow.extend(row[9:])
            printrow.append('\n')
            outFH.write(','.join(map(str, printrow)))
    
    
    if job == 'compound':
        
        #identify potential compound hets
        potentialCompounds = tools.compound(affection, family, genesHit)
        
        ### compound het:
        # compound het is a gene that is hit at least twice in het state in the affected
        # and is het in one of the unaffected parents at these positions
        # rewind file handle
        snpFH.seek(0)
        
        all_potentials = list()
        
        for fam in potentialCompounds:
            all_potentials.extend(potentialCompounds[fam])
    
        for row in csvfile:
            # find positions of some headers
            if row[0] == 'Chr':
                
                index_SNPhighqual = row.index('SNPcall:highQual')
                index_SNPlowqual = row.index('SNPcall:lowQual')
                index_Gene = row.index('Gene')
                index_HomHet = row.index('SampleDetails(SNP[0|1|2]Ref[0|1|2])')
                index_Refhom = row.index('Refcall:hom')
                index_Refhet = row.index('Refcall:het')
                
                try:
                    index_EVSMAF = row.index('EVS_MAF_TAC')# indels do not harbor Allele Frequencies
                    index_1000GMAF = row.index('MAF_1000G')
                except:
                    index_EVSMAF = 'NA'
                    index_1000GMAF = 'NA'
                
                continue
            
            # read values
            HomHet = row[index_HomHet]
            Gene = row[index_Gene]
            
            if index_EVSMAF != 'NA' and index_1000GMAF != 'NA':
                EVSMAF = float(row[index_EVSMAF])
                MAF1000G = float(row[index_1000GMAF])
                
                MAF = max(EVSMAF, MAF1000G) # use the maximum MAF of EVS or 1000Genomes
                
            else:
                MAF = 0 # indels do not harbor Allele Frequencies that's why index_MAF is 'NA'
            
            
            if Gene in all_potentials:
                printrow = row[0:9]
                if (MAF > cutoff):
                    printrow.append('0')
                else:
                    printrow.append('1')
                printrow.extend(row[9:])
                printrow.append('\n')
                outFH.write(','.join(map(str, printrow)))
            else:
                printrow = row[0:9]
                printrow.append('0')
                printrow.extend(row[9:])
                printrow.append('\n')
                outFH.write(','.join(map(str, printrow)))
        
    snpFH.close()
    famFH.close()
    outFH.close()
    
def usage():
    print "usage: \n python SNP_applyFamilyModel_2.py file pedigree model \n\n"
    print "# file: output of SNP_comparison_1.6.pl \n"
    print "# pedigree: file containing all family members following the pattern"
    print "@ family\tsample\taffected\n@ 1\tVH001\t1\n@ 1\tVH002\t1\n"
    print "# model: denovo|recfamily|domfamily|compound"
    print "### denovo will add a column containing 1 where the script thinks it identified a denovo variant - works in trios only"
    print "### homozygous will add a column containing 1 where the script thinks it identified a variant that is heterozygous in parents, but homozygous in offspring"
    print "### homozygous - accepts a MAF cutoff to exclude SNPs in compound finding with a MAF > cutoff."
    print "### recfamily will add a column containing 1 where the script thinks it identified a homozygous variant, which is not in unaffected individuals"
    print "### domfamily will add a column highlighting variants that are present in affected family members"
    print "### compound will add a column containing 1 where both parents supply a variant in the same gene, which is heterozygous in all investigated individuals - works in trios only\n\n"

main(sys.argv)
<filename>exomeCRG/SNP_applyFamilyModel_2.py
#!/usr/bin/python

import argparse
import csv
import sys
import tools

def main(args):
   
    #### argument reading
    try:
        snpFH = open(sys.argv[1])
    except IOError:
        sys.exit("ERROR: could not open SNP file")
    except IndexError:
        usage()
        sys.exit("ERROR: not enough arguments given")
    
    try:
        famFH = open(sys.argv[2])
    except IOError:
        sys.exit("ERROR: could not open Fam file")
    except IndexError:
        usage()
        sys.exit("ERROR: not enough arguments given")
    
    try:
        job = sys.argv[3]
    except IOError:
        usage()
        sys.exit("ERROR: not enough arguments given")
    
    if job == 'compound':
        try:
            cutoff = float(sys.argv[4])
        except IndexError:
            print 'setting MAF cutoff to default = 1'
            cutoff = float(1)
            
    
    try:
        outname = '%s.familyfiltered' % sys.argv[1]
        outFH = open(outname, 'w')
    except IOError:
        sys.exit("ERROR: could not open output file")
    
    #### read family file
    family = dict()
    affection = dict()
    
    for line in famFH:
        if line.startswith('family'):
            continue
        
        splitline = line.split()
        
        family[splitline[1]] = splitline[0]
        affection[splitline[1]] = splitline[2]
    
    #### read snp file and judge family info and SNP distribution
    csvfile = csv.reader(snpFH)
    
    index_SNPhighqual = 0
    index_SNPlowqual = 0
    index_Gene = 0
    index_HomHet = 0
    index_Refhom = 0
    index_Refhet = 0
    index_EVSMAF = 0
    index_1000GMAF = 0
    
    genesHit = dict()
    
    if job == 'denovo':
            print 'highlighting de novo variants'
    elif job == 'recfamily':
        print 'highlighting variants fitting to a recessive disease model'
    elif job == 'domfamily':
        print 'highlighting variants fitting to a dominant disease model'
    elif job == 'homozygous':
        print 'looking for variants that are het in parents and hom in offspring'
    elif job == 'compound':
        print 'looking for variants affecting the same gene, applying a cutoff of %f' % cutoff
    else:
        sys.exit('job id not understood, sorry. Select one of: denovo, domfamily, recfamily, homozygous, compound')

    
    for row in csvfile:
        # find positions of some headers
        if row[0] == 'Chr':
            header = row[0:9]
            header.append(job)
            #header.append('unfit')
            header.extend(row[9:])
            header.append('\n')
            outFH.write(','.join(header))
            
            index_SNPhighqual = row.index('SNPcall:highQual')
            index_SNPlowqual = row.index('SNPcall:lowQual')
            index_Gene = row.index('Gene')
            index_HomHet = row.index('SampleDetails(SNP[0|1|2]Ref[0|1|2])')
            index_Refhom = row.index('Refcall:hom')
            index_Refhet = row.index('Refcall:het')
            try:
                index_EVSMAF = row.index('EVS_MAF_TAC')# indels do not harbor Allele Frequencies
                index_1000GMAF = row.index('MAF_1000G')
            except:
                index_EVSMAF = 'NA'
                index_1000GMAF = 'NA'
            
            continue
        
        SNPhighqual = row[index_SNPhighqual]
        SNPlowqual = row[index_SNPlowqual]
        Refhom = row[index_Refhom]
        Refhet = row[index_Refhet]
        
        HomHet = row[index_HomHet]
        Gene = row[index_Gene]
        
        if index_EVSMAF != 'NA' and index_1000GMAF != 'NA':
            EVSMAF = float(row[index_EVSMAF])
            MAF1000G = float(row[index_1000GMAF])
            
            MAF = max(EVSMAF, MAF1000G) # use the maximum MAF of EVS or 1000Genomes
            
        else:
            MAF = 0 # indels do not harbor Allele Frequencies that's why index_MAF is 'NA'
        
        
        
        
        detailed_family = dict()
        
        if job == 'denovo':
            #print 'highlighting de novo variants'
            detailed_family = tools.denovo(SNPhighqual, SNPlowqual, Refhom, Refhet, affection, family)
        elif job == 'recfamily':
            #print 'highlighting variants fitting to a recessive disease model'
            #sys.exit('currently not implemented')
            detailed_family = tools.recfamily(SNPhighqual, SNPlowqual, Refhom, Refhet, HomHet, affection, family)
        elif job == 'domfamily':
            #print 'highlighting variants fitting to a dominant disease model'
            detailed_family = tools.domfamily(SNPhighqual, SNPlowqual, Refhom, Refhet, HomHet, affection, family)
        elif job == 'homozygous':
            #print 'looking for variants that are het in parents and hom in offspring'
            detailed_family = tools.homozygous(HomHet, affection, family)
        elif job == 'compound':
            #print 'looking for variants affecting the same gene'
            # this has to be done in another part of the script...
            genesHit = tools.compoundGenelist(HomHet, Gene, affection, family, genesHit, MAF, cutoff)
        else:
            sys.exit('job id not understood, sorry. Select one of: denovo, domfamily, recfamily, homozygous, compound')
        
        if job != 'compound':
            #write output to file
            fit = list()
            unfit = list()
            
            for element in detailed_family: # element is the family ID given in pedigree information
                if detailed_family[element] == 1:
                    fit.append(element)
                elif detailed_family[element] == 0:
                    unfit.append(element)
                else:
                    sys.stderr.write('error in assigning pattern status for %s %s \n' % (row[0], row[1]))
            
            fit = tools.fillempty(fit)
            unfit = tools.fillempty(unfit)
            
            printrow = row[0:9]
            printrow.append(' '.join(map(str, fit)))
            printrow.extend(row[9:])
            printrow.append('\n')
            outFH.write(','.join(map(str, printrow)))
    
    
    if job == 'compound':
        
        #identify potential compound hets
        potentialCompounds = tools.compound(affection, family, genesHit)
        
        ### compound het:
        # compound het is a gene that is hit at least twice in het state in the affected
        # and is het in one of the unaffected parents at these positions
        # rewind file handle
        snpFH.seek(0)
        
        all_potentials = list()
        
        for fam in potentialCompounds:
            all_potentials.extend(potentialCompounds[fam])
    
        for row in csvfile:
            # find positions of some headers
            if row[0] == 'Chr':
                
                index_SNPhighqual = row.index('SNPcall:highQual')
                index_SNPlowqual = row.index('SNPcall:lowQual')
                index_Gene = row.index('Gene')
                index_HomHet = row.index('SampleDetails(SNP[0|1|2]Ref[0|1|2])')
                index_Refhom = row.index('Refcall:hom')
                index_Refhet = row.index('Refcall:het')
                
                try:
                    index_EVSMAF = row.index('EVS_MAF_TAC')# indels do not harbor Allele Frequencies
                    index_1000GMAF = row.index('MAF_1000G')
                except:
                    index_EVSMAF = 'NA'
                    index_1000GMAF = 'NA'
                
                continue
            
            # read values
            HomHet = row[index_HomHet]
            Gene = row[index_Gene]
            
            if index_EVSMAF != 'NA' and index_1000GMAF != 'NA':
                EVSMAF = float(row[index_EVSMAF])
                MAF1000G = float(row[index_1000GMAF])
                
                MAF = max(EVSMAF, MAF1000G) # use the maximum MAF of EVS or 1000Genomes
                
            else:
                MAF = 0 # indels do not harbor Allele Frequencies that's why index_MAF is 'NA'
            
            
            if Gene in all_potentials:
                printrow = row[0:9]
                if (MAF > cutoff):
                    printrow.append('0')
                else:
                    printrow.append('1')
                printrow.extend(row[9:])
                printrow.append('\n')
                outFH.write(','.join(map(str, printrow)))
            else:
                printrow = row[0:9]
                printrow.append('0')
                printrow.extend(row[9:])
                printrow.append('\n')
                outFH.write(','.join(map(str, printrow)))
        
    snpFH.close()
    famFH.close()
    outFH.close()
    
def usage():
    print "usage: \n python SNP_applyFamilyModel_2.py file pedigree model \n\n"
    print "# file: output of SNP_comparison_1.6.pl \n"
    print "# pedigree: file containing all family members following the pattern"
    print "@ family\tsample\taffected\n@ 1\tVH001\t1\n@ 1\tVH002\t1\n"
    print "# model: denovo|recfamily|domfamily|compound"
    print "### denovo will add a column containing 1 where the script thinks it identified a denovo variant - works in trios only"
    print "### homozygous will add a column containing 1 where the script thinks it identified a variant that is heterozygous in parents, but homozygous in offspring"
    print "### homozygous - accepts a MAF cutoff to exclude SNPs in compound finding with a MAF > cutoff."
    print "### recfamily will add a column containing 1 where the script thinks it identified a homozygous variant, which is not in unaffected individuals"
    print "### domfamily will add a column highlighting variants that are present in affected family members"
    print "### compound will add a column containing 1 where both parents supply a variant in the same gene, which is heterozygous in all investigated individuals - works in trios only\n\n"

main(sys.argv)
#!/usr/bin/python #### argument reading #### read family file #### read snp file and judge family info and SNP distribution # find positions of some headers #header.append('unfit') # indels do not harbor Allele Frequencies # use the maximum MAF of EVS or 1000Genomes # indels do not harbor Allele Frequencies that's why index_MAF is 'NA' #print 'highlighting de novo variants' #print 'highlighting variants fitting to a recessive disease model' #sys.exit('currently not implemented') #print 'highlighting variants fitting to a dominant disease model' #print 'looking for variants that are het in parents and hom in offspring' #print 'looking for variants affecting the same gene' # this has to be done in another part of the script... #write output to file # element is the family ID given in pedigree information #identify potential compound hets ### compound het: # compound het is a gene that is hit at least twice in het state in the affected # and is het in one of the unaffected parents at these positions # rewind file handle # find positions of some headers # indels do not harbor Allele Frequencies # read values # use the maximum MAF of EVS or 1000Genomes # indels do not harbor Allele Frequencies that's why index_MAF is 'NA' ## denovo will add a column containing 1 where the script thinks it identified a denovo variant - works in trios only" ## homozygous will add a column containing 1 where the script thinks it identified a variant that is heterozygous in parents, but homozygous in offspring" ## homozygous - accepts a MAF cutoff to exclude SNPs in compound finding with a MAF > cutoff." ## recfamily will add a column containing 1 where the script thinks it identified a homozygous variant, which is not in unaffected individuals" ## domfamily will add a column highlighting variants that are present in affected family members" ## compound will add a column containing 1 where both parents supply a variant in the same gene, which is heterozygous in all investigated individuals - works in trios only\n\n"