pycivo.py

from __future__ import absolute_import
import numpy as np
import networkx as nx
import math, sys, os, time, subprocess, itertools, fastcluster, random
from scipy.sparse.csgraph import connected_components,csgraph_from_dense,dijkstra 
from scipy.stats import rankdata, hmean, zscore
from scipy.spatial.distance import pdist

from tempfile import NamedTemporaryFile
from collections import defaultdict
import lsqlin
from cvxopt.solvers import options

import warnings
from Bio import BiopythonWarning
warnings.simplefilter('ignore', BiopythonWarning)
from Bio import SeqIO

__author__ = "J Walker Gussler"

"""PYCIVO (INference of Super-spreaders through Computational Analysis of quasisPEcies), by J Walker Gussler
Given a number of .fas clinical HCV samples of HVR1 sequenced by the GHOST protocol, this program will determine the presence of any super-spreaders, if applicable
 
1. It requires as input the path to the aligned fasta files. All sequences are assumed to be unique.
2. From these sequences, 2 indices are computed for each clinical sample.
    a) A value based off evolution simulations which is meant to represent the likelihood of a person's HCV set evolving into the others listed
    b) A value of estimated chronicity based off several metrics (shannon entropy, nucleotide diversity).
3. If these two methods agree and the predicted source is different enough from the rest of the cases, an output is generated.
4. Output is generated by software similar to QUENTIN by Pavel Skums, implemented in Matlab
"""

def pre_align(f1,f2): #align with mafft and tempfiles 
    count=0
    with open(f1) as f:
        for record in SeqIO.parse(f,'fasta'):
            count+=1
    with NamedTemporaryFile(delete=False) as tmpcat:
        catname=tmpcat.name
        try:
            subprocess.check_call(['cat',f1,f2],stdout=tmpcat)
        except:
            print('Invalid file. Check -s argument if in csv mode')
            return 0
    
    with NamedTemporaryFile(delete=False) as aligned:
        alignname=aligned.name
        subprocess.check_call(['mafft', '--quiet', '--auto', '--thread', '20', '--preservecase', catname], stdout=aligned) 
    os.unlink(catname)
    seqs=SeqIO.parse(alignname,'fasta')
    seenSecond=False
    seqs1=[]
    c=0
    while not seenSecond:
        record=next(seqs)
        c+=1
        if c==count:
            seenSecond=True
        else:
            seqs1.append(record)
    # import ipdb
    # ipdb.set_trace()
    with NamedTemporaryFile(mode='w',delete=False) as align1:
        SeqIO.write(seqs1,align1,'fasta')
        f1name=align1.name
    seqs2=[]
    seqs2.append(record)#write record already accessed and ignored because it doesn't belong in first file
    
    done=False
    while not done:     
        try:
            record=next(seqs)
            seqs2.append(record)
        except StopIteration:
            done=True
    with NamedTemporaryFile(mode='w',delete=False) as align2:
        SeqIO.write(seqs2,align2,'fasta')
        f2name=align2.name
    os.unlink(alignname)    
    return(f1name,f2name)

def parse_input(input): #get sequences from a file
    seqs={}
    with open(input,'r') as f:
        for record in SeqIO.parse(f,'fasta'):
            freq = int(record.id.split('_')[-1])
            if record.seq not in seqs:
                seqs[record.seq]=freq
            else:
                seqs[record.seq]+=freq
    haploSize=len(list(seqs.keys())[0])
    for seq in seqs:
        if len(seq)!=haploSize:
            return False #all seqs must be same len. required by calc_map_val and calc_centroids
    return seqs

def calc_map_val(seqs): #calculates amount of heterogeneity in sample, used in evolution simulation
    numSeq=len(seqs)
    numPos=len(seqs[0])
    mapval=0;a=0;c=0;g=0;t=0;blank=0
    
    for pos in range(numPos,0,-1):
        for seq in seqs:
            try:
                if seq[pos-1].lower()=="a":
                    a+=1
                elif seq[pos-1].lower()=="c":
                    c+=1
                elif seq[pos-1].lower()=="t":
                    t+=1
                elif seq[pos-1].lower()=="g":
                    g+=1
                else:
                    blank+=1
            except IndexError:
                sys.exit("It looks like your data are not aligned - try again with -a!")
        consensus=max(a,c,t,g,blank)
        if a+c+t+g+blank<numSeq or consensus>=numSeq or blank>0:
            a=0;c=0;g=0;t=0;blank=0;
        else:
            mapval+=1;a=0;c=0;g=0;t=0;blank=0;
    return mapval
    
def calc_ordered_frequencies(haploNum,haploSize,seqs,byFreq): 
    freqCount = np.zeros((haploSize, 5))
    productVector = np.zeros((haploSize, 5))
    order={'A':0,'C':1,'G':2,'T':3,'-':4}
    try:
        total_reads=0
        for read in seqs:
            if byFreq:
                freq=seqs[read]
            else:
                freq=1
            total_reads+=freq
            for pos in range(haploSize):
                num=order[read[pos]]
                freqCount[pos, num] = freqCount[pos, num] + freq
        freqRel = np.divide(freqCount, float(total_reads), dtype = float)
    
    except IndexError:
        print("Your files are not aligned and it caused an error! Try again with -a")
    
    for pos in range(haploSize):
        for i in range(5):
            freqPos = freqRel[pos, i]
            if freqPos > 0:
                logFreqRel = math.log(freqPos, 2)
                productVector[pos, i] = -1*(np.multiply(freqPos, logFreqRel, dtype = float))                
    return np.sum(productVector, axis = 1)

def order_positions(hVector,seqs,haploSize): #order positions for faster building of k-step network
    invH =  np.multiply(-1, hVector, dtype = float)
    ordH = np.argsort(invH)
    # reorder the sequences by entropy
    ordSeqs = []

    for a in seqs:
        i=str(a)
        newOne = ''
        for p in range(haploSize):
            newOne = ''.join([newOne, i[ordH[p]]])
        ordSeqs.append(newOne)
    return ordSeqs

def calc_kstep(haploNum,haploSize,ordSeqs): #Calculate hamming distances between sequences
    compNum = haploSize
    compList = range(haploNum)
    t = 0
    adjMatrix = np.zeros((haploNum, haploNum))
    kStepList = []

    while compNum > 1:
        t = t + 1
        for r1 in range(haploNum-1):
            haplotype1 = ordSeqs[r1]
            for r2 in range(r1+1, haploNum):
                if compList[r1] != compList[r2]: 
                    haplotype2 = ordSeqs[r2]
                    tempDist = 0
                    for a, b in zip(haplotype1, haplotype2):
                        if a != b:
                            tempDist = tempDist + 1
                            if tempDist > t:
                                break
                                
                    if tempDist == t: 
                        adjMatrix[r1][r2] = 1
                        kStepList.append([r1, r2, t])
        sparseMatrix = csgraph_from_dense(adjMatrix)
        connected = connected_components(sparseMatrix, directed=False, connection='weak', return_labels=True)
        compNum = connected[0]
        compList = connected[1]
    return kStepList

