dyad_fold
system:sage


{{{id=1|
def hamming(s1,s2):
    counter=0
    l=len(s1)
    for i in range(l):
        if s1[i]!=s2[i]:
            counter+=1
    return counter
///
}}}

{{{id=2|
def hamming_map(s1,s2,dictionary):
    counter=0
    l=len(s1)
    for i in range(l):
        if s1[i]!='N' and s2[i]!='N':
            if s1[i]!=dictionary[s2[i]]:
                counter+=1
    return counter
///
}}}

{{{id=6|
def permute_alph(ab,n):
    l=len(ab)
    if l==1:
        list0=[]
        list0.append(ab[0])
        return list0
        #return dic
    else:
        fact=math.factorial(l-1)
        pos=floor(n/fact)
        #print pos
        rem=n%fact
        #print rem
        letter=ab.pop(pos)       
        temp=list(permute_alph(ab,rem))
        #print temp,'x'
        temp.append(letter)                
        #print temp,'y'
        return temp
    #dic={}
    #for i in range(l):
        #dic[ori[i]]=temp[i]        
    #return dic
///
}}}

{{{id=7|
def permute_dic(ab,n):
    l=len(ab)
    ori=[]
    for obj in ab:
        ori.append(obj)
    temp=permute_alph(ab,n)
    dic={}
    for i in range(l):
        dic[ori[i]]=temp[i]        
    return dic
///
}}}

{{{id=3|
def p_dis(s1,s2):
    dis=len(s1)
    for i in range(120):   
        alp=['A','C','G','T','x']     
        dic=permute_dic(alp,i)
        #print i
        #print dic
        d=hamming_map(s1,s2,dic)
        #print d
        if d<dis:
            dis=d
        #print dis
    return dis
///
}}}

{{{id=11|
def get_column(file):
    lines=file.readlines()
    c_list=[]
    for i in range(99):
        c_list.append('')
    for i in range(99):
        for j in range(100):
            c_list[i]+=lines[j][i]
    file.close()
    return c_list
///
}}}

{{{id=12|
def get_pdm(c_list,outfile):
    for i in range(10):
        outfile.write('c0'+str(i))
        for j in range(99):
            outfile.write(','+str(p_dis(c_list[i],c_list[j])))
        outfile.write('\n')
    for i in range(89):
        outfile.write('c'+str(i+10))
        for j in range(99):
            outfile.write(','+str(p_dis(c_list[i+10],c_list[j])))
        outfile.write('\n')
    outfile.close()
///
}}}

{{{id=4|
file=open('/Users/biocomplexity/Documents/dyadfold/sampled alignments/0.1/1.txt','r')
///
}}}

{{{id=16|
file.close()
///
}}}

{{{id=13|
outfile=open('/Users/biocomplexity/Documents/dyadfold/sampled alignments/0.1/pdm/1.txt','w')
///
}}}

{{{id=15|
for i in range(100):
    file=open('/Users/biocomplexity/Documents/dyadfold/sampled alignments/0.9/'+str(i+1)+'.txt','r')
    outfile=open('/Users/biocomplexity/Documents/dyadfold/sampled alignments/0.9/pdm/'+str(i+1)+'.csv','w')
    c_list=get_column(file)
    get_pdm(c_list,outfile)
///
}}}

{{{id=8|
c_list=get_column(file)
///
}}}

{{{id=14|
get_pdm(c_list,outfile)
///
}}}

{{{id=9|
len(lines[0])
///
100
}}}

{{{id=17|
for j in range(100):
    file=open('/Users/biocomplexity/Documents/dyadfold/sampled alignments/0.9/pdm/'+str(j+1)+'.csv','r')
    outfile=open('/Users/biocomplexity/Documents/dyadfold/tree/0.9/'+str(j+1)+'.txt','w')
    pdm=dendropy.PhylogeneticDistanceMatrix.from_csv(file,is_first_row_column_names=False,is_first_column_row_names=True,is_allow_new_taxa=True,delimiter=",",)
    upgma_tree = pdm.upgma_tree()
    outfile.write(upgma_tree.as_string("newick"))
    outfile.close()
    file.close()
///
}}}

{{{id=18|

///
}}}