import pandas as pd
from matplotlib import pyplot as plt
%matplotlib inline

"""
Explanation: Detecting horizontal gene transfer (HGT)
This sheet shows how we flag genes as potential HGT events.
End of explanation
"""


lossTaxa = pd.read_csv("lossTaxa_HUMAN.csv",index_col=0)
lossTaxa.head()

"""
Explanation: Number of taxa inferred to have been lost for each gene by each algorithm
End of explanation
"""


# Map genes to trimmed algorithms

outlierD = {}
with open("lossStats_HUMAN.csv") as f:
    f.readline() # skip header
    for line in f:
        line = line.strip().split(",")
        if line[3] == '':
            continue
        outlierD[line[0]] = line[3].split()
        
outlierD[outlierD.keys()[0]]

"""
Explanation: Get algorithms that have been removed as false positives
End of explanation
"""


dbsTrimmed = 0
avgs = pd.Series()
for index,row in lossTaxa.iterrows():
    if index in outlierD:
        dbsTrimmed += len(outlierD[index])
        dbs = [i for i in lossTaxa.columns if i not in outlierD[index]]
    else:
        dbs = lossTaxa.columns
    avgs[index] = row[dbs].mean() # only include algorithm that have not been trimmed
    
dbsTrimmed

lossTaxa["Avg"] = avgs
lossTaxa.head()

"""
Explanation: Calculate average number of lossTaxa for each gene and add this as a new column to lossTaxa
End of explanation
"""


def floatRange(start,stop,step):
    i = start
    while i <= stop:
        yield i
        i += step
        
quantile_steps = [i for i in floatRange(0,1,.05)]
    
quantiles = lossTaxa["Avg"].quantile(quantile_steps)
quantiles

quantiles.iloc[-1]

lossTaxa["Avg"].hist(bins=50,color='grey')
bline = plt.axvline(24.92,color='black',label="95th percentile")
plt.legend()

#plt.savefig("AvgLossTaxa_distribution.svg")

"""
Explanation: The average has a fat-tailed distribution
We will flag all the genes in the 95th percentile
End of explanation
"""


lossTaxa["HGT_flag"] = lossTaxa["Avg"] >= 24.92
lossTaxa.head()

lossTaxa.to_csv("HGTFlag_HUMAN.csv")

"""
Explanation: Flag genes
End of explanation
"""