def get_intermediate_sequences(kStepList,seqs,haploSize): #make 1-step list of sequences based off k-step list, return list of sequences
    fakeSeqs={}
    for element in kStepList:
        dist=element[2]
        if dist!=1:
            #element[1] -> id for seq1
            #seqs[element[1]] -> sequence for seq1
            s1=seqs[element[0]]
            s2=seqs[element[1]]

            id=0
            while dist>1:
                id+=1
                for nucl in range(haploSize):
                    if s1[nucl]!=s2[nucl]:
                        takeseq1=s1[:nucl+1]+s2[nucl+1:]
                        takeseq2=s1[:nucl]+s2[nucl:]
                        if takeseq1!=s1:
                            newseq=takeseq1
                        else:
                            newseq=takeseq2
                dist-=1
                s1=newseq
                
                if newseq not in seqs and newseq not in fakeSeqs.values():
                    newseqid=str(element[0])+"->"+str(element[1])+"_"+str(id)
                    fakeSeqs[newseqid]=newseq                

    #make sequence dictionary
    finalSeqs=[]
    for item in range(len(seqs)):
        finalSeqs.append(seqs[item])
    finalSeqs.extend(fakeSeqs.values())
    return finalSeqs

def simul_evol(D_all,nseq1,nseq2,len_eff,PERCENT,FAILTIME): #viral evolution simulator, works, probably could be changed
    #possible extension: take in starting frequencies instead of defaulting to 1
    ev=nseq2*float(PERCENT)/100.0
    maxPopl=10**12
    timeInter=FAILTIME-10 # -10 is somewhat arbitrary
    mutprob=.01
    evolved=False
    sz=len(D_all);
    Q=np.zeros([sz,sz],float)
    for u in range(sz):
        Q[u,u]=(mutprob/3)**D_all[u,u]*(1-mutprob)**(len_eff-D_all[u,u])
        for v in range(u+1,sz):
            Q[u,v]=D_all[u,v]*(mutprob/3)**(D_all[u,v])*(1-mutprob)**(len_eff-D_all[u,v])
            Q[v,u]=Q[u,v]
    x=np.zeros([sz,timeInter],float)
    x[0:nseq1,0]=1
    E=np.eye(sz)
    for t in range(1,timeInter):
        x[:,t]=(1-sum(x[:,t-1])/maxPopl) * np.dot((E+Q),x[:,t-1])
        where=np.where(x[nseq1:nseq1+nseq2,t]>=1)
        if len(where[0])>=ev:
            evolved=True
            break
    if evolved==True:
        time=t+1
    if evolved==False:
        time=FAILTIME
    return time

def get_evol_time(haploNum,haploSize,seqs,f1HaploNum,PERCENT,FAILTIME,NOGHOST): #get evolution time one way from already parsed sequences
    f2HaploNum=haploNum-f1HaploNum
    len_eff=calc_map_val(seqs) 
    hVector=calc_ordered_frequencies(haploNum,haploSize,seqs,False)
    ordSeqs=order_positions(hVector,seqs,haploSize)
    
    kStepList=calc_kstep(haploNum,haploSize,ordSeqs)
    
    finalSeqs=get_intermediate_sequences(kStepList,seqs,haploSize)
    if NOGHOST: # this part is commented out because it breaks my test for some reason
        D_all=calc_distance_matrix_slow(finalSeqs)
    else:
        D_all=calc_distance_matrix(finalSeqs)
    evolTime=simul_evol(D_all,f1HaploNum,f2HaploNum,len_eff,PERCENT,FAILTIME)
    return evolTime

def trimfh(file):
    return os.path.splitext(os.path.basename(file))[0]

def combine_seqs(seqs1,seqs2):
    outseqs2=[]
    for item in seqs2:
        if item not in seqs1:
            outseqs2.append(item)
    if len(seqs1+outseqs2)!=len(set(seqs1+outseqs2)):
        a=len(seqs1+outseqs2)-len(set(seqs1+outseqs2))
        sys.exit('error in combineseqs! '+str(a))
    f1HaploNum=len(seqs1)
    return (seqs1+outseqs2,f1HaploNum)
    
class get_evol_times():
    def __init__(self):
        self.ALI=ALI
        self.FAILTIME=FAILTIME
        self.PERCENT=PERCENT
        self.NOGHOST=NOGHOST     
        self.INPUTS=INPUTS
        self.CLUSTERS=CLUSTERS
    
    def __call__(self,i): #get evolution times both ways from fasta files (wrapper for 1-way evolution time func)
        f1=self.INPUTS[i[0]]
        f2=self.INPUTS[i[1]]
        if ALI:
            input1,input2=pre_align(f1,f2)
        else:
            input1=f1
            input2=f2
        seqs1=calc_centroids(self.CLUSTERS,input1,NOGHOST)
        if seqs1==False:
            return False
        seqs2=calc_centroids(self.CLUSTERS,input2,NOGHOST)
        
        seqs,f1HaploNum=combine_seqs(seqs1,seqs2) 
        haploSize=len(seqs[0])
        haploNum=len(seqs)
        forwards=get_evol_time(haploNum,haploSize,seqs,f1HaploNum,self.PERCENT,self.FAILTIME,self.NOGHOST)
        seqs,f2HaploNum=combine_seqs(seqs2,seqs1)
        backwards=get_evol_time(haploNum,haploSize,seqs,f2HaploNum,self.PERCENT,self.FAILTIME,self.NOGHOST)
        return forwards,backwards,i[0],i[1]

def get_evol_times_serial(i1,i2,INPUTS,CLUSTERS,PERCENT,FAILTIME,NOGHOST,ALI):
    f1=INPUTS[i1]
    f2=INPUTS[i2]
    if ALI:
        input1,input2=pre_align(f1,f2)
    else:
        input1=f1
        input2=f2
    seqs1=calc_centroids(CLUSTERS,input1,NOGHOST)
    if seqs1==False:
        return False
    seqs2=calc_centroids(CLUSTERS,input2,NOGHOST)
    seqs,f1HaploNum=combine_seqs(seqs1,seqs2)
    haploSize=len(seqs[0])
    haploNum=len(seqs)
    forwards=get_evol_time(haploNum,haploSize,seqs,f1HaploNum,PERCENT,FAILTIME,NOGHOST)
    seqs,f2HaploNum=combine_seqs(seqs2,seqs1)
    backwards=get_evol_time(haploNum,haploSize,seqs,f2HaploNum,PERCENT,FAILTIME,NOGHOST)
    return forwards,backwards

def calc_distance_matrix(finalSeqs): #calculate distance matrix from the 1-step list
    from pyseqdist import hamming
    l=len(finalSeqs)
    arr=np.zeros([l,l])
    hdist=hamming(finalSeqs,finalSeqs,ignore_gaps=False)
    for id in range(len(hdist)):
        item=hdist[id]
        arr[:,id]=item[:,0]
    return arr

