@merge( PARAMS["filename_gwas"], "gwas.import" )
def importGWAS( infile, outfile ):
    '''import GWAS intervals. 

    GWAS intervals are defined by their marking SNP. Each 
    snp is looked up in dbsnp and extended by ``gwas_interval_halfwidth``.

    The same SNP can be part of multiple gwas studies.

    The script checks not for overlap between annotations.
    '''

    reader = csv.DictReader( open(infile,"rU") )

    ids = collections.defaultdict( list )
    last_gwas = None

    for row in reader:

        snps, gwas, genes = row["snp"], row["class"], row["genes"]
        if gwas == "": gwas = last_gwas

        snps = snps.split(",")
        for snp in snps:
            s = snp.strip()
            if s == "na" or s == "": continue
            genes=re.sub(" ", "", genes.lower())
            ids[s.strip()].append( (gwas,genes) )

        last_gwas = gwas
            
    tmpf = P.getTempFile()
    tmpf.write( "%s\n" % "\t".join( ("roi_id","snp","class","contig","pos","start","end","genes")) )

    counter = E.Counter()
    
    found = set()
    inf = gzip.open( PARAMS["filename_dbsnp"], "r" )
    nsnps = 0
    halfwidth = PARAMS["gwas_interval_halfwidth"]
    for line in inf:
        data = line[:-1].split("\t")
        contig, snp, pos = data[1], data[4], data[2]
        if snp not in ids: continue
        if contig == None or contig == "": 
            counter.unmapped += 1
            continue
            
        counter.accepted += 1
        pos = int(pos)
        for gwas, genes in ids[snp]:
            nsnps += 1
            tmpf.write( "%s\n" % "\t".join( map(str, 
                                                (nsnps, 
                                                 snp, 
                                                 gwas,
                                                 contig, 
                                                 pos, 
                                                 pos - halfwidth,
                                                 pos + halfwidth,
                                                 genes
                                                 ) )  ))

            
        found.add( snp )

    inf.close()
    tmpf.close()

    table = outfile[:-len(".import")]
    
    tmpfilename = tmpf.name
    statement = '''
   python %(scriptsdir)s/csv2db.py %(csv2db_options)s \
               --index=roi_id \
               --index=class \
               --map=roi_id:str \
               --table=%(table)s \
    < %(tmpfilename)s > %(outfile)s
    '''
    P.run( **locals() )
    
    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( (len(ids), len(found), str(counter)) % "gwas import: input=%i, found=%i, %s")

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

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