@follows( importSNPsOfInterest, runMACS )
@transform(  "*.bam", 
             regex( r"(.*).bam"),
             inputs( (r"\1.bam", "snps_of_interest.import",) ), 
             r"\1.snpcoverage" )
def buildSNPCoverage( infiles, outfile ):

    infile, other = infiles
    dbhandle = sqlite3.connect( PARAMS["database"] )
    cc = dbhandle.cursor()

    tmpfilename = P.getTempFilename()

    outf = open(outfile, "w")
    fields = pysam.Pileup.PileupSubstitution._fields
    outf.write( "snp\tnreads\t%s\n" % "\t".join( fields ))

    samfile = pysam.Samfile( infile, "rb" )
    
    for snp, contig, pos, reference in cc.execute( "SELECT snp, contig, pos, reference FROM snps_of_interest"):

        # count reads in neighbourhood
        start = pos - PARAMS["snp_halfwidth"]
        end = pos + PARAMS["snp_halfwidth"] + 1
        nreads = len( list(samfile.fetch( contig, start, end ) ))

        # collect pileup stats
        start = pos
        end = pos + 1
        statement = '''
        samtools pileup -f genome.fa -c -S <(samtools view -h %(infile)s %(contig)s:%(start)i-%(end)i ) 2>/dev/null > %(tmpfilename)s
        '''
        P.run( **dict( locals().items() + PARAMS.items() ) )

        for s in pysam.Pileup.iterate(open( tmpfilename, "r")):
            if s.position == pos:
                if s.reference_base != reference:
                    E.warn((snp, reference, s.reference_base) % "mismatch for snp %s: expected %s, but got %s")
                outf.write( "%s\t%i\t%s\n" % (snp, 
                                              nreads,
                                              "\t".join( map(str,s)) ) )
                break
        else:
            outf.write( "%s\t%i\t%s\n" % (snp, 
                                          nreads,
                                          "\t".join( [""] * len(fields)) ) )

    os.unlink( tmpfilename )

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