def calc_distance_matrix_slow(finalSeqs): #calculate distance matrix from the 1-step list manually
    l=len(finalSeqs)
    arr=np.zeros([l,l])
    for id1,id2 in itertools.combinations(range(len(finalSeqs)),2):
        seq1=finalSeqs[id1]
        seq2=finalSeqs[id2]
        dist=sum(0 if a==b else 1 for a,b in zip(seq1,seq2))
        arr[id1,id2]=dist
        arr[id2,id1]=dist
    return arr

def outdegree_source(DSamp,files,QUIET): #reutrns dictionary with each file as key and evolution factor as value
    times={}
    for id in range(len(files)):
        evolutionFactor=sum(DSamp[:,id])/float(sum(DSamp[id,:]))
        times[files[id]]=evolutionFactor
        if not QUIET:
            print('%s %.2f' %(files[id],evolutionFactor))
    return times

def find_DSamp_dir(DSamp): #returns similar DSamp in which DSamp(i,j)==0 if DSamp(i,j)>DSamp(j,i)
    bools=DSamp<=np.transpose(DSamp)
    AMSamp_dir=bools.astype(int)
    return DSamp*AMSamp_dir

def one_way_outdegree(tmp, INPUTS): #returns dictionary with each file as key and outdegree as value
    DSamp=find_DSamp_dir(tmp)
    d={}
    for id in range(len(DSamp)):
        row=DSamp[id]
        input=INPUTS[id]
        c=0
        for item in row:
            if item!=0 and item !=2000:
                c+=1
        d[input]=c
    return d

def check_DSamp(DSamp): #check to see if any values of DSamp are below FAILTIME
    l=len(DSamp)
    for row in range(l):
        for col in range(l):
            if row!=col:
                if DSamp[row,col]!=FAILTIME:
                    return False
    return True

def window(seq, n): # https://docs.python.org/release/2.3.5/lib/itertools-example.html 
    "Returns a sliding window (of width n) over data from the iterable"
    "   s -> (s0,s1,...s[n-1]), (s1,s2,...,sn), ...                   "
    it = iter(seq)
    result = tuple(itertools.islice(it, n))
    if len(result) == n:
        yield result    
    for elem in it:
        result = result[1:] + (elem,)
        yield result

def get_kmers(seqs,kmer_len):
    kmers=defaultdict(int)
    for seq in seqs:
        for item in window(seq,kmer_len):
            kmers[''.join(item)]+=seqs[seq]
    return kmers

def who_won(mat): #compute ranks of column vectors (dense), finds harmonic mean of ranks over value set for each sample
    t=np.transpose(mat)
    tranks=[]
    tranks=np.zeros(np.shape(t))
    for id in range(len(t)):
        row=t[id]
        rankVec=rankdata(row*-1,method='dense')
        tranks[id]=rankVec
    hmens=[]
    for item in np.transpose(tranks):
        hmens.append(hmean(item))
    return hmens
    
def source_from_features(INPUTS,ALI,NOGHOST,QUIET): # determine source from list of features that we calculate
    predictions={}
    if not QUIET:
        print('file,kmerEntropy,nucDiv,nucEntropy,maxHamming,haploFreq')
    for file in INPUTS:
        #some setup, checks
        if ALI:
            with NamedTemporaryFile(delete=False) as aligned:
                alignname=aligned.name
                subprocess.check_call(['mafft', '--quiet', '--auto', '--thread', '20', '--preservecase', file], stdout=aligned) 
                aliName=aligned.name
                seqs=parse_input(aliName)
                haploSize=len(list(seqs.keys())[0])
        else:
            seqs=parse_input(file)
            haploSize=len(list(seqs.keys())[0])
        
        #calculate maximum hamming distance (maxHamming)
        seqs_nofreqs=list(seqs.keys())
        if NOGHOST:
            array=calc_distance_matrix_slow(seqs_nofreqs)
        else:
            array=calc_distance_matrix(seqs_nofreqs)
        maxHamming=np.amax(array)
        
        # calculate haplotide entropy (haploFreq)
        haploFreqs=[]
        totalReads=float(sum(seqs.values()))
        for i in seqs.values():
            freq=i/totalReads
            
            newval=freq*math.log(freq)
            haploFreqs.append(newval)
        haploFreq=sum(haploFreqs)*-1 
        # calculate kmer entropy (kmerEntropy)
        kmers=get_kmers(seqs,8)
        totalKmers=float(sum(kmers.values()))
        newkmers=[]
        for i in kmers.values():
            freq=i/totalKmers
            newval=freq*math.log(freq)
            newkmers.append(newval)
        kmerEntropy=sum(newkmers)*-1 
        #calculate nucleotide entropy (nucEntropy)
        haploNum=len(seqs)
        entropy=calc_ordered_frequencies(haploNum,haploSize,seqs,True)
        nucEntropy=sum(entropy)/len(entropy)
        #calculate nucleotide diversity (nucDiv)
        totalFreq=totalReads
        freqs=[x/totalFreq for x in seqs.values()]
        if NOGHOST:
            mat=calc_distance_matrix_slow(seqs_nofreqs)
        else:
            mat=calc_distance_matrix(seqs_nofreqs)
        nucDiv=0
        for a,b in itertools.combinations(range(haploNum),2):
            nucDiv+=freqs[a]*freqs[b]*mat[a,b]*2/float(haploSize)
        # record result
        plist=[kmerEntropy,nucDiv,nucEntropy,maxHamming,haploFreq]
        predictions[file]=plist
        if not QUIET:
            print(map(str,[file,kmerEntropy,nucDiv,nucEntropy,maxHamming,haploFreq]))
    X=np.zeros([len(INPUTS),len(plist)]) 
    for id in range(len(INPUTS)):
        X[id,:]=predictions[INPUTS[id]]
    hmens=who_won(X)
    retD={}
    for id in range(len(hmens)):
        retD[INPUTS[id]]=1/hmens[id]
    return retD

def max_dict_index(dict):
    #input: dictionary with numerical values
    #output: dictionary key with highest value
    #error if dictionary has 2 keys both having max value (on purpose)
    m=max(dict.values())
    candidates=[]
    for item in dict:
        if dict[item]==m:
            candidates.append(item)
    if len(candidates)>1:
        print(candidates)
        sys.exit('Tie in either evolution simulations or feature calculation! There is a problem somewhere, this output is no good')
    else:
        return candidates[0]

def figure(times,predictions,INPUTS,expect):
    #makes sure that the predicted source values are high enough to warrant a report
    zevo=zscore(list(times.values()))
    zfeatures=zscore(list(predictions.values()))
    # print('features, evolution')
    # print(','.join(map(str,map(max,[zfeatures,zevo]))))
    minz=min(map(max,[zfeatures,zevo]))
    if minz>2:
        return expect, True
    else:
        return expect, False
        
