diff --git a/Sex.DetERRmine.py b/Sex.DetERRmine.py index 28b4348..97c95a2 100755 --- a/Sex.DetERRmine.py +++ b/Sex.DetERRmine.py @@ -3,7 +3,7 @@ from math import sqrt from collections import OrderedDict -VERSION="1.1.1" +VERSION="1.1.2" def CalcErrors(AutSnps, XSnps, YSnps, NrAut, NrX, NrY): SNPs=[AutSnps, XSnps, YSnps] @@ -15,16 +15,21 @@ def CalcErrors(AutSnps, XSnps, YSnps, NrAut, NrX, NrY): rate={} rateErr={} - for Bin,Idx in zip(["Aut", "X", "Y"], range(3)): - Total=sum(Reads) - p[Bin]=Reads[Idx]/Total - ErrNr[Bin]=sqrt(Total*p[Bin]) - dp[Bin]=Reads[Idx]/SNPs[Idx] - Errdp[Bin]=ErrNr[Bin]/SNPs[Idx] + Total=sum(Reads) - for Bin in ["X","Y"]: - rate[Bin]=dp[Bin]/dp["Aut"] - rateErr[Bin]=sqrt((Errdp[Bin]/dp["Aut"])**2 + (Errdp["Aut"]*dp[Bin]/(dp["Aut"]**2))**2) + ## If no sites are covered at all, set error and rate to 0. + if Total == 0: + (rate, rateErr) = ({"Aut":0.0,"X":0.0,"Y":0.0}, {"Aut":0.0,"X":0.0,"Y":0.0}) + else : + for Bin,Idx in zip(["Aut", "X", "Y"], range(3)): + p[Bin]=Reads[Idx]/Total + ErrNr[Bin]=sqrt(Total*p[Bin]) + dp[Bin]=Reads[Idx]/SNPs[Idx] + Errdp[Bin]=ErrNr[Bin]/SNPs[Idx] + + for Bin in ["X","Y"]: + rate[Bin]=dp[Bin]/dp["Aut"] + rateErr[Bin]=sqrt((Errdp[Bin]/dp["Aut"])**2 + (Errdp["Aut"]*dp[Bin]/(dp["Aut"]**2))**2) return (rate, rateErr) @@ -91,6 +96,12 @@ def CalcErrors(AutSnps, XSnps, YSnps, NrAut, NrX, NrY): if Chrom == "X": NrX[Names[x]]+=depths[Names[x]] +if sum([AutSnps, XSnps, YSnps]) == 0: + raise IOError(""" +Input depth file is empty. Stopping execution. + +Please ensure you are using a bed file compatible with your reference genome coordinates.""") + # SortNames=OrderedDict(sorted(Names.items(), key=lambda t: t[1])) print ("#Sample", "#SnpsAut", "#SNPsX", "#SnpsY", "NrAut", "NrX", "NrY", "x-rate", "y-rate", "Err(x-rate)", "Err(y-rate)", sep="\t", file=sys.stdout) data=OrderedDict()