@merge( PARAMS["filename_regions_of_interest"], "snps.import" )
def importSNPs( infile, outfile ):
    '''import snps from :term:`roi` file.'''

    reader = csv.DictReader( open(infile,"rU") )
    ids = {}
    for row in reader:
        
        roi_id, contig, start, end, snps =\
                row["roi_id"], row["contig"], row["start"], row["end"], row["snp"]
        snps = snps.split(",")
        for snp in snps:
            s = snp.strip()
            if s == "na" or s == "": continue
            try:
                ids[snp.strip()] = (roi_id, contig, int(start), int(end))
            except ValueError:
                # ignore empty roi
                continue
            
    tmpf = P.getTempFile()
    tmpf.write("roi_id\tsnp\tcontig\tpos\n" )

    counter = E.Counter()
    
    found = set()
    inf = gzip.open( PARAMS["filename_dbsnp"], "r" )
    for line in inf:
        data = line[:-1].split("\t")
        contig, snp, pos = data[1], data[4], data[2]
        if snp not in ids: continue

        roi_id, roi_contig, roi_start, roi_end = ids[snp]

        if contig != roi_contig or not roi_start <= int(pos) < roi_end:
            counter.out_of_range += 1
            continue

        counter.accepted += 1

        tmpf.write( "%s\t%s\t%s\t%s\n" % (roi_id, snp, contig, pos ) )
        found.add( snp )

    inf.close()

    table = outfile[:-len(".import")]
    
    tmpfilename = tmpf.name
    statement = '''
   python %(scriptsdir)s/csv2db.py %(csv2db_options)s \
               --index=roi_id \
               --map=roi_id:str \
               --table=%(table)s \
    < %(tmpfilename)s > %(outfile)s
    '''
    P.run( **dict( locals().items() + PARAMS.items() ) )
    
    os.unlink( tmpfilename )

    outfile = open( outfile, "a" )
    outfile.write( "input=%i, found=%i, %s\n" % (len(ids), len(found), str(counter)))
    outfile.write( "snps not found: %s\n" % ",".join( set(ids.keys()).difference( found) ) )
    outfile.close()
    
    E.info( "snp import: input=%i, found=%i, %s" % (len(ids), len(found), str(counter)))

    if len(ids) != len(found):
        E.warn( ",".join( set(ids.keys()).difference( found) ) % "snps not found: %s" )

###################################################################
###################################################################
###################################################################
##
###################################################################
dataset/ETHPy150Open CGATOxford/cgat/obsolete/pipeline_vitaminD.py/importSNPs/swapped operands of `%`