def consensus_seq(seqs):
    if len(seqs)==1:
        return seqs[0]
    else:
        q=len(seqs[0])
        arr=np.zeros([4,q])
        order={'A':0,'T':1,'C':2,'G':3}
        border={0:'A',1:'T',2:'C',3:'G'}
        for seq in seqs:
            for id in range(len(seq)):
                char=seq[id]
                if char in 'ATCG':
                    arr[order[char],id]+=1
        out=[]
        for a in range(len(seqs[0])):
            slice=list(arr[:,a])
            out.append(border[slice.index(max(slice))])
        return ''.join(out)
        
def calc_centroids(nclust,filename,NOGHOST):
    dict=parse_input(filename)
    seqs=list(dict.keys())
    freq=list(dict.values())
    if NOGHOST:
        DM=calc_distance_matrix_slow(seqs)
    else:    
        DM=calc_distance_matrix(seqs)
    nseq = len(seqs)
    centroids = []
    if nseq>nclust:
        Z = fastcluster.linkage(DM,method='weighted')
        T = cluster(Z,nclust)
        for c in range(nclust):
            clust = np.where(T == c)[0]
            seqs_clust = []
            for u in range(np.shape(clust)[0]):
                for v in range(freq[clust[u]]):
                    seqs_clust.append(seqs[clust[u]])
            cons = consensus_seq(seqs_clust)
            if cons not in centroids:
                centroids.append(cons)
        return centroids
    else:
        return seqs
    
def matrix_to_weighted_NX(transnet,fileList,probD): #go from a transNet to a weighted networkx graph
    g=nx.DiGraph()
    numRows=len(transnet)
    numCols=len(transnet[0])
    if type(probD)!=bool:
        for row in range(numRows):
            for col in range(numCols):
                if row!=col:
                    if probD[row,col]>.20:
                        g.add_edge(fileList[row],fileList[col],label='{0:.2f}'.format(probD[row,col]))
    else:
        for row in range(numRows):
            for col in range(numCols):
                if row!=col:
                    if transnet[row,col]!=0:
                        g.add_edge(fileList[row],fileList[col])
    return g

def render_with_graphviz(transnet,fileList,OUTPUT,strength,probD,expect): 
    try:
        shortNames=[]
        for item in fileList:
            shortNames.append(trimfh(item))
        gr=matrix_to_weighted_NX(transnet,shortNames,probD)
    except IndexError:
        gr=matrix_to_weighted_NX(transnet,shortNames,probD)
    with NamedTemporaryFile(delete=False) as tmpdot:
        tmpdotname=tmpdot.name
        nx.drawing.nx_agraph.write_dot(gr,tmpdotname)
    if strength==False:
        lineco=0
        with open(tmpdotname,'r') as f:
            lines=f.readlines()
        l=len(lines)
        with NamedTemporaryFile(delete=False) as outdot: 
            outdotname=outdot.name
            for line in lines:
                lineco+=1
                if lineco!=l:
                    outdot.write(line)
                else:
                    outdot.write("{\n")
                    outdot.write("    Legend [shape=none, margin=0, colorscheme=set23 label=<\n")
                    outdot.write("""    <TABLE BORDER="0" CELLBORDER="0" CELLSPACING="0" CELLPADDING="4">\n""")
                    outdot.write("     <TR>\n")
                    outdot.write("""      <TD>%s is a low-confidence source candidate</TD>\n""" %expect)
                    outdot.write("     </TR>\n    </TABLE>\n   >];\n  }\n} ")
        subprocess.check_call(['dot','-Tpng',outdotname,'-o',OUTPUT+'.png'])
        os.unlink(outdotname)
    elif strength==True:
        lineco=0
        with open(tmpdotname,'r') as f:
            lines=f.readlines()
        l=len(lines)
        with NamedTemporaryFile(delete=False) as outdot:
            outdotname=outdot.name
            for line in lines:
                lineco+=1
                if lineco!=l:
                    outdot.write(line)
                else:
                    outdot.write("{\n")
                    outdot.write("    Legend [shape=none, margin=0, colorscheme=set23 label=<\n")
                    outdot.write("""    <TABLE BORDER="0" CELLBORDER="0" CELLSPACING="0" CELLPADDING="4">\n""")
                    outdot.write("     <TR>\n")
                    outdot.write("""      <TD>Among the samples listed, %s appears to be the most likely source candidate</TD>\n""" %expect )
                    outdot.write("     </TR>\n    </TABLE>\n   >];\n  }\n} ")
        subprocess.check_call(['dot','-Tpng',outdotname,'-o',OUTPUT+'.png'])
        os.unlink(outdotname)
    os.unlink(tmpdotname)

def obj_transnet_phylo_fit(AMtree,tree,DM): #calculates likelihood of given genetic distance given a possible transmission tree using a least-square approach
    rec=0;n=len(DM);am=len(AMtree);E=[]; nPairs=0
    for row in range(am):
        for col in range(am):
            if row>=col and AMtree[row,col]!=0:
                E.append([col,row])
    nedges=len(E)
    for row in range(n):
        for col in range(n):
            if DM[row,col]!=0:
                nPairs+=1
    C=np.zeros([nPairs,nedges])
    d=np.zeros([nPairs]) 
    pathsPairs=[]; p=-1
    for i in range(n):
        for j in range(n):
            if DM[i,j]>0:
                d[p]=1
                p+=1
                sparse_amtree=csgraph_from_dense(AMtree)
                distances=dijkstra(sparse_amtree)
                path=graph_shortest_path(distances,i,j)
                pathsPairs.append(path)
                for v in range(len(path)-1):
                    if path[v]<path[v+1]:
                        edge=[path[v],path[v+1]]
                    else:
                        edge=[path[v+1],path[v]]
                    for ind in range(len(E)):
                        ele=E[ind]
                        if edge==ele:
                            C[p,ind]=1/(float(DM[i,j]))
    deg=np.sum(AMtree,axis=0)
    root_tmp=np.where(deg==2)
    leaves_tmp=np.where(deg==1)
    root=int(root_tmp[0])
    leaves=leaves_tmp[0]
    nleaves=len(leaves)
    nbranches=nleaves-1
    treeLen=nleaves+nbranches
    charVectPaths=np.zeros([nleaves,nedges])
    for leaf in leaves:
        path=graph_shortest_path(distances,root,leaf)
        for v in range(len(path)-1):
            if path[v]<path[v+1]:
                edge=[path[v],path[v+1]]
            else:
                edge=[path[v+1],path[v]]
            for ind in range(len(E)):
                ele=E[ind]
                if edge==ele:
                    charVectPaths[leaf,ind]=1
    Aeq=np.zeros([nbranches,nedges])
    beq=np.zeros([nbranches,1])
    for i in range(nbranches):
        Aeq[i]=charVectPaths[i]-charVectPaths[i+1]
    lsq=lsqlin.lsqlin(C,d,0,None,None,Aeq,beq,0)
    x=lsqlin.cvxopt_to_numpy_matrix(lsq['x'])
    DSamp_tree=np.zeros([n,n])
    pathID=-1; vec1id=-1
    vec1=np.zeros(n*n)
    vec2=np.zeros(n*n)
    for i in range(n):
        for j in range(n):
            vec1id+=1
            vec1[vec1id]=DM[j][i] #supposed to be analagous to vec1=DM(:);
            if DM[i][j]>0:
                pathID+=1
                pathEdges=np.where(C[pathID]>0)
                DSamp_tree[i][j]=np.sum(x[pathEdges])
    vec2id=-1
    for i in range(n):
        for j in range(n):
            vec2id+=1
            vec2[vec2id]=DSamp_tree[j][i]
    ind=np.where(vec1>0)
    vec1=vec1[ind]
    vec2=vec2[ind]
    val=np.corrcoef(vec1,vec2)
    if val.ndim==2:
        v=val[0][1]
        if v>rec:
            rec=v
    return rec

