Skip to content

Commit

Permalink
IOError when input is empty. Rates and Errors set to 0 when no sites …
Browse files Browse the repository at this point in the history
…are covered.
  • Loading branch information
TCLamnidis committed Jun 11, 2020
1 parent bb8acb5 commit 89ebec1
Showing 1 changed file with 21 additions and 10 deletions.
31 changes: 21 additions & 10 deletions Sex.DetERRmine.py
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand All @@ -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)

Expand Down Expand Up @@ -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()
Expand Down

0 comments on commit 89ebec1

Please sign in to comment.