kgori

import scipy
import scipy.stats as ss
import numpy as np

def discretize(alpha, ncat, dist=ss.gamma):
    if dist == ss.gamma:
        dist = dist(alpha, scale=1 / alpha)
    elif dist == ss.lognorm:
        dist = dist(s=alpha, scale=np.exp(0.5 * alpha**2))
    quantiles = dist.ppf(np.arange(0, ncat) / ncat)

kgori

# Replace these with the full path to the executables, if necessary
BWA=bwa
SAMTOOLS=samtools
BCFTOOLS=bcftools

# Makefile targets
all: ref.fa reads.fq aln.bam aln.bam.bai pileup.vcf filtered_calls.bcf consensus.fa
.PHONY: all

simulation: ref.fa reads.fq

kgori

# Standardisation formula
z <- function(x) (x - mean(x)) / sd(x)

#' Convert regression parameters obtained from standardised variables
#' back to the scale of the original data.
#'
#' For a regression y ~ x*, where x* denotes the standardised form of x:
#' Result:      y = m(x*) + c
#' Transformed: y = (m / sd(x))(x) + (c - (m*mean(x)/sd(x))
#'

kgori

def revcomp(s):
    complements = dict(
        A='T', C='G', G='C', T='A',
        M='K', R='Y', W='W', S='S',
        Y='R', K='M', V='B', H='D',
        D='H', B='V', N='N'
    )
    return ''.join(complements[char] for char in s.upper()[::-1])

kgori

import pysam
import math
from itertools import islice

def filterfn(read):
    return (read.is_proper_pair and
            read.is_paired and
            read.tlen > 0 and
            not read.is_supplementary and
            not read.is_duplicate and

kgori

import sys

class Progress():
    def __init__(self, capacity, label = '', width=80, fill='#', blank=' '):
        self.capacity = capacity
        self.label = label
        self.width = width
        self.fill = fill
        self.blank = blank
        self.completed = 0

kgori

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "discrete_gamma.h"
// PAML code for discrete gamma

// Forward decs
double QuantileChi2(double prob, double v);
double LnGammaFunction(double alpha);
double IncompleteGamma(double x, double alpha, double ln_gamma_alpha);

kgori

"""
PYMC3 model code to do signature fitting using NUTS
"""
import pymc3
from theano import tensor

mytest = pymc3.Model()
counts = data
with mytest:
    mixing_proportions = pymc3.Dirichlet('mixing_proportions', a = np.ones(30), shape = 30)

kgori

import re

# Fasta one-liners
# Parse fasta from a string
def fasta_parse_string(s):
    return dict([(lambda x:(x[0].strip(), re.sub(r'\s+','',x[2])))(l.partition('\n')) for l in s.split('>')[1:]])

# Parse fasta from a file
def fasta_parse_file(f):
    with open(f) as h:return fasta_parse_string(h.read())

kgori