def treeAM(tree): #creates adjacency matrix from tree 
    #AMtree===numLabels,numLabels structure with {0,1} as entries which represents same tree as 'tree'
    numLabels=len(tree)
    n=len(tree)
    AMtree=np.zeros([n,n])
    for i in range(n):
        child1=int(tree[i,0])
        child2=int(tree[i,1])
        if child1+child2>0:
            AMtree[i,child1]=1
            AMtree[i,child2]=1
            AMtree[child1,i]=1
            AMtree[child2,i]=1
    return AMtree

def tree_to_transnet(tree,DSamp,nSamp): #Turns tree into transmission network
    nTreeVert=2*nSamp-1
    TransNet=np.zeros([nSamp,nSamp])
    for i in range(nTreeVert):
        l=int(tree[i,2])
        child1=int(tree[i,0])
        child2=int(tree[i,1])
        if child1+child2!=0:
            l1=int(tree[child1,2])
            l2=int(tree[child2,2])
            if l==l1:
                TransNet[l,l2]=DSamp[l,l2]
            if l==l2:
                TransNet[l,l1]=DSamp[l,l1]
    return TransNet  

def calc_labels_phylo_tree(tree,DM,nSamp,FAILTIME): #calculates values for bottom right portion of tree (internal nodes' parents)
    nVert=len(tree)
    s=np.zeros(nVert-1)
    t=np.zeros(nVert-1)
    e=0
    G=nx.Graph()
    for i in range(nVert): 
        G.add_node(i)
        if tree[i,0]>0:
            s[e]=i
            t[e]=tree[i,0]# +-1's were here
            e=e+1
        if tree[i,1]>0:
            s[e]=i
            t[e]=tree[i,1]# +-1's were here
            e=e+1
    for i in range(len(s)):
        G.add_edge(t[i],s[i])
    backOrder=[]
    for node in nx.dfs_preorder_nodes(G): #does this do the same depth first search as matlab's dfsearch? probably.
        backOrder.append(node)
    order=[]
    for ele in reversed(backOrder):
        order.append(ele)
    internal=[]
    for ele in range(nVert):
        zero=tree[ele,0]>0
        one=tree[ele,1]>0
        internal.append(any([zero,one]))
    w=1
    for i in range(nVert):
        v=int(order[i])
        if internal[v]==True:
            child1=int(tree[v,0])
            child2=int(tree[v,1])
            l1=int(tree[child1,2])
            l2=int(tree[child2,2])
            if DM[l1,l2]<DM[l2,l1]:
                tree[v,2]=l1
            else:
                tree[v,2]=l2
            if DM[l1,l2]==FAILTIME and DM[l2,l1]==FAILTIME:
                w=0
    return(tree,w)

def cluster(D,maxclust):
    Z=D+1
    m = len(Z)+1
    T = np.zeros(m) 
    if m <= maxclust:
        return range(m)
    elif maxclust==1:
        return np.zeros(m)
    else:
        clsnum = 1
        for k in range(m-maxclust+1,m):
            i = int(Z[k-1,0]) # left tree
            if i <= m: # original node, no leafs
                T[i-1] = clsnum
                clsnum = clsnum + 1
            elif i < (2*m-maxclust+1): # created before cutoff, search down the tree
                T = cluster_num(Z, T, i-m, clsnum)
                clsnum = clsnum + 1
            i = int(Z[k-1,1]) # right tree
            if i <= m:  # original node, no leafs
                T[i-1] = clsnum
                clsnum = clsnum + 1
            elif i < (2*m-maxclust+1): # created before cutoff, search down the tree
                T = cluster_num(Z, T, i-m, clsnum)
                clsnum = clsnum + 1
        return T-1

def cluster_num(X,T,k,c): #takes in matrix already incremented by +1 versus matlab 
    m=len(X)+1
    while not np.all(k==0): 
        if isinstance(k,int):
            children=X[k-1,0:2]
        else:
            children=[]
            for item in k:
                ind=int(item)
                children.append(int(X[[ind-1],0]))
                children.append(int(X[[ind-1],1]))
            
        # t=children<=m
        # f=children>m
        klen=0
        for ch in children:
            if ch>m:
                klen+=1
        k=np.zeros(klen) 
        kid=0
        for ch in children:
            child=int(ch)
            if child<=m:
                T[child-1]=c
            else:
                k[kid]=child-m
                kid+=1
    return T
    
def flatten(inarr): #mimics behavior of matlab's out=in(:) with the python call out=flatten(in)
    numRows=len(inarr)
    try:
        numCols=len(inarr[0])
    except TypeError:
        out=inarr
        return out
    out=np.zeros(numRows*numCols)
    it=-1
    for col in range(numCols):
        for row in range(numRows):
            it+=1
            out[it]=inarr[row][col]
    return out

def flatten_horiz(inarr): #is this more robust than calling flatten on the transpose? Maybe.
    numRows=len(inarr)
    try:
        numCols=len(inarr[0])
    except TypeError:
        out=inarr
        return out
    out=np.zeros(numRows*numCols)
    it=-1
    for row in range(numRows):
        for col in range(numCols):
            it+=1
            out[it]=inarr[row][col]
    return out

def estimate_hubs(TransNet): #estimates how many hubs we have in the network (strange pavel code)
    AM=(TransNet+np.transpose(TransNet))>0
    deg=sum(AM)
    deg=sorted(deg,reverse=True)
    transposed=[]
    for item in deg:
        transposed.append([item])
    dists=pdist(transposed)
    linkage_out=fastcluster.linkage(dists,method='single')
    clustered=cluster(linkage_out,2)
    c=clustered[0]
    k=len(np.where(clustered==c))
    return k

def obj_transnet_quad_deg(DM,nhubs): #calculates prior probabilty of a possible transmission tree using an s-metric parameter (strange pavel code)
    n=len(DM)
    AM=np.zeros([n,n])
    for i in range(n):
        for j in range(n):
            if DM[i][j]+DM[j][i]>0:
                AM[i][j]=1
    degseq=np.zeros(n)
    for id in range(len(AM)):
        degseq[id]=np.sum(AM[id])
    val=0
    for i in range(n):
        for j in range(i+1,n):
            if AM[i][j]==1:
                val=val+degseq[i]*degseq[j]
    s=sCore(n,nhubs)
    val=math.exp(-abs(val-s)/s)
    return val

