comparison assign_clades.py @ 3:bb1cdfafee59 draft default tip

planemo upload for repository https://github.com/Public-Health-Bioinformatics/flu_classification_suite commit 6f09c69c51ec3d6bd7487f55384b97155355c456

2:0d3dad155413

cladeList = [] #empty list to hold clade tuples i.e. ("3C.3a", 1 ,{"3":"I", "9":"V"..})

'''Searches record for required amino acids at defined positions. If found, assigns
clade name to sequence name by appending underscore and clade name to record id.'''
def call_clade(record):
    print "---------------------------------------------------------------------"
    print "Parsing %s for matching flu clade definitions..." % (record.id)
    matchList = [] #empty list to hold clades that match 100%
    #iterate over each tuple in the clade list
    for clade in cladeList:
        cladeName = clade[0] #temp variable for name
        depth = clade[1] #temp variable for depth
        sites = clade[2] #temp variable for aa def dictionary
        shouldFind = len(sites) #number of sites that should match
        found = 0 #a counter to hold matches to antigenic sites
        #iterate over each position in sites dictionary
        for pos, aa in sites.iteritems():
            #translate pos to corresponding index in target sequence
            index = int(pos) - 1
            #if record at index has same amino acid as 'aa', increment 'found'
            if record[index] == aa:
                found += 1

            aa = new_list[i + 1]
            sites[pos] = aa
        #add the clade info as a tuple to the cladeList[]
        oneClade =(name, depth, sites)
        cladeList.append(oneClade)
    print "The List of Clades:"
    for clade in cladeList:
        print "Clade Name: %s Depth: %i Antigenic Sites: %i" % (clade[0], clade[1], len(clade[2]))
        for pos, aa in clade[2].iteritems():
            print "Pos: %s\tAA: %s" % (pos,aa)

'''opens readable input file of sequences to parse using filename from cmd line,
    instantiates as AA Sequence objects, with ppercase sequences'''
with open(inFileHandle1,'r') as inFile:
    #read in Sequences from fasta file, uppercase and add to seqList
    for record in SeqIO.parse(inFile, "fasta", alphabet=IUPAC.protein):
        record = record.upper()
        seqList.append(record) #add Seq to list of Sequences
    print "\n%i flu HA sequences will be compared to current clade definitions..." % len(seqList)
    #parse each target sequence object
    for record in seqList:
        clade_call = '' #empty variale for final clade call on sequence
        matchingCladeList = call_clade(record) #holds matching clade tuples
        #if there is more than one clade match

            clade = matchingCladeList[0] #take the first tuple in the list
            clade_call = clade[0] #clade name is the first item in the tuple
        #empty list return, no matches
        else:
            clade_call = "No_Match"
        print clade_call
        seq_name = record.id
        mod_name = seq_name + "_" + clade_call
        print "New Sequence Name: " + mod_name
        record.id = mod_name


#output fasta file with clade calls appended to sequence names
SeqIO.write(seqList,outFile,"fasta")

3:bb1cdfafee59

cladeList = [] #empty list to hold clade tuples i.e. ("3C.3a", 1 ,{"3":"I", "9":"V"..})

'''Searches record for required amino acids at defined positions. If found, assigns
clade name to sequence name by appending underscore and clade name to record id.'''
def call_clade(record):
    print("---------------------------------------------------------------------")
    print("Parsing %s for matching flu clade definitions..." % (record.id))
    matchList = [] #empty list to hold clades that match 100%
    #iterate over each tuple in the clade list
    for clade in cladeList:
        cladeName = clade[0] #temp variable for name
        depth = clade[1] #temp variable for depth
        sites = clade[2] #temp variable for aa def dictionary
        shouldFind = len(sites) #number of sites that should match
        found = 0 #a counter to hold matches to antigenic sites
        #iterate over each position in sites dictionary
        for pos, aa in sites.items():
            #translate pos to corresponding index in target sequence
            index = int(pos) - 1
            #if record at index has same amino acid as 'aa', increment 'found'
            if record[index] == aa:
                found += 1

            aa = new_list[i + 1]
            sites[pos] = aa
        #add the clade info as a tuple to the cladeList[]
        oneClade =(name, depth, sites)
        cladeList.append(oneClade)
    print("The List of Clades:")
    for clade in cladeList:
        print("Clade Name: %s Depth: %i Antigenic Sites: %i" % (clade[0], clade[1], len(clade[2])))
        for pos, aa in clade[2].items():
            print("Pos: %s\tAA: %s" % (pos,aa))

'''opens readable input file of sequences to parse using filename from cmd line,
    instantiates as AA Sequence objects, with ppercase sequences'''
with open(inFileHandle1,'r') as inFile:
    #read in Sequences from fasta file, uppercase and add to seqList
    for record in SeqIO.parse(inFile, "fasta", alphabet=IUPAC.protein):
        record = record.upper()
        seqList.append(record) #add Seq to list of Sequences
    print("\n%i flu HA sequences will be compared to current clade definitions..." % len(seqList))
    #parse each target sequence object
    for record in seqList:
        clade_call = '' #empty variale for final clade call on sequence
        matchingCladeList = call_clade(record) #holds matching clade tuples
        #if there is more than one clade match

            clade = matchingCladeList[0] #take the first tuple in the list
            clade_call = clade[0] #clade name is the first item in the tuple
        #empty list return, no matches
        else:
            clade_call = "No_Match"
        print(clade_call)
        seq_name = record.id
        mod_name = seq_name + "_" + clade_call
        print("New Sequence Name: " + mod_name)
        record.id = mod_name


#output fasta file with clade calls appended to sequence names
SeqIO.write(seqList,outFile,"fasta")