[glu-genetics] 2 new revisions pushed by bioinformed@gmail.com on 2012-07-25 12:13 GMT
6 views
Skip to first unread message
glu-ge...@googlecode.com
Jul 25, 2012, 8:14:10 AM7/25/12
to glu...@googlegroups.com
2 new revisions:
Revision: d5a41fc190ce
Author: Kevin B. Jacobs <jac...@bioinformed.com>
Date: Wed Jul 25 05:07:46 2012
Log: Major clean-up of normalization. Added storage of cluster
centers and...
http://code.google.com/p/glu-genetics/source/detail?r=d5a41fc190ce
Revision: 90b601b8a94c
Author: Kevin B. Jacobs <jac...@bioinformed.com>
Date: Wed Jul 25 05:12:51 2012
Log: Fix imports, null chromosome checks, naming issues
http://code.google.com/p/glu-genetics/source/detail?r=90b601b8a94c
==============================================================================
Revision: d5a41fc190ce
Author: Kevin B. Jacobs <jac...@bioinformed.com>
Date: Wed Jul 25 05:07:46 2012
Log: Major clean-up of normalization. Added storage of cluster
centers and
(temporarily) disabled compression on re-estimated LRR and BAF. This
progress was mostly inspired by trying to fix a memory leak in SciPy, but
the end result was useful and productive: the leak is fixed and the code
looks much better.
http://code.google.com/p/glu-genetics/source/detail?r=d5a41fc190ce
Modified:
/glu/modules/cnv/gdat.py
/glu/modules/cnv/normalize.py
=======================================
--- /glu/modules/cnv/gdat.py Fri Jul 20 09:54:30 2012
+++ /glu/modules/cnv/gdat.py Wed Jul 25 05:07:46 2012
@@ -225,6 +225,9 @@
def __len__(self):
return self.sample_count
+ def flush(self):
+ self.gdat.flush()
+
def close(self):
self.gdat.close()
@@ -243,21 +246,25 @@
def block_iter(table):
- if not hasattr(table,'chunks'):
- for row in table:
- yield row
- return
-
- chunksize = table.chunks[0]
- start = 0
- last = len(table)
-
- while start<last:
- stop = min(last,start+chunksize)
- chunk = table[start:stop]
- for row in chunk:
+ chunksize = table.chunks[0] if getattr(table,'chunks',0) and
table.chunks else 1
+
+ assert chunksize>=0
+
+ # Fall back to normal iteration
+ if chunksize==1:
+ for row in table:
yield row
- start = stop
+
+ # Explicitly read each chunk to avoid caching issues
+ else:
+ start = 0
+ last = len(table)
+ while start<last:
+ stop = min(last,start+chunksize)
+ chunk = table[start:stop]
+ for row in chunk:
+ yield row
+ start = stop
def parallel_gdat_iter(*tables):
@@ -266,15 +273,15 @@
class BatchTableWriter(object):
- def __init__(self,table):
+ def __init__(self,table,start=0):
self.table = table
- self.batchsize = table.chunks[0]
+ self.batchsize = table.chunks[0] if table.chunks else 1
self.scale = table.attrs['SCALE']
self.nan = table.attrs['NAN']
self.min = table.attrs['MIN']
self.max = table.attrs['MAX']
self.batch = []
- self.index = 0
+ self.index = start
def write(self, value):
batch = self.batch
@@ -291,9 +298,11 @@
table = self.table
+ #print ' Writing chunk %d (size=%d)...' %
(self.index//self.batchsize+1,self.batchsize)
+
if len(batch)==1:
chunk = gdat_encode_f(batch[0], self.scale, self.min,
self.max, self.nan)
- table[self.index] = chunk
+ table[self.index] = np.array(chunk, dtype=table.dtype)
self.index += 1
else:
chunk = [ gdat_encode_f(c, self.scale, self.min,
self.max, self.nan) for c in batch ]
@@ -306,6 +315,9 @@
self.index = end
+ #print ' Flushing chunk %d (size=%d)...' %
(self.index//self.batchsize,self.batchsize)
+
+ table.file.flush()
batch[:] = []
def close(self):
@@ -313,30 +325,48 @@
self.table = None
-def create_gdat_qn(gdatobject,s,n):
- comp = dict(compression='gzip',compression_opts=5)
+def create_gdat_cluster_model(gdatobject):
+ s = gdatobject.snp_count
+ gdat = gdatobject.gdat
+ gdat.require_dataset('CLUSTER_R', (3,s), np.float64, fillvalue=np.nan)
+ gdat.require_dataset('CLUSTER_T', (3,s), np.float64, fillvalue=np.nan)
+
+
+def create_gdat_qn(gdatobject):
+ n = gdatobject.sample_count
+ s = gdatobject.snp_count
+
+ #comp = dict(compression='gzip',compression_opts=5)
chunks = (1,s)
+ comp = {}
+ #chunks = None
shape = (n,s)
shuffle = False
gdat = gdatobject.gdat
+ #BAF_QN = gdat.require_dataset('BAF_QN', shape, BAF_TYPE,
+ #
maxshape=shape,chunks=chunks,shuffle=shuffle,
+ # fillvalue=BAF_NAN,**comp)
BAF_QN = gdat.require_dataset('BAF_QN', shape, BAF_TYPE,
-
maxshape=shape,chunks=chunks,shuffle=shuffle,
- fillvalue=BAF_NAN,**comp)
+
shuffle=shuffle,fillvalue=BAF_NAN,**comp)
BAF_QN.attrs['SCALE'] = BAF_SCALE
BAF_QN.attrs['NAN'] = BAF_NAN
BAF_QN.attrs['MIN'] = BAF_MIN
BAF_QN.attrs['MAX'] = BAF_MAX
+ #LRR_QN = gdat.require_dataset('LRR_QN', shape, LRR_TYPE,
+ #
maxshape=shape,chunks=chunks,shuffle=shuffle,
+ # fillvalue=LRR_NAN,**comp)
LRR_QN = gdat.require_dataset('LRR_QN', shape, LRR_TYPE,
-
maxshape=shape,chunks=chunks,shuffle=shuffle,
- fillvalue=LRR_NAN,**comp)
+
shuffle=shuffle,fillvalue=LRR_NAN,**comp)
LRR_QN.attrs['SCALE'] = LRR_SCALE
LRR_QN.attrs['NAN'] = LRR_NAN
LRR_QN.attrs['MIN'] = LRR_MIN
LRR_QN.attrs['MAX'] = LRR_MAX
+ gdat.flush()
+
def sqlite_magic(con):
con.execute('PRAGMA synchronous=OFF;')
=======================================
--- /glu/modules/cnv/normalize.py Wed Feb 1 08:31:41 2012
+++ /glu/modules/cnv/normalize.py Wed Jul 25 05:07:46 2012
@@ -10,6 +10,7 @@
import os
+import sys
from math import modf
@@ -18,7 +19,8 @@
from glu.lib.glm import Linear
from glu.lib.progressbar import progress_loop
-from glu.modules.cnv.gdat import GDATFile, BatchTableWriter,
parallel_gdat_iter, create_gdat_qn, get_gcmodel
+from glu.modules.cnv.gdat import GDATFile, BatchTableWriter,
parallel_gdat_iter, \
+ create_gdat_qn,
create_gdat_cluster_model, get_gcmodel
from glu.lib.genolib.transform import _intersect_options, _union_options
@@ -295,6 +297,7 @@
#lm0 = Linear(r_geno_fit, design_fit0)
#lm1 = Linear(r_geno_fit, design_fit1)
#lm2 = Linear(r_geno_fit, design_fit2)
+
lm = Linear(r_geno_fit, design_fit)
#lm0.fit()
@@ -332,65 +335,11 @@
return r
-def option_parser():
- from glu.lib.glu_argparse import GLUArgumentParser
-
- parser = GLUArgumentParser(description=__abstract__)
-
- parser.add_argument('gdat', help='Input GDAT
file')
-
- parser.add_argument('--includemodel', metavar='FILE', action='append',
- help='List of samples to use to recalibrate model')
- parser.add_argument('--excludemodel', metavar='FILE', action='append',
- help='List of samples not to use to recalibrate model')
- parser.add_argument('--prenorm', metavar='NAME', default='quantile',
choices=['none','quantile'],
- help='Intensity pre-normalization: none or
quantile (default)')
- parser.add_argument('--rmodel', metavar='NAME', default='quadratic',
choices=['linear','quadratic'],
- help='Intensity correction model: linear or
quadratic (default)')
- parser.add_argument('--gcmodel', metavar='GCM', help='GC model file
to use')
- parser.add_argument('--gcmodeldir', metavar='DIR', help='GC models
directory')
- parser.add_argument('--maxmissing', metavar='RATE', default=0.05,
type=float,
- help='Maximum missing rate for assays to estimate
cluster centers (default=0.05)')
- parser.add_argument('--minqual', metavar='N', default=0.01,
type=float,
- help='Minimum genotype quality score (GC)
(default=0.01)')
- parser.add_argument('--minqqual', metavar='Q', default=0, type=float,
- help='Minimum genotype quality score (GC)
quantile (default=0, disabled)')
- parser.add_argument('-P', '--progress', action='store_true',
- help='Show analysis progress bar, if possible')
-
- return parser
-
-
-def main():
- parser = option_parser()
- options = parser.parse_args()
-
- if options.rmodel=='linear':
- extra_terms = 2
- elif options.rmodel=='quadratic':
- extra_terms = 4
- else:
- raise ValueError('Invalid intensity correction model')
-
- gccorrect = bool(options.gcmodel or options.gcmodeldir)
-
- gdat = GDATFile(options.gdat,'r+')
+def pass1(gdat,options,design,dmask,r_AA,t_AA,r_AB,t_AB,r_BB,t_BB):
n = gdat.sample_count
s = gdat.snp_count
qnorm = options.prenorm=='quantile'
- create_gdat_qn(gdat,s,n)
-
- if gccorrect:
- manifest =
os.path.splitext(os.path.basename(gdat.attrs['ManifestName']))[0]
- print 'Loading GC/CpG model for %s...' % manifest
- filename = options.gcmodel or '%s/%s.gcm' %
(options.gcmodeldir,manifest)
-
- design,dmask = get_gcmodel(filename, gdat.chromosome_index,
ploidy=False, extra_terms=extra_terms)
- else:
- design = np.ones( (s,1+extra_terms), dtype=float )
- dmask = design[:,0]==1
-
X = gdat['X']
Y = gdat['Y']
GC = gdat['GC']
@@ -399,23 +348,31 @@
print 'SNPs=%d, samples=%d' % (s,n)
- np.set_printoptions(linewidth=120,precision=2,suppress=True)
- np.seterr(all='ignore')
+ if not options.force and 'CLUSTER_R' in gdat and 'CLUSTER_T' in gdat:
+ r = gdat['CLUSTER_R']
+
+ r_AA[:] = r[0]
+ r_AB[:] = r[1]
+ r_BB[:] = r[2]
+
+ t = gdat['CLUSTER_T']
+
+ t_AA[:] = t[0]
+ t_AB[:] = t[1]
+ t_BB[:] = t[2]
+
+ print 'PASS 1: Read existing cluster data... (skipped re-estimation)'
+
+ return
print 'PASS 1: Re-estimate cluster centers...'
include = _intersect_options(options.includemodel or [])
- exclude = _union_options(options.excludemodel or [])
-
- n_AA = np.zeros(s, dtype=int )
- r_AA = np.zeros(s, dtype=float)
- t_AA = np.zeros(s, dtype=float)
- n_AB = np.zeros(s, dtype=int )
- r_AB = np.zeros(s, dtype=float)
- t_AB = np.zeros(s, dtype=float)
- n_BB = np.zeros(s, dtype=int )
- r_BB = np.zeros(s, dtype=float)
- t_BB = np.zeros(s, dtype=float)
+ exclude = _union_options(options.excludemodel or [])
+
+ n_AA = np.zeros(s, dtype=int)
+ n_AB = np.zeros(s, dtype=int)
+ n_BB = np.zeros(s, dtype=int)
skipped_missing = 0
skipped_exclude = 0
@@ -482,38 +439,147 @@
print ' ... skipped %d samples (%.2f%%) for missing rate > %.2f%%' %
(skipped_missing,skipped_missing/n*100,options.maxmissing*100)
print ' ... removing %d loci with nonsensical allelic ratio (theta)' %
bad.sum()
+ try:
+ create_gdat_cluster_model(gdat)
+
+ r = gdat['CLUSTER_R']
+
+ r[0] = r_AA
+ r[1] = r_AB
+ r[2] = r_BB
+
+ t = gdat['CLUSTER_T']
+
+ t[0] = t_AA
+ t[1] = t_AB
+ t[2] = t_BB
+
+ finally:
+ gdat.flush()
+
+
+def pass2(gdat,options,design,dmask,r_AA,t_AA,r_AB,t_AB,r_BB,t_BB):
print 'PASS 2: Updating LRR and BAF...'
+ n = gdat.sample_count
+ qnorm = options.prenorm=='quantile'
+
+ X = gdat['X']
+ Y = gdat['Y']
+ samples = gdat['Samples']
+ genotypes = gdat['Genotype']
+
+ create_gdat_qn(gdat)
+
LRR_QN = BatchTableWriter(gdat['LRR_QN'])
BAF_QN = BatchTableWriter(gdat['BAF_QN'])
- pass2 = enumerate(parallel_gdat_iter(samples,X,Y,genotypes))
-
- if options.progress:
- pass2 = progress_loop(pass2, length=n, units='samples',
label='PASS 2: ')
-
- for i,(sample,x,y,genos) in pass2:
- if not options.progress:
- print ' Sample %5d / %d: %s' % (i+1,n,sample)
-
- if qnorm:
- x,y = quantile_normalize(x,y,max_threshold=1.5)
-
- t = (2/np.pi)*np.arctan2(y,x)
-
- r =
regress_intensity(x,y,design,dmask,genos,r_AA,r_AB,r_BB,rmodel=options.rmodel,thin=7)
-
- lrr,baf = compute_lrr_baf(t,r,r_AA,r_AB,r_BB,t_AA,t_AB,t_BB)
-
- mask = np.isfinite(lrr)&np.isfinite(baf)
- print ' ... sigmaLRR=%.2f, sigmaBAFhet=%.3f' %
(lrr[mask].std(),baf[mask&(genos=='AB')].std())
-
- LRR_QN.write(lrr)
- BAF_QN.write(baf)
-
- # Close and flush re-normalize tables
- LRR_QN.close()
- BAF_QN.close()
+ try:
+ pass2 = enumerate(parallel_gdat_iter(samples,X,Y,genotypes))
+
+ if options.progress:
+ pass2 = progress_loop(pass2, length=n, units='samples',
label='PASS 2: ')
+
+ for i,(sample,x,y,genos) in pass2:
+ if not options.progress:
+ print ' Sample %5d / %d: %s' % (i+1,n,sample)
+
+ if qnorm:
+ x,y = quantile_normalize(x,y,max_threshold=1.5)
+
+ t = (2/np.pi)*np.arctan2(y,x)
+
+ r =
regress_intensity(x,y,design,dmask,genos,r_AA,r_AB,r_BB,rmodel=options.rmodel,thin=7)
+
+ lrr,baf = compute_lrr_baf(t,r,r_AA,r_AB,r_BB,t_AA,t_AB,t_BB)
+
+ mask = np.isfinite(lrr)&np.isfinite(baf)
+ print ' ... sigmaLRR=%.2f, sigmaBAFhet=%.3f' %
(lrr[mask].std(),baf[mask&(genos=='AB')].std())
+
+ LRR_QN.write(lrr)
+ BAF_QN.write(baf)
+
+ finally:
+ # Close and flush re-normalize tables
+ LRR_QN.close()
+ BAF_QN.close()
+ sys.stderr.flush()
+ sys.stdout.flush()
+
+
+def option_parser():
+ from glu.lib.glu_argparse import GLUArgumentParser
+
+ parser = GLUArgumentParser(description=__abstract__)
+
+ parser.add_argument('gdat', help='Input GDAT
file')
+
+ parser.add_argument('--includemodel', metavar='FILE', action='append',
+ help='List of samples to use to recalibrate model')
+ parser.add_argument('--excludemodel', metavar='FILE', action='append',
+ help='List of samples not to use to recalibrate model')
+ parser.add_argument('--prenorm', metavar='NAME', default='quantile',
choices=['none','quantile'],
+ help='Intensity pre-normalization: none or
quantile (default)')
+ parser.add_argument('--rmodel', metavar='NAME', default='quadratic',
choices=['linear','quadratic'],
+ help='Intensity correction model: linear or
quadratic (default)')
+ parser.add_argument('--gcmodel', metavar='GCM', help='GC model file
to use')
+ parser.add_argument('--gcmodeldir', metavar='DIR', help='GC models
directory')
+ parser.add_argument('--maxmissing', metavar='RATE', default=0.05,
type=float,
+ help='Maximum missing rate for assays to estimate
cluster centers (default=0.05)')
+ parser.add_argument('--minqual', metavar='N', default=0.01,
type=float,
+ help='Minimum genotype quality score (GC)
(default=0.01)')
+ parser.add_argument('--minqqual', metavar='Q', default=0, type=float,
+ help='Minimum genotype quality score (GC)
quantile (default=0, disabled)')
+ parser.add_argument('-f', '--force', action='store_true',
+ help='Force cluster re-estimation')
+ parser.add_argument('-P', '--progress', action='store_true',
+ help='Show analysis progress bar, if possible')
+
+ return parser
+
+
+def main():
+ parser = option_parser()
+ options = parser.parse_args()
+
+ if options.rmodel=='linear':
+ extra_terms = 2
+ elif options.rmodel=='quadratic':
+ extra_terms = 4
+ else:
+ raise ValueError('Invalid intensity correction model')
+
+ np.set_printoptions(linewidth=120,precision=2,suppress=True)
+ np.seterr(all='ignore')
+
+ gdat = GDATFile(options.gdat,'r+')
+ s = gdat.snp_count
+
+ r_AA = np.zeros(s, dtype=float)
+ t_AA = np.zeros(s, dtype=float)
+ r_AB = np.zeros(s, dtype=float)
+ t_AB = np.zeros(s, dtype=float)
+ r_BB = np.zeros(s, dtype=float)
+ t_BB = np.zeros(s, dtype=float)
+
+ gccorrect = bool(options.gcmodel or options.gcmodeldir)
+
+ if gccorrect:
+ manifest =
os.path.splitext(os.path.basename(gdat.attrs['ManifestName']))[0]
+ print 'Loading GC/CpG model for %s...' % manifest
+ filename = options.gcmodel or '%s/%s.gcm' %
(options.gcmodeldir,manifest)
+
+ design,dmask = get_gcmodel(filename, gdat.chromosome_index,
ploidy=False, extra_terms=extra_terms)
+ else:
+ design = np.ones( (s,1+extra_terms), dtype=float )
+ dmask = design[:,0]==1
+
+ try:
+ pass1(gdat,options,design,dmask,r_AA,t_AA,r_AB,t_AB,r_BB,t_BB)
+ pass2(gdat,options,design,dmask,r_AA,t_AA,r_AB,t_AB,r_BB,t_BB)
+
+ finally:
+ gdat.close()
if __name__ == '__main__':
==============================================================================
Revision: 90b601b8a94c
Author: Kevin B. Jacobs <jac...@bioinformed.com>
Date: Wed Jul 25 05:12:51 2012
Log: Fix imports, null chromosome checks, naming issues
http://code.google.com/p/glu-genetics/source/detail?r=90b601b8a94c
Modified:
/glu/lib/illumina.py
=======================================
--- /glu/lib/illumina.py Fri Jul 20 09:47:41 2012
+++ /glu/lib/illumina.py Wed Jul 25 05:12:51 2012
@@ -11,7 +11,7 @@
import struct
from glu.lib.utils import chunk, namedtuple
-from glu.lib.fileutils import
autofile,parse_augmented_filename,guess_format,get_arg
+from glu.lib.fileutils import
autofile,parse_augmented_filename,guess_format,get_arg,namefile
from glu.lib.sections import read_sections
from glu.lib.genolib.locus import STRAND_UNKNOWN
@@ -238,10 +238,8 @@
snp_count = self.snp_count
norm_ids = self.norm_ids
- unpack4B = struct.Struct('<4B').unpack
unpackL = struct.Struct('<L').unpack
unpackLL = struct.Struct('<LL').unpack
- unpack4L = struct.Struct('<4L').unpack
unpack3BLB = struct.Struct('<3BLB').unpack
unpack4B4L = struct.Struct('<4B4L').unpack
@@ -404,7 +402,7 @@
assayid_idx = find_index(header,['IlmnID','Ilmn ID'])
name_idx = find_index(header,['Name'])
alleles_idx = find_index(header,['SNP'])
- chr_idx = find_index(header,['Chr','CHR'])
+ chrom_idx = find_index(header,['Chr','CHR'])
loc_idx = find_index(header,['MapInfo'])
cstrand_idx = find_index(header,['CustomerStrand','SourceStrand'])
dstrand_idx = find_index(header,['IlmnStrand','Ilmn Strand'])
@@ -416,7 +414,7 @@
for assay in manifest:
lname = assay[name_idx]
alleles = assay[alleles_idx]
- chr = assay[chr_idx] or None
+ chrom = assay[chrom_idx] or None
loc = assay[loc_idx]
cstrand = assay[cstrand_idx].lower()
dstrand = assay[dstrand_idx].lower()
@@ -426,29 +424,29 @@
if targetstrand in ('real_forward','real_reverse'):
gstrand = gstrandmap.get(assay[gstrand_idx]) or 'U'
- if chr.startswith('chr'):
- chr = chr[3:]
- if chr=='Mt':
- chr='M'
+ if chrom and chrom.startswith('chr'):
+ chrom = chrom[3:]
+ if chrom=='Mt':
+ chrom='M'
if cstrand not in validstrand:
errorhandler('Invalid customer strand %s for %s' % (cstrand,lname))
- yield lname,chr,loc,STRAND_UNKNOWN,None
+ yield lname,chrom,loc,STRAND_UNKNOWN,None
continue
if dstrand not in validstrand:
errorhandler('Invalid design strand %s for %s' % (dstrand,lname))
- yield lname,chr,loc,STRAND_UNKNOWN,None
+ yield lname,chrom,loc,STRAND_UNKNOWN,None
continue
if gstrand not in '+-FRU':
errorhandler('Unknown gstrand %s for %s' % (gstrand,lname))
- yield lname,chr,loc,STRAND_UNKNOWN,None
+ yield lname,chrom,loc,STRAND_UNKNOWN,None
continue
if len(alleles) != 5 or (alleles[0],alleles[2],alleles[4]) !=
('[','/',']'):
errorhandler('Invalid SNP alleles %s for %s' % (alleles,lname))
- yield lname,chr,loc,STRAND_UNKNOWN,None
+ yield lname,chrom,loc,STRAND_UNKNOWN,None
continue
a,b = alleles[1],alleles[3]
@@ -458,14 +456,14 @@
# CNV probes can simply be skipped
if (a,b) == NA:
- yield lname,chr,loc,STRAND_UNKNOWN,cnv
+ yield lname,chrom,loc,STRAND_UNKNOWN,cnv
continue
# Indels are strand neutral
elif cstrand in plusminus or dstrand in plusminus:
a = indelmap[a]
b = indelmap[b]
- yield lname,chr,loc,None,(a,b)
+ yield lname,chrom,loc,None,(a,b)
continue
# SNP with A and B alleles on the design strand
@@ -475,7 +473,7 @@
if probea and probeb and (a,b) != (probea[-1],probeb[-1]):
errorhandler('Design alleles do not match probes (%s,%s) !=
(%s,%s) for %s'
% (a,b,probea[-1],probeb[-1],lname))
- yield lname,chr,loc,STRAND_UNKNOWN,None
+ yield lname,chrom,loc,STRAND_UNKNOWN,None
continue
try:
@@ -483,7 +481,7 @@
tstrand = tstrand.lower()
if tstrand != 'top':
errorhandler('Supplied sequence is not correctly normalize to top
strand for %s' % lname)
- yield lname,chr,loc,STRAND_UNKNOWN,None
+ yield lname,chrom,loc,STRAND_UNKNOWN,None
continue
except ValueError:
@@ -495,7 +493,7 @@
# a,b and aa,bb should both be normalized to tstrand and must be equal
if (a,b) != (aa,bb):
errorhandler('Assay alleles do not match sequence alleles
(%s/%s != %s/%s) for %s' % (a,b,aa,bb,lname))
- yield lname,chr,loc,STRAND_UNKNOWN,None
+ yield lname,chrom,loc,STRAND_UNKNOWN,None
continue
if gstrand == '+':
@@ -533,7 +531,7 @@
a,b = complement_base(a),complement_base(b)
strand = strandmap[strand]
- yield lname,chr,loc,strand,(a,b)
+ yield lname,chrom,loc,strand,(a,b)
def create_abmap(manifest,genome):
@@ -542,8 +540,8 @@
entries and updating the genome metadata for each locus.
'''
abmap = {}
- for lname,chr,loc,strand,alleles in manifest:
- genome.merge_locus(lname, chromosome=chr, location=loc, strand=strand)
+ for lname,chrom,loc,strand,alleles in manifest:
+ genome.merge_locus(lname, chromosome=chrom, location=loc,
strand=strand)
if alleles:
abmap[lname] = alleles
return abmap
Revision: d5a41fc190ce
Author: Kevin B. Jacobs <jac...@bioinformed.com>
Date: Wed Jul 25 05:07:46 2012
Log: Major clean-up of normalization. Added storage of cluster
centers and...
http://code.google.com/p/glu-genetics/source/detail?r=d5a41fc190ce
Revision: 90b601b8a94c
Author: Kevin B. Jacobs <jac...@bioinformed.com>
Date: Wed Jul 25 05:12:51 2012
Log: Fix imports, null chromosome checks, naming issues
http://code.google.com/p/glu-genetics/source/detail?r=90b601b8a94c
==============================================================================
Revision: d5a41fc190ce
Author: Kevin B. Jacobs <jac...@bioinformed.com>
Date: Wed Jul 25 05:07:46 2012
Log: Major clean-up of normalization. Added storage of cluster
centers and
(temporarily) disabled compression on re-estimated LRR and BAF. This
progress was mostly inspired by trying to fix a memory leak in SciPy, but
the end result was useful and productive: the leak is fixed and the code
looks much better.
http://code.google.com/p/glu-genetics/source/detail?r=d5a41fc190ce
Modified:
/glu/modules/cnv/gdat.py
/glu/modules/cnv/normalize.py
=======================================
--- /glu/modules/cnv/gdat.py Fri Jul 20 09:54:30 2012
+++ /glu/modules/cnv/gdat.py Wed Jul 25 05:07:46 2012
@@ -225,6 +225,9 @@
def __len__(self):
return self.sample_count
+ def flush(self):
+ self.gdat.flush()
+
def close(self):
self.gdat.close()
@@ -243,21 +246,25 @@
def block_iter(table):
- if not hasattr(table,'chunks'):
- for row in table:
- yield row
- return
-
- chunksize = table.chunks[0]
- start = 0
- last = len(table)
-
- while start<last:
- stop = min(last,start+chunksize)
- chunk = table[start:stop]
- for row in chunk:
+ chunksize = table.chunks[0] if getattr(table,'chunks',0) and
table.chunks else 1
+
+ assert chunksize>=0
+
+ # Fall back to normal iteration
+ if chunksize==1:
+ for row in table:
yield row
- start = stop
+
+ # Explicitly read each chunk to avoid caching issues
+ else:
+ start = 0
+ last = len(table)
+ while start<last:
+ stop = min(last,start+chunksize)
+ chunk = table[start:stop]
+ for row in chunk:
+ yield row
+ start = stop
def parallel_gdat_iter(*tables):
@@ -266,15 +273,15 @@
class BatchTableWriter(object):
- def __init__(self,table):
+ def __init__(self,table,start=0):
self.table = table
- self.batchsize = table.chunks[0]
+ self.batchsize = table.chunks[0] if table.chunks else 1
self.scale = table.attrs['SCALE']
self.nan = table.attrs['NAN']
self.min = table.attrs['MIN']
self.max = table.attrs['MAX']
self.batch = []
- self.index = 0
+ self.index = start
def write(self, value):
batch = self.batch
@@ -291,9 +298,11 @@
table = self.table
+ #print ' Writing chunk %d (size=%d)...' %
(self.index//self.batchsize+1,self.batchsize)
+
if len(batch)==1:
chunk = gdat_encode_f(batch[0], self.scale, self.min,
self.max, self.nan)
- table[self.index] = chunk
+ table[self.index] = np.array(chunk, dtype=table.dtype)
self.index += 1
else:
chunk = [ gdat_encode_f(c, self.scale, self.min,
self.max, self.nan) for c in batch ]
@@ -306,6 +315,9 @@
self.index = end
+ #print ' Flushing chunk %d (size=%d)...' %
(self.index//self.batchsize,self.batchsize)
+
+ table.file.flush()
batch[:] = []
def close(self):
@@ -313,30 +325,48 @@
self.table = None
-def create_gdat_qn(gdatobject,s,n):
- comp = dict(compression='gzip',compression_opts=5)
+def create_gdat_cluster_model(gdatobject):
+ s = gdatobject.snp_count
+ gdat = gdatobject.gdat
+ gdat.require_dataset('CLUSTER_R', (3,s), np.float64, fillvalue=np.nan)
+ gdat.require_dataset('CLUSTER_T', (3,s), np.float64, fillvalue=np.nan)
+
+
+def create_gdat_qn(gdatobject):
+ n = gdatobject.sample_count
+ s = gdatobject.snp_count
+
+ #comp = dict(compression='gzip',compression_opts=5)
chunks = (1,s)
+ comp = {}
+ #chunks = None
shape = (n,s)
shuffle = False
gdat = gdatobject.gdat
+ #BAF_QN = gdat.require_dataset('BAF_QN', shape, BAF_TYPE,
+ #
maxshape=shape,chunks=chunks,shuffle=shuffle,
+ # fillvalue=BAF_NAN,**comp)
BAF_QN = gdat.require_dataset('BAF_QN', shape, BAF_TYPE,
-
maxshape=shape,chunks=chunks,shuffle=shuffle,
- fillvalue=BAF_NAN,**comp)
+
shuffle=shuffle,fillvalue=BAF_NAN,**comp)
BAF_QN.attrs['SCALE'] = BAF_SCALE
BAF_QN.attrs['NAN'] = BAF_NAN
BAF_QN.attrs['MIN'] = BAF_MIN
BAF_QN.attrs['MAX'] = BAF_MAX
+ #LRR_QN = gdat.require_dataset('LRR_QN', shape, LRR_TYPE,
+ #
maxshape=shape,chunks=chunks,shuffle=shuffle,
+ # fillvalue=LRR_NAN,**comp)
LRR_QN = gdat.require_dataset('LRR_QN', shape, LRR_TYPE,
-
maxshape=shape,chunks=chunks,shuffle=shuffle,
- fillvalue=LRR_NAN,**comp)
+
shuffle=shuffle,fillvalue=LRR_NAN,**comp)
LRR_QN.attrs['SCALE'] = LRR_SCALE
LRR_QN.attrs['NAN'] = LRR_NAN
LRR_QN.attrs['MIN'] = LRR_MIN
LRR_QN.attrs['MAX'] = LRR_MAX
+ gdat.flush()
+
def sqlite_magic(con):
con.execute('PRAGMA synchronous=OFF;')
=======================================
--- /glu/modules/cnv/normalize.py Wed Feb 1 08:31:41 2012
+++ /glu/modules/cnv/normalize.py Wed Jul 25 05:07:46 2012
@@ -10,6 +10,7 @@
import os
+import sys
from math import modf
@@ -18,7 +19,8 @@
from glu.lib.glm import Linear
from glu.lib.progressbar import progress_loop
-from glu.modules.cnv.gdat import GDATFile, BatchTableWriter,
parallel_gdat_iter, create_gdat_qn, get_gcmodel
+from glu.modules.cnv.gdat import GDATFile, BatchTableWriter,
parallel_gdat_iter, \
+ create_gdat_qn,
create_gdat_cluster_model, get_gcmodel
from glu.lib.genolib.transform import _intersect_options, _union_options
@@ -295,6 +297,7 @@
#lm0 = Linear(r_geno_fit, design_fit0)
#lm1 = Linear(r_geno_fit, design_fit1)
#lm2 = Linear(r_geno_fit, design_fit2)
+
lm = Linear(r_geno_fit, design_fit)
#lm0.fit()
@@ -332,65 +335,11 @@
return r
-def option_parser():
- from glu.lib.glu_argparse import GLUArgumentParser
-
- parser = GLUArgumentParser(description=__abstract__)
-
- parser.add_argument('gdat', help='Input GDAT
file')
-
- parser.add_argument('--includemodel', metavar='FILE', action='append',
- help='List of samples to use to recalibrate model')
- parser.add_argument('--excludemodel', metavar='FILE', action='append',
- help='List of samples not to use to recalibrate model')
- parser.add_argument('--prenorm', metavar='NAME', default='quantile',
choices=['none','quantile'],
- help='Intensity pre-normalization: none or
quantile (default)')
- parser.add_argument('--rmodel', metavar='NAME', default='quadratic',
choices=['linear','quadratic'],
- help='Intensity correction model: linear or
quadratic (default)')
- parser.add_argument('--gcmodel', metavar='GCM', help='GC model file
to use')
- parser.add_argument('--gcmodeldir', metavar='DIR', help='GC models
directory')
- parser.add_argument('--maxmissing', metavar='RATE', default=0.05,
type=float,
- help='Maximum missing rate for assays to estimate
cluster centers (default=0.05)')
- parser.add_argument('--minqual', metavar='N', default=0.01,
type=float,
- help='Minimum genotype quality score (GC)
(default=0.01)')
- parser.add_argument('--minqqual', metavar='Q', default=0, type=float,
- help='Minimum genotype quality score (GC)
quantile (default=0, disabled)')
- parser.add_argument('-P', '--progress', action='store_true',
- help='Show analysis progress bar, if possible')
-
- return parser
-
-
-def main():
- parser = option_parser()
- options = parser.parse_args()
-
- if options.rmodel=='linear':
- extra_terms = 2
- elif options.rmodel=='quadratic':
- extra_terms = 4
- else:
- raise ValueError('Invalid intensity correction model')
-
- gccorrect = bool(options.gcmodel or options.gcmodeldir)
-
- gdat = GDATFile(options.gdat,'r+')
+def pass1(gdat,options,design,dmask,r_AA,t_AA,r_AB,t_AB,r_BB,t_BB):
n = gdat.sample_count
s = gdat.snp_count
qnorm = options.prenorm=='quantile'
- create_gdat_qn(gdat,s,n)
-
- if gccorrect:
- manifest =
os.path.splitext(os.path.basename(gdat.attrs['ManifestName']))[0]
- print 'Loading GC/CpG model for %s...' % manifest
- filename = options.gcmodel or '%s/%s.gcm' %
(options.gcmodeldir,manifest)
-
- design,dmask = get_gcmodel(filename, gdat.chromosome_index,
ploidy=False, extra_terms=extra_terms)
- else:
- design = np.ones( (s,1+extra_terms), dtype=float )
- dmask = design[:,0]==1
-
X = gdat['X']
Y = gdat['Y']
GC = gdat['GC']
@@ -399,23 +348,31 @@
print 'SNPs=%d, samples=%d' % (s,n)
- np.set_printoptions(linewidth=120,precision=2,suppress=True)
- np.seterr(all='ignore')
+ if not options.force and 'CLUSTER_R' in gdat and 'CLUSTER_T' in gdat:
+ r = gdat['CLUSTER_R']
+
+ r_AA[:] = r[0]
+ r_AB[:] = r[1]
+ r_BB[:] = r[2]
+
+ t = gdat['CLUSTER_T']
+
+ t_AA[:] = t[0]
+ t_AB[:] = t[1]
+ t_BB[:] = t[2]
+
+ print 'PASS 1: Read existing cluster data... (skipped re-estimation)'
+
+ return
print 'PASS 1: Re-estimate cluster centers...'
include = _intersect_options(options.includemodel or [])
- exclude = _union_options(options.excludemodel or [])
-
- n_AA = np.zeros(s, dtype=int )
- r_AA = np.zeros(s, dtype=float)
- t_AA = np.zeros(s, dtype=float)
- n_AB = np.zeros(s, dtype=int )
- r_AB = np.zeros(s, dtype=float)
- t_AB = np.zeros(s, dtype=float)
- n_BB = np.zeros(s, dtype=int )
- r_BB = np.zeros(s, dtype=float)
- t_BB = np.zeros(s, dtype=float)
+ exclude = _union_options(options.excludemodel or [])
+
+ n_AA = np.zeros(s, dtype=int)
+ n_AB = np.zeros(s, dtype=int)
+ n_BB = np.zeros(s, dtype=int)
skipped_missing = 0
skipped_exclude = 0
@@ -482,38 +439,147 @@
print ' ... skipped %d samples (%.2f%%) for missing rate > %.2f%%' %
(skipped_missing,skipped_missing/n*100,options.maxmissing*100)
print ' ... removing %d loci with nonsensical allelic ratio (theta)' %
bad.sum()
+ try:
+ create_gdat_cluster_model(gdat)
+
+ r = gdat['CLUSTER_R']
+
+ r[0] = r_AA
+ r[1] = r_AB
+ r[2] = r_BB
+
+ t = gdat['CLUSTER_T']
+
+ t[0] = t_AA
+ t[1] = t_AB
+ t[2] = t_BB
+
+ finally:
+ gdat.flush()
+
+
+def pass2(gdat,options,design,dmask,r_AA,t_AA,r_AB,t_AB,r_BB,t_BB):
print 'PASS 2: Updating LRR and BAF...'
+ n = gdat.sample_count
+ qnorm = options.prenorm=='quantile'
+
+ X = gdat['X']
+ Y = gdat['Y']
+ samples = gdat['Samples']
+ genotypes = gdat['Genotype']
+
+ create_gdat_qn(gdat)
+
LRR_QN = BatchTableWriter(gdat['LRR_QN'])
BAF_QN = BatchTableWriter(gdat['BAF_QN'])
- pass2 = enumerate(parallel_gdat_iter(samples,X,Y,genotypes))
-
- if options.progress:
- pass2 = progress_loop(pass2, length=n, units='samples',
label='PASS 2: ')
-
- for i,(sample,x,y,genos) in pass2:
- if not options.progress:
- print ' Sample %5d / %d: %s' % (i+1,n,sample)
-
- if qnorm:
- x,y = quantile_normalize(x,y,max_threshold=1.5)
-
- t = (2/np.pi)*np.arctan2(y,x)
-
- r =
regress_intensity(x,y,design,dmask,genos,r_AA,r_AB,r_BB,rmodel=options.rmodel,thin=7)
-
- lrr,baf = compute_lrr_baf(t,r,r_AA,r_AB,r_BB,t_AA,t_AB,t_BB)
-
- mask = np.isfinite(lrr)&np.isfinite(baf)
- print ' ... sigmaLRR=%.2f, sigmaBAFhet=%.3f' %
(lrr[mask].std(),baf[mask&(genos=='AB')].std())
-
- LRR_QN.write(lrr)
- BAF_QN.write(baf)
-
- # Close and flush re-normalize tables
- LRR_QN.close()
- BAF_QN.close()
+ try:
+ pass2 = enumerate(parallel_gdat_iter(samples,X,Y,genotypes))
+
+ if options.progress:
+ pass2 = progress_loop(pass2, length=n, units='samples',
label='PASS 2: ')
+
+ for i,(sample,x,y,genos) in pass2:
+ if not options.progress:
+ print ' Sample %5d / %d: %s' % (i+1,n,sample)
+
+ if qnorm:
+ x,y = quantile_normalize(x,y,max_threshold=1.5)
+
+ t = (2/np.pi)*np.arctan2(y,x)
+
+ r =
regress_intensity(x,y,design,dmask,genos,r_AA,r_AB,r_BB,rmodel=options.rmodel,thin=7)
+
+ lrr,baf = compute_lrr_baf(t,r,r_AA,r_AB,r_BB,t_AA,t_AB,t_BB)
+
+ mask = np.isfinite(lrr)&np.isfinite(baf)
+ print ' ... sigmaLRR=%.2f, sigmaBAFhet=%.3f' %
(lrr[mask].std(),baf[mask&(genos=='AB')].std())
+
+ LRR_QN.write(lrr)
+ BAF_QN.write(baf)
+
+ finally:
+ # Close and flush re-normalize tables
+ LRR_QN.close()
+ BAF_QN.close()
+ sys.stderr.flush()
+ sys.stdout.flush()
+
+
+def option_parser():
+ from glu.lib.glu_argparse import GLUArgumentParser
+
+ parser = GLUArgumentParser(description=__abstract__)
+
+ parser.add_argument('gdat', help='Input GDAT
file')
+
+ parser.add_argument('--includemodel', metavar='FILE', action='append',
+ help='List of samples to use to recalibrate model')
+ parser.add_argument('--excludemodel', metavar='FILE', action='append',
+ help='List of samples not to use to recalibrate model')
+ parser.add_argument('--prenorm', metavar='NAME', default='quantile',
choices=['none','quantile'],
+ help='Intensity pre-normalization: none or
quantile (default)')
+ parser.add_argument('--rmodel', metavar='NAME', default='quadratic',
choices=['linear','quadratic'],
+ help='Intensity correction model: linear or
quadratic (default)')
+ parser.add_argument('--gcmodel', metavar='GCM', help='GC model file
to use')
+ parser.add_argument('--gcmodeldir', metavar='DIR', help='GC models
directory')
+ parser.add_argument('--maxmissing', metavar='RATE', default=0.05,
type=float,
+ help='Maximum missing rate for assays to estimate
cluster centers (default=0.05)')
+ parser.add_argument('--minqual', metavar='N', default=0.01,
type=float,
+ help='Minimum genotype quality score (GC)
(default=0.01)')
+ parser.add_argument('--minqqual', metavar='Q', default=0, type=float,
+ help='Minimum genotype quality score (GC)
quantile (default=0, disabled)')
+ parser.add_argument('-f', '--force', action='store_true',
+ help='Force cluster re-estimation')
+ parser.add_argument('-P', '--progress', action='store_true',
+ help='Show analysis progress bar, if possible')
+
+ return parser
+
+
+def main():
+ parser = option_parser()
+ options = parser.parse_args()
+
+ if options.rmodel=='linear':
+ extra_terms = 2
+ elif options.rmodel=='quadratic':
+ extra_terms = 4
+ else:
+ raise ValueError('Invalid intensity correction model')
+
+ np.set_printoptions(linewidth=120,precision=2,suppress=True)
+ np.seterr(all='ignore')
+
+ gdat = GDATFile(options.gdat,'r+')
+ s = gdat.snp_count
+
+ r_AA = np.zeros(s, dtype=float)
+ t_AA = np.zeros(s, dtype=float)
+ r_AB = np.zeros(s, dtype=float)
+ t_AB = np.zeros(s, dtype=float)
+ r_BB = np.zeros(s, dtype=float)
+ t_BB = np.zeros(s, dtype=float)
+
+ gccorrect = bool(options.gcmodel or options.gcmodeldir)
+
+ if gccorrect:
+ manifest =
os.path.splitext(os.path.basename(gdat.attrs['ManifestName']))[0]
+ print 'Loading GC/CpG model for %s...' % manifest
+ filename = options.gcmodel or '%s/%s.gcm' %
(options.gcmodeldir,manifest)
+
+ design,dmask = get_gcmodel(filename, gdat.chromosome_index,
ploidy=False, extra_terms=extra_terms)
+ else:
+ design = np.ones( (s,1+extra_terms), dtype=float )
+ dmask = design[:,0]==1
+
+ try:
+ pass1(gdat,options,design,dmask,r_AA,t_AA,r_AB,t_AB,r_BB,t_BB)
+ pass2(gdat,options,design,dmask,r_AA,t_AA,r_AB,t_AB,r_BB,t_BB)
+
+ finally:
+ gdat.close()
if __name__ == '__main__':
==============================================================================
Revision: 90b601b8a94c
Author: Kevin B. Jacobs <jac...@bioinformed.com>
Date: Wed Jul 25 05:12:51 2012
Log: Fix imports, null chromosome checks, naming issues
http://code.google.com/p/glu-genetics/source/detail?r=90b601b8a94c
Modified:
/glu/lib/illumina.py
=======================================
--- /glu/lib/illumina.py Fri Jul 20 09:47:41 2012
+++ /glu/lib/illumina.py Wed Jul 25 05:12:51 2012
@@ -11,7 +11,7 @@
import struct
from glu.lib.utils import chunk, namedtuple
-from glu.lib.fileutils import
autofile,parse_augmented_filename,guess_format,get_arg
+from glu.lib.fileutils import
autofile,parse_augmented_filename,guess_format,get_arg,namefile
from glu.lib.sections import read_sections
from glu.lib.genolib.locus import STRAND_UNKNOWN
@@ -238,10 +238,8 @@
snp_count = self.snp_count
norm_ids = self.norm_ids
- unpack4B = struct.Struct('<4B').unpack
unpackL = struct.Struct('<L').unpack
unpackLL = struct.Struct('<LL').unpack
- unpack4L = struct.Struct('<4L').unpack
unpack3BLB = struct.Struct('<3BLB').unpack
unpack4B4L = struct.Struct('<4B4L').unpack
@@ -404,7 +402,7 @@
assayid_idx = find_index(header,['IlmnID','Ilmn ID'])
name_idx = find_index(header,['Name'])
alleles_idx = find_index(header,['SNP'])
- chr_idx = find_index(header,['Chr','CHR'])
+ chrom_idx = find_index(header,['Chr','CHR'])
loc_idx = find_index(header,['MapInfo'])
cstrand_idx = find_index(header,['CustomerStrand','SourceStrand'])
dstrand_idx = find_index(header,['IlmnStrand','Ilmn Strand'])
@@ -416,7 +414,7 @@
for assay in manifest:
lname = assay[name_idx]
alleles = assay[alleles_idx]
- chr = assay[chr_idx] or None
+ chrom = assay[chrom_idx] or None
loc = assay[loc_idx]
cstrand = assay[cstrand_idx].lower()
dstrand = assay[dstrand_idx].lower()
@@ -426,29 +424,29 @@
if targetstrand in ('real_forward','real_reverse'):
gstrand = gstrandmap.get(assay[gstrand_idx]) or 'U'
- if chr.startswith('chr'):
- chr = chr[3:]
- if chr=='Mt':
- chr='M'
+ if chrom and chrom.startswith('chr'):
+ chrom = chrom[3:]
+ if chrom=='Mt':
+ chrom='M'
if cstrand not in validstrand:
errorhandler('Invalid customer strand %s for %s' % (cstrand,lname))
- yield lname,chr,loc,STRAND_UNKNOWN,None
+ yield lname,chrom,loc,STRAND_UNKNOWN,None
continue
if dstrand not in validstrand:
errorhandler('Invalid design strand %s for %s' % (dstrand,lname))
- yield lname,chr,loc,STRAND_UNKNOWN,None
+ yield lname,chrom,loc,STRAND_UNKNOWN,None
continue
if gstrand not in '+-FRU':
errorhandler('Unknown gstrand %s for %s' % (gstrand,lname))
- yield lname,chr,loc,STRAND_UNKNOWN,None
+ yield lname,chrom,loc,STRAND_UNKNOWN,None
continue
if len(alleles) != 5 or (alleles[0],alleles[2],alleles[4]) !=
('[','/',']'):
errorhandler('Invalid SNP alleles %s for %s' % (alleles,lname))
- yield lname,chr,loc,STRAND_UNKNOWN,None
+ yield lname,chrom,loc,STRAND_UNKNOWN,None
continue
a,b = alleles[1],alleles[3]
@@ -458,14 +456,14 @@
# CNV probes can simply be skipped
if (a,b) == NA:
- yield lname,chr,loc,STRAND_UNKNOWN,cnv
+ yield lname,chrom,loc,STRAND_UNKNOWN,cnv
continue
# Indels are strand neutral
elif cstrand in plusminus or dstrand in plusminus:
a = indelmap[a]
b = indelmap[b]
- yield lname,chr,loc,None,(a,b)
+ yield lname,chrom,loc,None,(a,b)
continue
# SNP with A and B alleles on the design strand
@@ -475,7 +473,7 @@
if probea and probeb and (a,b) != (probea[-1],probeb[-1]):
errorhandler('Design alleles do not match probes (%s,%s) !=
(%s,%s) for %s'
% (a,b,probea[-1],probeb[-1],lname))
- yield lname,chr,loc,STRAND_UNKNOWN,None
+ yield lname,chrom,loc,STRAND_UNKNOWN,None
continue
try:
@@ -483,7 +481,7 @@
tstrand = tstrand.lower()
if tstrand != 'top':
errorhandler('Supplied sequence is not correctly normalize to top
strand for %s' % lname)
- yield lname,chr,loc,STRAND_UNKNOWN,None
+ yield lname,chrom,loc,STRAND_UNKNOWN,None
continue
except ValueError:
@@ -495,7 +493,7 @@
# a,b and aa,bb should both be normalized to tstrand and must be equal
if (a,b) != (aa,bb):
errorhandler('Assay alleles do not match sequence alleles
(%s/%s != %s/%s) for %s' % (a,b,aa,bb,lname))
- yield lname,chr,loc,STRAND_UNKNOWN,None
+ yield lname,chrom,loc,STRAND_UNKNOWN,None
continue
if gstrand == '+':
@@ -533,7 +531,7 @@
a,b = complement_base(a),complement_base(b)
strand = strandmap[strand]
- yield lname,chr,loc,strand,(a,b)
+ yield lname,chrom,loc,strand,(a,b)
def create_abmap(manifest,genome):
@@ -542,8 +540,8 @@
entries and updating the genome metadata for each locus.
'''
abmap = {}
- for lname,chr,loc,strand,alleles in manifest:
- genome.merge_locus(lname, chromosome=chr, location=loc, strand=strand)
+ for lname,chrom,loc,strand,alleles in manifest:
+ genome.merge_locus(lname, chromosome=chrom, location=loc,
strand=strand)
if alleles:
abmap[lname] = alleles
return abmap
Reply all
Reply to author
Forward
0 new messages