def sCore(n,k): # (strange pavel code)
    n=float(n)
    k=float(k)
    d1=(n-k)/k+k-1
    d2=(n-k)/k+1
    return (n-k)/k*d1+(k-1)*(n-k)/k*d2+(k-1)*d1*d2    
    
def modify_tree(tree,steps,nSamp): #central to the whole altorithm: modifies nonzero internal node children stochastically (bottom left portion) 
    tree[nSamp:len(tree),0:2]+=1
    for i in range(steps):
        u=random.randint(nSamp,nSamp*2-2)
        j=random.randint(0,1)
        jinv=1-j
        v=int(tree[u][j])-1
        newtree=tree
        k=random.randint(0,1)
        w=newtree[u][jinv]
        if newtree[v,k]==0:
            continue
        newtree[u][jinv]=newtree[v][k]
        newtree[v][k]=w
        tree=newtree
    tree[nSamp:len(tree),0:2]-=1
    return tree

def find_DSamp_dir(DSamp): #returns similar DSamp in which DSamp(i,j)==0 if DSamp(i,j)>DSamp(j,i)
    bools=DSamp<=np.transpose(DSamp)
    AMSamp_dir=bools.astype(int)
    return DSamp*AMSamp_dir

def phy_tree(linkage_out,nSamp):
    #custom implementation of matlabs phytree
    #creates an ultrametric phylogenetic tree object
    numLeaves=nSamp
    numBranches=nSamp-1
    numLabels=numLeaves+numBranches
    B=linkage_out[:,0:2]
    B=B.astype(int)
    D=np.zeros([numLabels,1])
    names=[]
    for id in range(numLeaves):
        names.append("Leaf "+str(id))
    for id in range(numBranches):
        names.append('Branch '+str(id))
        D[id+numLeaves]=linkage_out[id,2]
    tmp1=np.squeeze(D[B])
    tmp2=np.squeeze(D[numLeaves+(np.arange(numBranches))*[[1],[1]]])-np.transpose(tmp1)
    t2=flatten(tmp2)
    b=flatten_horiz(B)
    if len(tmp2.shape)==1:
        for ind in range(numLabels-1):
            D[ind]=tmp2[ind]
    elif len(tmp2.shape)==2:
        t2=flatten(tmp2)
        b=flatten_horiz(B)
        for ind in range(numLabels-1):
            id=int(b[ind])
            D[ind]=t2[id]
    else:
        raise ValueError('tmp2 is neither 1 nor 2 dimensional, abort')
    D[-1]=0
    needsReorder=False
    mi=np.zeros(numLabels)
    ma=np.zeros(numLabels)
    sz=np.ones(numLabels)
    for ind in range(numLeaves):
        mi[ind]=ind
        ma[ind]=ind
    for ind in range(numBranches):
        mi[ind+numLeaves]=0
        ma[ind+numLeaves]=0 
    for i in range(numBranches):
        v=B[i,:]
        j=i+numLeaves
        mi[j]=min(mi[v])
        ma[j]=max(ma[v])
        sz[j]=sum(sz[v])
        if ma[j]-mi[j]+1!=sz[j]:
            needsReorder=True
    if needsReorder:
        [B,names]=pretty_order(B,D,names,numBranches,numLeaves,numLabels)
    return(B,names)

def pretty_order(B,D,names,numBranches,numLeaves,numLabels): #reorders the leaf nodes to avoid branch crossings
    #custom implementation of matlabs PRETTYORDER
    L=np.zeros(numLabels)
    X=np.zeros(numLabels)
    for a in range(numLeaves):
        L[a]=1
    for ind in range(numBranches):
        v=B[ind,:]
        L[ind+numLeaves]=sum(L[v])
    for ind in range(numBranches-1,-1,-1):
        v=B[ind,:]
        tmp=D[v].reshape(2,)
        # X[v]=D[v]+X[ind+numLeaves]
        X[v]=tmp+X[ind+numLeaves]
    Li=np.zeros(numLabels)
    Ls=np.zeros(numLabels)
    Ls[-1]=numLeaves
    for ind in range(numBranches-1,-1,-1): 
        v=B[ind,:]
        Ls[v]=Ls[ind+numLeaves]
        Li[v]=Li[ind+numLeaves]
        if X[v[1]]-X[v[0]]>=0:
            Ls[B[ind,0]]=Li[B[ind,0]]+L[B[ind,0]]
            Li[B[ind,1]]=Ls[B[ind,1]]-L[B[ind,1]]
        else:
            Ls[B[ind,1]]=Li[B[ind,1]]+L[B[ind,1]]
            Li[B[ind,0]]=Ls[B[ind,0]]-L[B[ind,0]]
    Ls=Ls.astype(int)
   
    qw=np.arange(numLeaves,numLabels)
    for ind in range(numLeaves,numLabels):
        Ls[ind]=ind+1
    namesFin=[]
    for (a,b) in sorted(zip(Ls,names)):
        namesFin.append(b)
    B=Ls[B]
    B=B-1
    return[B,namesFin]

def graph_shortest_path(dijkstra_output,source,target): #walks through output of dijkstra algorithm to find shortest path in unweighted, undirected graph 
    numNodes=len(dijkstra_output)
    if source<0 or source>=numNodes or target<0 or target>=numNodes:
        raise ValueError("Source and target must be integers within [0,numNodes-1]")
    theList=[source]
    dist=int(dijkstra_output[source,target])
    for nodesVisited in range(dist):
        sourceRow=dijkstra_output[source,:]
        neighbors=np.where(sourceRow==1)[0]
        if len(neighbors)==1:
            source=int(neighbors)
            theList.append(source)
            if source==target:
                return theList
        else:
            possibilities=[]
            for item in neighbors:
                if item==target:
                    theList.append(target)
                    return theList
                possibilities.append(dijkstra_output[target,item])
            tmp=np.argmin(possibilities)
            source=neighbors[tmp]
            if dijkstra_output[source,target]==1:
                theList.append(source)
                theList.append(target)
                return theList
            else:
                theList.append(source)

def calc_record(DSamp_dir,tree,nSamp,FAILTIME,k,DSamp_mcmc):
    newtree,w = calc_labels_phylo_tree(tree,DSamp_mcmc,nSamp,FAILTIME)
    AMtree=treeAM(tree) 
    TransNet=tree_to_transnet(tree,DSamp_mcmc,nSamp)
    obj1=obj_transnet_phylo_fit(AMtree,tree,DSamp_dir)
    obj2=obj_transnet_quad_deg(TransNet,k)
    record=w*obj1*obj2
    return record

