Commit 1caadb43 authored by Motazedi's avatar Motazedi

Haplotype estimation algorithm for trios using NGS reads

parents
TriPoly, version 1.5, Ehsan Motazedi (ehsan.motazedi@gmail.com), last modified: March 28 2018
Introduction
============
**TriPoly** \(executable *main.py*\) is an efficient tool for haplotype estimation in polyploid F1 populations with two parents and one or several offspring, using next generation sequence reads. The method considers each offspring with the parents as a trio and tries to simultaneously maximize the likelihood of all of the trio haplotypes in the population, using a Bayesian framework that incorporates sequence data and recombination probability in the likelihood function. It can handle multi-allelic markers and discontinuous haplotype blocks.
The main parameters of the algorithm: variant calling error, recombination probability between adjacent markers, τ, branching threshold, ρ, and pruning rate, κ, \(*Motazedi et al.* 2017\) can be passed as optional arguments to TriPoly. TriPoly accepts a multi-sample BAM file (assuming the first two sample names are that of the parents and the rest of the progeny), from which the so-called SNP-fragment matrix is constructed for haplotype estimation. To build this matrix from the sequence reads, the minimum acceptable mapping quality and base quality scores, as well as the maixmum acceptable insert-size to consider the reads as paired, could be passed as parameters.
It is important that the parent samples appear before the offspring in the BAM header, as the first two detected sample names in the read group lines of the BAM header are considered parents and the rest are considered as offspring by TriPoly. Also, the VCF file must contain only SNPs, that is to say no complex variants and indels. This can be achieved by running the *break_vcf.sh* script on the original VCF file.
An extra option: -*a*, --*aggressive* has also been provided to specify an aggressive pruning rate, α, so that at each step, (100×α)%; of the candidate extensions are pruned if the original pruning scheme is unable to rule-out at least (100×α)% of these extensions. Setting a high α can greatly increase the speed of the algorithm. Finally, a warm-up length in base-pairs could be set along which, from the start of a haplotype block, no branching or pruning is performed. With --nogen option, it is also possible to estimate the haplotypes from the reads without restricting the dosages to those detected in the VCF file.
By default, TriPoly reports all of the haplotypes that have survived the successive pruning steps, together with their relative likelihoods. The command line option -t, --top changes this behavior to report only the most likely haplotype (or a random haplotype if several haplotypes have a likelihood equal to the maximum).
TriPoly is written in *Python* and is compatible with both Python 2 and Python 3.
Citation:
=====================
To cite TriPoly, please refer to *TriPoly: a haplotype estimation approach for polyploids using sequencing data of related individuals*, Ehsan Motazedi, Dick de Ridder, Richard Finkers and Chris Maliepaard, 2017, bioRxiv 163162; *doi: https://doi.org/10.1101/163162*
Input parameters:
=====================
For the input parameters, see the help of the software:
./main.py --help, -h
***Example:***
./main.py bamfile vcffile output_dirname --kappa 0.85 -a 0.8 --rho 0.2 --top
Download:
=====================
The software is available at gitlab, Wageningen UR:
git clone git@git.wageningenur.nl:motaz001/TriPoly.git —recursive
Copyright and License:
=====================
Copyright (c) 2017, Ehsan Motazedi, Wageningen UR.
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files \(the "Software"\), to deal in the Software without restriction, including without limitation the rights to use, copy, share and modify the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. This permission is only valid for academic and personal use, i.e. not for commercial purposes.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
This diff is collapsed.
File added
This diff is collapsed.
File added
# Break complex variants and throw out indels from the input vcf file. i
# EX: ./break_vcf.sh my.vcf out.vcf
# Written by Ehsan Motazedi, Wageningen UR, 17-10-2017
# Last Modified: 27-07-2018
file=$1
outfile=$2
tmpfile=$(mktemp /tmp/dummy.XXXXXX)
exec 3> "$tmpfile" #open tmpfile for writing on file descpritor 3
python << EOS
import re
with open("$file", 'rU') as vcf, open("$tmpfile", "w") as recoded:
detected_SNPs = dict()
for line in vcf:
if line.lstrip()[0]=='#':
recoded.write(line)
else:
fields = line.split()
contig = fields[0]
if contig not in detected_SNPs.keys():
detected_SNPs[contig] = []
if '.' in fields[3] or '.' in fields[4]: #indels
pass
elif ',' in fields[4]: # multi-allelic variants
_temp_alts = fields[4].split(',')
if not all(len(_alt)==len(fields[3]) for _alt in _temp_alts): # do not parse complex variants that contain indels
pass
else: # break complex variants into bi-allelic SNPs if possible
difflocation = []
diff_ref_alts = []
diff_ref_alts_dic = []
for _base in range(0, len(fields[3])):
ref_alt = [fields[3][_base]] + [_x for _x in set(_alt[_base] for _alt in _temp_alts)-set([fields[3][_base]])] # check at each position within the complex variant if ref, alt corresponds to a SNP
if len(ref_alt)>1:
ref_alt_dic = {ref_alt[_n]:str(_n) for _n in range(0, len(ref_alt))}
ref_alt_dic['.']='.' # missing genotype added
difflocation.append(_base)
diff_ref_alts.append(ref_alt)
diff_ref_alts_dic.append(ref_alt_dic)
for _n in range(0, len(difflocation)):
_pos = str(int(fields[1])+difflocation[_n])
if (_pos not in detected_SNPs[contig]):
recoded.write(fields[0]+'\t'+_pos+'\t'+fields[2]+'\t'+diff_ref_alts[_n][0]+'\t'+','.join(diff_ref_alts[_n][1:])+'\t'+'\t'.join(fields[5:9]))
for _genotype in fields[9:]:
_dosage = re.split('/|\|', _genotype.split(':')[0])
_old_dic = {'0':diff_ref_alts[_n][0], '.':'.'} #'.' stands for missing genotypes
_old_dic.update({str(_m+1):_temp_alts[_m][difflocation[_n]] for _m in range(0, len(_temp_alts))})
_new_dosage = '/'.join(diff_ref_alts_dic[_n][_old_dic[_allele]] for _allele in _dosage)
recoded.write('\t'+_new_dosage+':'+':'.join(_genotype.split(':')[1:]))
recoded.write('\n')
detected_SNPs[contig].append(_pos)
else:
pass
elif len(fields[3])!=len(fields[4]): #indels
pass
elif len(fields[3])==1: #biallelic SNPs
if (fields[1] not in detected_SNPs[contig]):
recoded.write(line)
detected_SNPs[contig].append(fields[1])
else:
pass
else: # MNPs
difflocation = []
for _n in range(0, len(fields[3])):
if fields[3][_n]!=fields[4][_n]:
difflocation.append(_n)
for _n in range(0, len(difflocation)):
_pos = str(int(fields[1])+difflocation[_n])
if (_pos not in detected_SNPs[contig]):
recoded.write(fields[0]+'\t'+_pos+'\t'+fields[2]+'\t'+fields[3][difflocation[_n]]+'\t'+fields[4][difflocation[_n]]+'\t'+'\t'.join(fields[5:])+'\n')
detected_SNPs[contig].append(_pos)
else:
pass
EOS
cp $tmpfile $outfile
status=$?
if [ $status -eq 0 ]; then
rm $tmpfile
fi
This diff is collapsed.
File added
This diff is collapsed.
File added
This diff is collapsed.
# Written by Ehsan Motazedi, Wageningen UR, 05-02-2017.
# Last Updated: 01-09-2017
import math
import sys
cpdef int _xor(x, y):
return 0 if str(x)==str(y) else 1
cpdef Hamming(a, b, null_lst = [".", "-"]):
""" Hamming distance and similarity between two strings or lists or tuples.
(Hamming, Richard W., Error detecting and error correcting codes,
Bell System technical journal 29.2 (1950), pp. 147-160)."""
if isinstance(a, int):
a = [a]
if isinstance(b, int):
b = [b]
if len(a) != len(b):
raise ValueError('The length of the two inputs must be equal for comparison!')
if isinstance(a, str):
a = list(a)
if isinstance(b, str):
b = list(b)
hamvec = [_xor(x, y) for x, y in zip(a,b) if (str(x) not in null_lst) and (str(y) not in null_lst)]
hamgelijk = hamvec.count(0)
hamafstand = hamvec.count(1)
return hamafstand, hamgelijk
cpdef diffVec(a, b, null_lst = [".","-"]):
""" Obtain differences between two lists. Return two vectors: one specifying positions of similarity and one positions of difference."""
if isinstance(a, int):
a = [a]
if isinstance(b, int):
b = [b]
if len(a) != len(b):
raise ValueError('The length of the two inputs must be equal for comparison!')
if isinstance(a, str):
a = list(a)
if isinstance(b, str):
b = list(b)
ongelijk, gelijk = [], []
for n, xy in enumerate(zip(a,b)):
if (str(xy[0]) in null_lst) or (str(xy[1]) in null_lst):
continue
if _xor(*xy)==1:
ongelijk.append(n)
else:
gelijk.append(n)
return ongelijk, gelijk
File added
This diff is collapsed.
# Written by Ehsan Motazedi, Wageningen UR, 06-02-2017.
# Last Updated: 05-07-2017
import sys
from logprob import diffVec, Hamming
from haplotypes import Haplotypes
from libc.math cimport log, exp
from reads import Read
cdef extern from "math.h":
float INFINITY
cpdef double loge(double x):
""" Return log x in base e."""
try:
return log(x)
except ValueError as e:
garbage = sys.stderr.write("WARNING: {0}\n".format(e.args[0]))
return -1*INFINITY
def veclog(vec):
""" Logarithm function defined on vector space."""
try:
return map(loge, vec)
except TypeError:
return loge(vec)
def GetlogProb(r, Vset, eps, qual = None, counts = False, min_r_length=2):
""" Return the log probability of a read conditional on a presumed haplotype and
variant error rate, i.e. P(r|Vset, eps) (Berger et al. 2014, p. 4). If counts if True,
assign the read to the most compatible homologue (count 1 for that homologue and zero of
the others) to later calculate the number of reads compatible with each homologue for a set
of reads by the upstream functions."""
ar = r.GetalleleSet(Vset.GetStart(), Vset.GetStop())
if len([_allele for _allele in ar if _allele not in set(['-','.'])])<min_r_length:
if counts:
return [0, [0 for _h in range(0, len(Vset.GetVS()))]]
return 0 # Corresponding to P[r|Vset, eps] = 1
sr='-'*(r.GetBegin()-Vset.GetStart())+ar+'-'*(Vset.GetStop()-r.GetEnd()) # Make the lengths equal to calculate the distance
#warn, veclogprobs = False, []
warn, probs = False, []
if qual:
try:
er = []
escores = []
for _pos in r.GetPos():
escores.append(10**(-float(qual[_pos])/10))
_n = 0
for _x in ar:
if _x=='-':
er.append('-')
else:
er.append(escores[_n])
_n+=1
er1 = []
er2 = []
for _x in range(0, r.GetBegin()-Vset.GetStart()):
er1.append('-')
for _x in range(0, Vset.GetStop()-r.GetEnd()):
er2.append('-')
er = er1 + er
er = er + er2
except ValueError as e:
e.args=(e.args[0]+"\nInvalid quality scores detected in the alignment file!",)+e.args[1:]
raise
except KeyError as e:
warn = True
qual = None
#garbage=sys.stderr.write("******\nER:\n "+str(er)+"\n")
#garbage=sys.stderr.write("\nSR:\n "+str(sr)+"\n")
if not qual:
for _homolo in Vset.GetVS(): # _homolo corresponds to v in Vset in Berger et al. 2014 p. 4
D, A = Hamming(sr, _homolo) # D corresponds to D(r,v) and A correspond to A(r,v) in Berger et al. 2014 p. 4
#veclogprobs.append(A*veclog((1-eps)/(1-2./3*eps))+D*veclog(eps/3.)) # in case of error, any of the three other bases could have been called. So the probability for a specific erroneous base-call is eps./3. As the probability of no error is 1-eps, we have to normalize by 1-eps+eps/3 to get probs for dichotomy no-error/specific erroneous base-call.
probs.append(((1-eps)/(1-2./3*eps))**A*(eps/3.)**D) # in case of error, any of the three other bases could have been called. So the probability for a specific erroneous base-call is eps./3. As the probability of no error is 1-eps, we have to normalize by 1-eps+eps/3 to get probs for dichotomy no-error/specific erroneous base-call.
else: # Use eps(s) as a function of SNP position in the read instead of a fixed eps.
for _homolo in Vset.GetVS():
D, A = diffVec(sr, _homolo)
#garbage=sys.stderr.write("\nD:\n "+str(D)+"\n")
#garbage=sys.stderr.write("\nA:\n "+str(A)+"\n")
#_veclogprob = 0
_prob = 1
for _s in A:
#_veclogprob+=loge((1-er[_s])/(1-2./3*er[_s]))
_prob*= ((1-er[_s])/(1-2./3*er[_s]))
for _s in D:
#_veclogprob+=loge(er[_s]/3.)
_prob *= (er[_s]/3.)
#veclogprobs.append(_veclogprob)
probs.append(_prob)
#garbage=sys.stderr.write("******\n")
if warn:
garbage = sys.stderr.write("WARNING: Quality scores specified for none or only some positions of the read! The fixed given or deafult error rate was used for all the positions!\n")
unique_homolos = set(Vset.GetVS())
dict_homolos = {_dh:1./len([_h for _h in Vset.GetVS() if _h == _dh]) for _dh in unique_homolos} # the counts are set equally to 1/j for j similar homologues
h_counts=[0 for _h in Vset.GetVS()]
for _k, _h in enumerate(Vset.GetVS()):
#if _h==Vset.GetVS()[veclogprobs.index(min(veclogprobs))]:
if _h==Vset.GetVS()[probs.index(min(probs))]:
h_counts[_k]=dict_homolos[_h]
if counts:
#print([loge(sum(probs)/float(len(Vset.GetVS()))), h_counts])# corresponds to log(P[r|Vset, eps]) = log(1/k * sum(P(r|v in Vset))) in Berger et al. 2014 p. 4, assign the read to the most compatible homologue if counts of the reads assigned to each homologue are needed.
return [loge(sum(probs)/float(len(Vset.GetVS()))), h_counts] # corresponds to log(P[r|Vset, eps]) = log(1/k * sum(P(r|v in Vset))) in Berger et al. 2014 p. 4, assign the read to the most compatible homologue if counts of the reads assigned to each homologue are needed.
#return [loge(sum(map(exp,veclogprobs)))-loge(len(Vset.GetVS())), h_counts] # corresponds to log(P[r|Vset, eps]) = log(1/k * sum(P(r|v in Vset))) in Berger et al. 2014 p. 4, assign the read to the most compatible homologue if counts of the reads assigned to each homologue are needed.
return loge(sum(probs)/float(len(Vset.GetVS()))) # corresponds to log(P[r|Vset, eps]) = log(1/k * sum(P(r|v in Vset))) in Berger et al. 2014 p. 4
#return veclog(sum(map(exp,veclogprobs)))-veclog(len(Vset.GetVS())) # corresponds to veclog(P[r|Vset, eps]) = veclog(1/k * sum(P(r|v in Vset))) in Berger et al. 2014 p. 4
File added
This diff is collapsed.
# Written by Ehsan Motazedi, Wageningen UR, 27-07-2016.
# Last Updated: 22-12-2017
import ast
import networkx as nx
def remove_last_SNP(r):
""" Remove the last SNP from a read."""
d = r.GetDict() #ast.literal_eval(repr(r))
try:
del d[max(d.keys())]
except ValueError:
return Read(d)
return Read(d)
def Extract_Components(fragmentlist):
""" Make the connected Read_Graphs from the list of Read objects, so that each connected Read_Graph is phased separately."""
Read_Graph = nx.Graph() # The SNP-fragment graph to be grown from the input fragment file
_fragdict, _edges, _ls = (dict(), [], []) # Dictionary and lists to add each fragment to the SNP-fragment graph
for _frag in fragmentlist: # Add the fragments to the graph nodes, and add an edge between every two variants that lie on the same fragment
_fragdict = _frag.GetDict() #dict(ast.literal_eval(repr(_frag))) # Convert each fragment Read object to a conventional dictionary
Read_Graph.add_nodes_from(_fragdict.keys()) # Adding the fragments as nodes to the graph
_ls = list((x,y) for x in _fragdict.keys() for y in _fragdict.keys()) # Considering an edge between each two variants within the same fragment
_edges = list(edge for edge in _ls if edge[0]!=edge[1]) # Adding the edges to the graph, avoiding self-loops
Read_Graph.add_edges_from(_edges)
_edges, _ls, _fragdict = (None, None, None)
return sorted(nx.connected_components(Read_Graph), key=lambda x: min(x))
def getReads(fragfile, type=0):
"""Extract the reads from a HapTree fragment file containing a dictionary for each (paired-end) fragment (type=0),\n
as produced by 'fragmentpoly.py' or a HapCut (Bansal & Bafina 2008) format file containing a tab-delimited row for each\n
fragment (type=1). Return the reads as Read objects in a list."""
Rout=[]
with open(fragfile, 'rU') as ffile:
if type==0:
for _strdict in ffile:
_dict = ast.literal_eval(_strdict)
if len(_dict)>1:
Rout.append(Read(_dict)) # Add the read is it has more than two variants
else:
for _rec in ffile:
_reclst = _rec.split()
del _reclst[1], _reclst[-1] # get rid of the fragment name and quality scores
if len(''.join(_reclst[_x] for _x in range(2, len(_reclst), 2)))>1:
Rout.append(Read([int(_reclst[_x])-1 for _x in range(1, len(_reclst)-1, 2)],# convert 1, 2, ... numbering in HapCut to 0, 1, ...
[_reclst[_x] for _x in range(2, len(_reclst), 2)]))
return Rout
class Read(object):
def __init__(self, *args):
""" Two ways to construct a read object:\n 1) Give a fragment dictionary as input.\n 2) Give a vector input:\n
[start positions (single int or list/tuple: first part, second part(maybe None)),\n
alleles (single string or list/tuple: first part, second part (maybe None))].\nNo input corresponds to an empty dictionary."""
if len(args)==0:
args = [dict()]
if len(args)==1:
if not isinstance(args[0], dict):
raise TypeError("The single input must be a dictionary!")
rdict = args[0]
self.pos = sorted((int(_x) for _x in rdict.keys()), key=lambda p: int(p))
self.alleles = []
else:
self.start, self.alleles = args
if (not isinstance(self.start, int)) and (not all(isinstance(_start, int) for _start in self.start)):
raise TypeError("The given start positions of the (sub)reads must be integers!")
self.pos = []
if isinstance(self.start, int):
self.start = [self.start,]
if isinstance(self.alleles, str):
self.alleles = [self.alleles,]
for _n, _start in enumerate(self.start):
self.pos+=range(_start, _start+len(self.alleles[_n]))
rdict = dict(zip(self.pos, list(''.join(self.alleles))))
self.alleles=[]
try:
self.start = self.pos[0]
except IndexError: # This exception occurs if the position list, self.pos, is empty
self.start = 0
try:
self.end = self.pos[-1]
except IndexError:
self.end = -1
for _pos in range(self.start, self.end+1):
if _pos==self.pos[0]:
self.alleles.append(rdict[_pos])
del(self.pos[0])
else:
self.alleles.append('-')
self.pos = list(range(self.start, self.end+1))
self.alleles = self.alleles
def __eq__(self, other):
return isinstance(other, Read) and self.alleles == other.allele and self.pos == other.pos
def __repr__(self):
rdict = dict()
for _pos, _allele in zip(self.pos, self.alleles):
if _allele!='-':
rdict[_pos] = _allele
return str(rdict)
def __str__(self):
if self.pos==[]:
return 'None'
return "begin={0}, end={1}, alleles={2}".format(self.start, self.end, ''.join(str(_x) for _x in self.alleles))
def __len__(self):
return len([_allele for _allele in self.alleles if _allele not in set(['.','-'])])
def isNULL(self):
return self.GetDict()=={}
#return ast.literal_eval(repr(self))=={}
def GetalleleSet(self, start, stop):
if self.pos==[]:
return ''
return ''.join(str(_x) for _x in self.alleles[min(self.end-self.start+1, max(0,start-self.start)):max(0, min(stop+1-self.start, self.end-self.start+1))])
def GetBegin(self):
if self.pos==[]:
return 'None'
return self.start
def GetEnd(self):
if self.pos==[]:
return 'None'
return self.end
def GetDict(self):
rdict = dict()
for _pos, _allele in zip(self.pos, self.alleles):
if _allele!='-':
rdict[_pos] = _allele
return rdict
def GetPos(self):
return sorted(_pos for _pos in (self.GetDict()).keys())
#return sorted(_pos for _pos in dict(ast.literal_eval(repr(self))).keys())
def SemiRead(R, s):
""" Make a semi-read relative to the SNP position 's' from a read 'R' (Berger et al. 2014 p. 4)."""
if not isinstance(R, Read):
raise TypeError("The argument to SemiRead must be a Read object!")
d = R.GetDict()
if s not in d.keys():
return Read(dict())
rdict = dict()
for _s, _allele in d.iteritems():
if _s<s+1:
rdict[_s]=str(_allele)
if len(rdict)<2:
return Read(dict())
else:
return Read(rdict)
def SemiRead_old(R, s):
""" Make a semi-read relative to the SNP position 's' from a read 'R' (Berger et al. 2014 p. 4)."""
if not isinstance(R, Read):
raise TypeError("The argument to SemiRead must be a Read object!")
if s<R.GetBegin()+1 or R.GetPos()==[]:
return Read(dict())
rdict = dict()
for _n, _allele in enumerate(R.GetalleleSet(R.start, s)):
if _allele!='-':
rdict[R.start+_n]=_allele
if s not in rdict.keys() or len(rdict)<2:
return Read(dict())
else:
return Read(rdict)
def SubReads(Rset, Sset, min_subread_length=2):
"""Obtain the sub-reads of a read set, i.e. those reads used
for phasing of Sset and only those read positions relevant to phasing Sset
(Berger et al. 2014 p. 4)."""
Rout=[]
if not all(isinstance(r, Read) for r in Rset):
raise TypeError("The argument to SubReads must be a list of Read objects!")
for r in Rset:
_Sset=set(_s for _s in Sset)
nrdict = dict()
rdict = r.GetDict() #dict(ast.literal_eval(repr(r)))
for s in rdict.keys():
if s in _Sset:
_Sset.remove(s)
nrdict[s]=rdict[s]
if len(nrdict) >= min_subread_length:
Rout.append(Read(nrdict))
else:
Rout.append(Read(dict()))
return Rout
#diclst=[{1:1, 2:1, 3:1, 4:1}, {3:1, 4:1, 5:0, 6:0}, {4:0, 5:1, 6:1},
#{4:0, 5:1,6:1, 7:0}, {5:0, 6:0, 7:1}, {5:1, 6:1, 7:0}]
#a=Read({14:1, 15:1, 16:0, 21:0, 22:1})
#b=Read({14:1, 15:1, 16:0})
#Read(ast.literal_eval(repr(a)))
#str(a)
#str(b)
#a.GetalleleSet(16,22)
#c=Read(dict())
#nlst=map(lambda rdict: Read(rdict), diclst)
#for n in range(1,8):
# for r in nlst:
# print("{0}->{1}".format(n, repr(SemiRead(r, n))))
#SubReads(nlst, {1,2,3,4,5})
#a='1 HWI-ST163_0386:6:27:3566:6857#P7PEM12p 14 01010 ;:;78'.split()
#b='2 HWI-ST163_0386:6:6:10650:188444#P7PEM12p_MP 5 0100 17 000 ;;7<779'.split()
#del a[1], b[1], a[-1], b[-1]
#Read(a[0], int(a[1]), a[2])
#Read(b[0], (int(b[1]), int(b[3])), (b[2], b[4]))
#Read(b[0], (b[1], b[3]), (b[2], b[4]))
File added
#!/usr/bin/env python
import numpy
from distutils.core import setup
from Cython.Build import cythonize
#from distutils.extension import Extension
#from Cython.Distutils import build_ext
#ext_modules = [
# Extension("c_haplotypes", ["c_haplotypes.pyx"]),
# Extension("c_reads", ["c_reads.pyx"]),
# Extension("c_branchprune", ["c_branchprune.pyx"]),
# Extension("HapTreePop", ["HapTreePop.pyx"])
#]
#setup(
#name = 'HapTreePop',
# ext_modules = ext_modules,
# cmdclass = {'build_ext': build_ext}
#)
setup(
ext_modules = cythonize("logprob2.pyx"),
include_dirs=[numpy.get_include()]
)
#python setup.py build_ext --inplace
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment