From 2f51cacc4bb33658bbb6423df370e651d571d125 Mon Sep 17 00:00:00 2001 From: rlyu <rlyu@svi.edu.au> Date: Sat, 16 Jan 2021 17:48:28 +1100 Subject: [PATCH] distr deps --- README.md | 16 +++++++++++++++- src/sscocaller.nim | 24 ++++++++++++------------ sscocaller.nimble | 4 ++-- 3 files changed, 29 insertions(+), 15 deletions(-) diff --git a/README.md b/README.md index 8eef208..5e84fe2 100644 --- a/README.md +++ b/README.md @@ -43,6 +43,7 @@ Examples ## Setup/installation +### Install with nimble `sscocaller` uses `hts-nim`(https://github.com/brentp/hts-nim) that requires the `hts-lib` library. If you are building the `sscocaller` from source, you would need to install `hts-lib` @@ -58,12 +59,25 @@ ls -lh $HOME/htslib/*.so Then, `sscocaller` can be installed using `nimble` -### Install with nimble + `nimble install https://gitlab.svi.edu.au/biocellgen-public/sscocaller.git` The built binary in $HOME/.nimble/bin/sscocaller +### Using docker + +`sscocaller` is also available in docker image [`svirlyu/sscocaller`] https://hub.docker.com/r/svirlyu/sscocaller + +``` +docker run -it svirlyu/sscocaller + +## execute sscocaller +/usr/bin/sscocaller -h + +``` + + ### Static builds The static bianry can be simply downloaded which works for GNU/Linux type OS: \ No newline at end of file diff --git a/src/sscocaller.nim b/src/sscocaller.nim index 88756df..196f360 100755 --- a/src/sscocaller.nim +++ b/src/sscocaller.nim @@ -13,7 +13,7 @@ import hts import tables import sequtils import math -from rmath import dbinom +import distributions/rmath import streams #echo paramCount(), " ", paramStr(1) @@ -201,7 +201,7 @@ proc sscocaller(argv: var seq[string]): int = proc getTrans(pos1:int64,pos2:int64,cmPmb=0.1): float = if pos2<pos1: quit "Wrong order of snps" - var rec = 1-exp(-float(pos2-pos1)*1e-8*cmPmb) # for autosomes 1*cmPbm cM per Mb + var rec = 1-math.exp(-float(pos2-pos1)*1e-8*cmPmb) # for autosomes 1*cmPbm cM per Mb return rec #maxTotalReads:int proc countAllele(rec:Variant, ibam:Bam, chrom:string, mapq: int, @@ -327,8 +327,8 @@ proc sscocaller(argv: var seq[string]): int = if scSpermSeq[ithSperm].viNodeseq.len==0: echo "first node to ", ithSperm var currentViNode = ViNode() - currentViNode.pathScore[stateRef] = ln(initProb[stateRef])+emissionArray[stateRef] - currentViNode.pathScore[stateAlt] = ln(initProb[stateAlt])+emissionArray[stateAlt] + currentViNode.pathScore[stateRef] = math.ln(initProb[stateRef])+emissionArray[stateRef] + currentViNode.pathScore[stateAlt] = math.ln(initProb[stateAlt])+emissionArray[stateAlt] currentViNode.pathState[stateRef] = stateN currentViNode.pathState[stateAlt] = stateN currentViNode.state = stateN @@ -342,8 +342,8 @@ proc sscocaller(argv: var seq[string]): int = else: #echo "additional node to a cell" let preVNode = scSpermSeq[ithSperm].viNodeseq[high(scSpermSeq[ithSperm].viNodeseq)] - var ltransProb = ln(getTrans(preVNode.pos,int(rec.POS),cmPmb=cmPmb)) - var lnoTransProb = ln(1-getTrans(preVNode.pos,int(rec.POS),cmPmb=cmPmb)) + var ltransProb = math.ln(getTrans(preVNode.pos,int(rec.POS),cmPmb=cmPmb)) + var lnoTransProb = math.ln(1-getTrans(preVNode.pos,int(rec.POS),cmPmb=cmPmb)) var currentViNode = ViNode() # ref/alt -> ref var refTref = preVNode.pathScore[stateRef] + lnoTransProb @@ -431,7 +431,7 @@ proc sscocaller(argv: var seq[string]): int = else: # there is transition to different state: ref(start) to alt(end) now output the segment info posStart = scSpermSeq[ithSperm].viNodeseq[^(i)].pos #leftGapSize = scSpermSeq[ithSperm].viNodeseq[^(i)].pos - scSpermSeq[ithSperm].viNodeseq[^(i+1)].pos - leftGap=[ln(transProb),ln(1-transProb)] + leftGap=[math.ln(transProb),math.ln(1-transProb)] inferProb+=leftGap[0] reverseProb+=leftGap[1] viSegmentInfo.writeLine("ithSperm" & $ithSperm & " " & $posStart & " " & $posEnd & " " & $(inferProb-reverseProb) & " " & $cSNP & " 2") @@ -454,7 +454,7 @@ proc sscocaller(argv: var seq[string]): int = else: ## state transit posStart = scSpermSeq[ithSperm].viNodeseq[^(i)].pos - leftGap=[ln(transProb),ln(1-transProb)] + leftGap=[math.ln(transProb),math.ln(1-transProb)] inferProb += leftGap[0] reverseProb += leftGap[1] viSegmentInfo.writeLine("ithSperm" & $ithSperm & " " & $posStart & " " & $posEnd & " " & $(inferProb-reverseProb) & " " & $cSNP & " 1" ) @@ -485,12 +485,12 @@ proc sscocaller(argv: var seq[string]): int = scSpermSeq[ithSperm].viNodeseq[1].pos, cmPmb=cmPmb) if scSpermSeq[ithSperm].viNodeseq[0].state == stateRef: - inferProb = currentEm[stateRef]+ln(transProb) - reverseProb = currentEm[stateAlt]+ln(1-transProb) + inferProb = currentEm[stateRef]+math.ln(transProb) + reverseProb = currentEm[stateAlt]+math.ln(1-transProb) viSegmentInfo.writeLine("ithSperm" & $ithSperm & " " & $posStart & " " & $posEnd & " " & $(inferProb-reverseProb) & " " & $cSNP & " 1" ) else: - inferProb = currentEm[stateAlt]+ln(transProb) - reverseProb = currentEm[stateRef]+ln(1-transProb) + inferProb = currentEm[stateAlt]+math.ln(transProb) + reverseProb = currentEm[stateRef]+math.ln(1-transProb) viSegmentInfo.writeLine("ithSperm" & $ithSperm & " " & $posStart & " " & $posEnd & " " & $(inferProb-reverseProb) & " " & $cSNP & " 2" ) transitFlag = false outFileTotalCountMtx.setPosition(49) diff --git a/sscocaller.nimble b/sscocaller.nimble index 90dfafa..288932a 100644 --- a/sscocaller.nimble +++ b/sscocaller.nimble @@ -10,8 +10,8 @@ bin = @["sscocaller"] # Dependencies -requires "nim >= 1.4.0","docopt","hts >= 0.3.4" -requires "https://github.com/sdwfrost/distributions.git" +requires "nim >= 1.2.0","docopt","hts >= 0.3.4" +requires "https://github.com/ruqianl/distributions.git" # task test, "run the tests": # exec "nim c -d:useSysAssert -d:useGcAssert --lineDir:on --debuginfo -r tests/test_groups" -- GitLab