def find_transnet_MCMC(DSamp_mcmc,FAILTIME): #this is the second 'wrapper' part of the program - iteratively constructs most likely transmission network
    nSamp=len(DSamp_mcmc)
    nVertTree=2*nSamp-1
    bsum=0
    branches_seen=np.zeros([nVertTree,nVertTree])
    tDict=defaultdict(int)
    DSamp_dir=find_DSamp_dir(DSamp_mcmc)
    DSamp_symm=np.zeros([nSamp,nSamp])
    for u in range (nSamp):
        for v in range(u+1,nSamp):
            if DSamp_mcmc[u,v]<=DSamp_mcmc[v,u]:
                DSamp_symm[u,v]=DSamp_mcmc[u,v]
                DSamp_symm[v,u]=DSamp_mcmc[u,v]
            else:
                DSamp_symm[u,v]=DSamp_mcmc[v,u]
                DSamp_symm[v,u]=DSamp_mcmc[v,u]
    l=fastcluster.linkage(DSamp_symm,method='average') #~10x faster than scipy's linkage
    #linkage function: builds hierarchical clustering on distance matrix 
    [branchLists,nodeNames]=phy_tree(l,nSamp)
    leafID=np.zeros(nSamp)
    tree=np.zeros([nVertTree,3])
    for i in range(nSamp):
        name=nodeNames[i]
        leafID[i]=int(name[5:])
        tree[i,2]=i
    for i in range(nSamp-1):
        for j in range(2):
            if branchLists[i,j]<nSamp:
                tree[nSamp+i,j]=leafID[branchLists[i,j]] # +-1's were here
            else:
                tree[nSamp+i,j]=branchLists[i,j] # +-1's were here
    TransNet=tree_to_transnet(tree,DSamp_mcmc,nSamp)
    k=estimate_hubs(TransNet)
    record=calc_record(DSamp_dir,tree,nSamp,FAILTIME,k,DSamp_mcmc)
    TransNetTrees=[]
    TransNetTrees.append(tree)
    bsum+=1
    delta=.0001;iEdgeMotif=1;nEdgeMotifCurr=2;nIter=4837
    i=-1
    while i<=nIter*2:
        i+=1
        if iEdgeMotif==nIter+1:
            iEdgeMotif=1
            nEdgeMotifCurr-=1
        newtree=modify_tree(tree,nEdgeMotifCurr,nSamp)
        if np.all(newtree==0):
            continue
        val=calc_record(DSamp_dir,tree,nSamp,FAILTIME,k,DSamp_mcmc)
        try: #record is sometimes zero
            prob=float(val)/float(record)
        except ZeroDivisionError:
            prob=0
        if 1<=prob:
            tree=newtree
            if i*2>nIter:
                bsum+=1
                stringTree=tree.tostring()
                tDict[stringTree]+=1
        if val > record+delta: #I have never seen this happen
            record=val
            TransNetTrees.append(newtree)
        if abs(val-record)<delta:
            iscont=False
            nopt=len(TransNetTrees)
            for c in range(nopt):
                if np.all(np.equal(TransNetTrees[c],newtree)):
                    iscont=True
                    break
            if not iscont:
                TransNetTrees.append(newtree)
        iEdgeMotif+=1
    nopt=len(TransNetTrees)
    TransNets=[]
    for i in  range(nopt):
        try:
            tree=TransNetTrees[i]
        except UnboundLocalerror:
            tree=TransNetTrees
        AM_tree=tree_to_transnet(tree,DSamp_mcmc,nSamp)
        isExist=False
        for j in range(len(TransNets)):
            if np.all(np.equal(TransNets[j],AM_tree)):
                isExist=True
                break
        if not isExist:
            TransNets.append(AM_tree)
    for item in tDict:
        tDict[item]=tDict[item]/float(bsum)
    if len(TransNets)==1:
        return TransNets[0],tDict
    else:
        return TransNets,tDict

def determine_components(DSamp): # 
    numFiles=len(DSamp)
    AMSamp=np.eye(numFiles)
    for u in range (numFiles):
        for v in range(u+1,numFiles):
            if DSamp[u,v]<=DSamp[v,u]:
                if DSamp[u,v]<FAILTIME-20: 
                    AMSamp[u,v]=1
            else:
                if DSamp[v,u]<FAILTIME-20: 
                    AMSamp[v,u]=1
    sparseMatrix = csgraph_from_dense(AMSamp)
    return connected_components(sparseMatrix, directed=False, connection='weak', return_labels=True)
    
def print_components(DSamp,INPUTS,S,C): #complain and print (more than two connected components via evolutionary metric inside outbreak) (rare for true outbreaks)
    #QUENTIN was made to be able to handle files which arent related, this is sort of a remnant of that
    for a in range(S):
        print('component number: '+str(a+1))
        runlist=[]
        for i in range(len(C)):
            if C[i]==a:
                runlist.append(INPUTS[i])
        if len(runlist)>1:
            rerun='python '+os.getcwd()+'quentin.py'
            for item in runlist:
                rerun=rerun+" "+item
            print(rerun)
        else:
            print("%s does not appear to be connected to any of the other files" % runlist[0])
    sys.exit("exiting PYCIVO...")        

def render_text_output(INPUTS,OUTPUT,strength,workTime,source):
    with open(OUTPUT+'.txt','w') as f:
        f.write("PYCIVO completed successfully on the following FASTA files!\n")
        for i in INPUTS:
            f.write(i+'\n')
        if strength==True:
            f.write('Among the given samples, %s appears to be the most likely candidate for the source' % source)
        else:
            f.write("PYCIVO analysis: LOW CONFIDENCE PREDICTION! %s is the predicted source. Analysis took %.3f seconds. Thank you for using PYCIVO!" % (source, workTime))
        f.write('\n')

def main(INPUTS,ALI,FAILTIME,PERCENT,NOGHOST,OUTPUT,CLUSTERS,QUIET,INFER_NETWORK): #overall wrapper - wraps the wrappers, checks some things
    if not QUIET:
        print('Welcome to PYCIVO!')
        print('reducing files to '+str(CLUSTERS)+' seqs')
        if ALI:
            print('files will be aligned')
        else:
            print('files should already be aligned')
    options['show_progress'] = False #silences output of lsqlin (very annoying, little information contained)
    startTime = time.time()
    numFiles=len(INPUTS)
    DSamp=np.zeros([numFiles,numFiles],int)
    for input in INPUTS:
        if not input.endswith("fas") and not input.endswith('fasta') and not input.endswith('fa'):
            print(input)
            sys.exit("Warning! The above file may not be in fasta format, this will almost certainly cause an error later. Exiting")
#####for serial processing####################################################################################################################
    if NOGHOST: 
        if not QUIET:
            print('Running evolutionary simulations without ghost...')
        for i1,i2 in itertools.combinations(range(len(INPUTS)),2):
            a=get_evol_times_serial(i1,i2,INPUTS,CLUSTERS,PERCENT,FAILTIME,NOGHOST,ALI)
            if a==False:
                sys.exit('alignment error')
            f=a[0]
            b=a[1]
            DSamp[i1,i2]=f
            DSamp[i2,i1]=b
