[glu-genetics] push by bioinformed@gmail.com - Fix bug in convert.from_gfr where forward alleles were being stored as... on 2013-01-02 13:49 GMT

[glu-ge...@googlecode.com]

Revision: 727ec0b03a64
Branch: default
Author: Kevin B. Jacobs <jac...@bioinformed.com>
Date: Wed Jan 2 05:47:29 2013
Log: Fix bug in convert.from_gfr where forward alleles were being
stored as AB
when filestrand!='forward'.
http://code.google.com/p/glu-genetics/source/detail?r=727ec0b03a64

Modified:
/glu/lib/illumina.py
/glu/modules/convert/from_gfr.py

=======================================
--- /glu/lib/illumina.py Wed Oct 31 15:52:00 2012
+++ /glu/lib/illumina.py Wed Jan 2 05:47:29 2013
@@ -539,20 +539,23 @@
yield lname,chrom,loc,strand,(a,b)


-def create_abmap(manifest,genome):
+def create_abmap(manifest,genome=None):
'''
Create mapping from A/B probes to final alleles from the oriented
manifest
entries and updating the genome metadata for each locus.
'''
abmap = {}
for lname,chrom,loc,strand,alleles in manifest:
- genome.merge_locus(lname, chromosome=chrom, location=loc,
strand=strand)
+ if genome:
+ genome.merge_locus(lname, chromosome=chrom, location=loc,
strand=strand)
+
if alleles:
abmap[lname] = alleles
+
return abmap


-def
create_Illumina_abmap(filename,genome,targetstrand='customer',errorhandler=None):
+def
create_Illumina_abmap(filename,genome=None,targetstrand='customer',errorhandler=None):
manifest = IlluminaManifest(filename)
manifest =
orient_manifest(manifest,targetstrand=targetstrand,errorhandler=errorhandler)
return create_abmap(manifest,genome)
=======================================
--- /glu/modules/convert/from_gfr.py Thu Aug 30 13:38:22 2012
+++ /glu/modules/convert/from_gfr.py Wed Jan 2 05:47:29 2013
@@ -263,7 +263,7 @@
return a


-def gdat_writer(gdat, gfr_data, genome, transform, abmap=None):
+def gdat_writer(gdat, gfr_data, genome, transform, abmap=None,
fwdmap=None):
sample_chunk = 16

gdat_SNPs = gdat['SNPs']
@@ -354,10 +354,13 @@
snp_recs = []

for i,name in enumerate(snp_names):
- ab = abmap[name] if abmap is not None else 'AB'
+ # Store forward alleles
+ ab = fwdmap[name] if fwdmap is not None else 'AB'
loc = locusmap[name]
snp_recs.append( (name,loc.chromosome,loc.location,''.join(ab)) )

+ # Remap file alleles to AB
+ ab = abmap[name] if abmap is not None else 'AB'
gmap = mapcache.get(ab)
if gmap is None:
a = fix_allele(ab[0])
@@ -375,6 +378,7 @@
for name,gmap,g in izip(snp_names,genomap,geno):
if g not in gmap:
print 'Invalid allele mapping for locus %s, genotype=%s,
mapping=%s' % (name,g,gmap)
+ raise

gc_chunk[j] = gc
x_chunk[j] = x
@@ -438,12 +442,18 @@
abmap =
create_Illumina_abmap(options.manifest,genome,targetstrand=filestrand,
errorhandler=errorhandler)

+ if filestrand!='forward':
+ fwdmap =
create_Illumina_abmap(options.manifest,genome,targetstrand='forward',
+ errorhandler=errorhandler)
+ else:
+ fwdmap = abmap
+
sys.stderr.write('done.\n')

gdat = create_gdat(options.output, num_snps, num_samples)

#gfr_data = progress_bar(gfr_data,num_samples)
- gdat_writer(gdat, gfr_data, genome, transform, abmap)
+ gdat_writer(gdat, gfr_data, genome, transform, abmap, fwdmap)

attrs = gdat.attrs
attrs['GLU_FORMAT'] = 'gdat'