@files( ((None, "mouseC57BL.recalls.gz" ), ) )
def recallGenomicSNPs( infiles, outfile ):
    '''validate the genomic SNP calls for C57BL6 using the
    same methodology as for the RNASeq calls.

    This tool requires a file
    
    xgenome.fa/xgenome.fa.fai in the build directory with
    samtools conform contig names (i.e. without the ``chr``
    prefix).
    '''
    filename_bam = "ftp://ftp.sanger.ac.uk/pub/mouse_genomes/current_bams/C57BL.bam"

    dbhandle = sqlite3.connect( PARAMS["database"] )
    cc = dbhandle.cursor()
    
    statement = "SELECT contig, pos, reference, genotype, status, genotypes FROM mouseC57BL_validated"
    
    samfile = pysam.Samfile( filename_bam, "rb")  
    fastafile = pysam.Fastafile( "xgenome.fa" )
    i = samfile.pileup( select = "snpcalls", fastafile = fastafile )
    caller = pysam.SNPCaller( i )

    outf= IOTools.openFile( outfile, "w" )
    outf.write( "\t".join( ("contig", "pos", "ref", 
                           "orig_genotype", "rnaseq_status", "rnaseq_genotypes", 
                           "recall_genotype", "recall_consensus_quality", 
                           "recall_snp_quality", "recall_mapping_quality", "recall_coverage" ) ) + "\n" )

    for contig, pos, ref, genotype, status, genotypes in cc.execute(statement):

        contig = re.sub("chr", "", contig)
        try:
            call = caller.call( contig, pos )
        except __HOLE__, msg:
            E.warn( "could not call %s:%i: msg=%s" % (contig, pos, msg) )

        outf.write( "\t".join( map(str, (contig, pos, ref, genotype, status, genotypes, 
                                         call.genotype,
                                         call.consensus_quality,
                                         call.snp_quality,
                                         call.mapping_quality,
                                         call.coverage ) ) ) + "\n" )
        outf.flush()

    outf.close()
        
###################################################################
###################################################################
###################################################################
## MAIN PIPELINE
###################################################################
###################################################################
###################################################################

###################################################################
###################################################################
###################################################################
## Targets for prepare
###################################################################

###################################################################
###################################################################
###################################################################
dataset/ETHPy150Open CGATOxford/cgat/obsolete/pipeline_snps.py/recallGenomicSNPs