@merge( PARAMS["filename_gwas"], "gwas_merged.import" )
def importMergedGWAS( 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.

    This method checks 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
            
    counter = E.Counter()
    
    found = set()
    inf = gzip.open( PARAMS["filename_dbsnp"], "r" )
    halfwidth = PARAMS["gwas_interval_halfwidth"]
    # collecat all snp coordinates for each gwas interval
    intervals = collections.defaultdict( list )
    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]:
            intervals[gwas].append( (contig, pos, snp) )
            
        found.add( snp )

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


    roi_id = 0
    for gwas, snps in intervals.iteritems():
        
        def iter_overlaps( snps, halfwidth ):

            snps.sort()
            last_contig, last_pos, last_snp = snps[0]
            width = halfwidth * 2
            r = [(last_pos,last_snp)]
            for contig, pos, snp in snps[1:]:
                if last_contig != contig or \
                        pos - last_pos > width:
                    yield last_contig, r
                    r = []
                    last_contig = contig
                    last_pos = pos
                r.append( (pos, snp) )
            yield last_contig, r
                
            
        for contig, overlaps in iter_overlaps( snps, halfwidth ):
            pos = [x[0] for x in overlaps]
            pos.sort()
            start = pos[0] - halfwidth
            end = pos[-1] + halfwidth
            genes = []
            xsnps = [x[1] for x in overlaps]
            for snp in xsnps:
                for gwas, ggenes in ids[snp]:
                    if ggenes: 
                        genes.append( ggenes )
            genes = sorted(set(genes))

            roi_id += 1
            tmpf.write( "%s\n" % "\t".join( map(str, 
                                                (roi_id, 
                                                 ",".join(xsnps), 
                                                 gwas,
                                                 contig, 
                                                 ",".join(map(str,pos)), 
                                                 start,
                                                 end,
                                                 ",".join(genes),
                                                 ) )  ))



    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( "gwas import: input=%i, found=%i, %s" % (len(ids), len(found), str(counter)))

    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/importMergedGWAS/original