#####for parallel processing##################################################################################################################
    else:
        if not QUIET:
            print('Running evolutionary simulations with ghost...')
        from melnitz import parallel
        parallelFunc=get_evol_times()
        with parallel() as parallel_executor:
            for (f,b,i1,i2) in parallel_executor.map(parallelFunc,itertools.combinations(range(len(INPUTS)),2)):
                DSamp[i1,i2]=f
                DSamp[i2,i1]=b
##############################################################################################################################################
    if not QUIET:
        print('Evolutionary simulations complete! Here are the times')
    if not QUIET:
        workTime =  time.time() - startTime
        print('so far, we took %.3f seconds' %workTime)
    totalFail=check_DSamp(DSamp)
    if totalFail==True:
        print("None of your samples are to be related to any of your samples, according to me! (samples too genetically distant) Exiting")
        print(DSamp)
        endTime = time.time()
        workTime =  endTime - startTime
        statement_time='Analysis took %.3f seconds' % workTime
        sys.exit(statement_time)
    [S,C]=determine_components(DSamp) #S==number of connected components in the graph (number of outbreaks by PYCIVO metric)        
    #if S==1, we have a normal outbreak with all cases related. Software can automatically remove one guy if hes apart but otherwise we just complain and print
    if S>1:
        if S==2: #2 components, see if we can just remove that one guy and continue
            if sum(C)==1:
                offender=np.where(C==1)[0][0]
            elif sum(C)==len(C)-1:
                offender=np.where(C==0)[0][0]
            else:
                print_components(DSamp,INPUTS,S,C)
            DSamp=np.delete(np.delete(DSamp,offender,axis=1),offender,axis=0)
            print('One of your files, '+INPUTS[offender]+' is all by itself, but the others are still related. I deleted it from the list of inputs, and will continue processing')
            del INPUTS[offender]
        else:
            print_components(DSamp,INPUTS,S,C)
    #we have 3 dictionaries in which keys are INPUTS and values are scalar numerical
    #type(times|predictions.values()[0])===float (no ties)
    #type(outdegree.values()[0])===int (ties ok)
    #make sure they all agree, otherwise error out
    # if not QUIET:
    times=outdegree_source(DSamp,INPUTS,QUIET)
    timeSource=max_dict_index(times)
    
    outdegrees=one_way_outdegree(DSamp,INPUTS)
    maxout=max(outdegrees.values())
    
    predictions=source_from_features(INPUTS,ALI,NOGHOST,QUIET)
    featureSource=max_dict_index(predictions)
    #should we not care about this?
    # fail=True
    # for item in outdegrees:
        # if outdegrees[item]==maxout:
            # if item==timeSource:
                # fail=False
    # if fail==True:
        # sys.exit('Error with evolutionary simulations: highest outdegree != highest evolutionary factor. Exiting') #should we not care about this?
    # print(timeSource,featureSource)
    workTime =  time.time() - startTime
    # print('wtf')
    if timeSource==featureSource:
        expect=timeSource
        source,strength=figure(times,predictions,INPUTS,timeSource)
        if strength==True:
            print('Among the given samples, %s appears to be the most likely candidate for the source' % source)
        else:
            print("PYCIVO analysis: LOW CONFIDENCE PREDICTION! %s is the predicted source. Analysis took %.3f seconds. Thank you for using PYCIVO!" % (source, workTime))
        if INFER_NETWORK:
            transnet,tDict=find_transnet_MCMC(DSamp,FAILTIME)
            probD=np.zeros(np.shape(DSamp))
            treeLen=len(DSamp)*2-1
            print("hi")
            if tDict!={}:
                for a in tDict:
                    tr=np.fromstring(a).reshape(treeLen,3)
                    tn=tree_to_transnet(tr,DSamp,len(DSamp))
                    probD+=tn.astype(bool)*tDict[a]
                render_with_graphviz(transnet,INPUTS,OUTPUT,strength,probD,expect)
            else:
                render_with_graphviz(transnet,INPUTS,OUTPUT,strength,False,expect)
        else:
            workTime =  time.time() - startTime
            render_text_output(INPUTS,OUTPUT,strength,workTime,source)
    else:
        print(featureSource,timeSource)
        sys.exit('PYCIVO analysis: source undeterminable or not present - feature source != evolution source! Analysis took %.3f seconds. Thank you for using PYCIVO!' % workTime)
    # os.system('firefox tmp.png')
    sys.exit('PYCIVO analyis was successful! Analysis took %.3f seconds. Thank you for using PYCIVO!' % workTime)
    
###STATIC VARIABLES###
ALI_DEFAULT=False
FAILTIME_DEFAULT=2000
PERCENT_DEFAULT=100
NOGHOST_DEFAULT=True
CLUSTERS_DEFAULT=100
QUIET_DEFAULT=True
INFER_NETWORK_DEFAULT=False
OUTPUT_DEFAULT="PYCIVO_OUTPUT"
###STATIC VARIABLES###

if __name__ == '__main__':
    import argparse # possible arguments to add: delta, nIter
    parser = argparse.ArgumentParser(description='PYCIVO: predict directionality of HCV infection')
    parser.add_argument('files',
        nargs='+',
        help="List  of  files  to  be  analyzed,  order  of  files  does  not  matter")
    parser.add_argument('-o','--output', 
        type=str, required=False, default=OUTPUT_DEFAULT,
        help="The name of the output .png file")
    parser.add_argument('-a',  '--align', 
        action='store_true',  default=ALI_DEFAULT,
        help="Pass this as an argument to align your files. Highly recommended if they are not already aligned")
    parser.add_argument('-t', '--time', 
        type=int, required=False, default=FAILTIME_DEFAULT,
        help="Time allowed for file1 to evolve into file2")
    parser.add_argument('-p',  '--percent', 
        type=int, required=False, default=PERCENT_DEFAULT,
        help="percent evolution needed for you to say that file1 has evolved into file2")
    parser.add_argument('-n',  '--numClusters', 
        type=int, required=False, default=CLUSTERS_DEFAULT,
        help="number of sequences each file is compressed into")
    parser.add_argument('-s', '--slow', 
        action='store_true', default=NOGHOST_DEFAULT,
        help="Pass this flag to process files without the use of ghost")
    parser.add_argument('-q', '--quiet', 
        action='store_true', default=QUIET_DEFAULT,
        help="Pass this flag to process files without the use of ghost")
    parser.add_argument('-in', '--inferNetwork', 
        action='store_true', default=INFER_NETWORK_DEFAULT,
        help="Pass this flag to process files without the use of ghost")
    
    args = parser.parse_args()
    INPUTS=args.files
    ALI=args.align
    FAILTIME=args.time
    PERCENT=args.percent
    NOGHOST=args.slow
    OUTPUT=args.output
    CLUSTERS=args.numClusters
    QUIET=args.quiet
    INFER_NETWORK=args.inferNetwork
    main(INPUTS,ALI,FAILTIME,PERCENT,NOGHOST,OUTPUT,CLUSTERS,QUIET,INFER_NETWORK)