@split( PARAMS["filename_snps"], "*.pileup.gz" )
    def importSNPs( infile, outfile ):

        outfiles = IOTools.FilePool( "mouse%s.pileup" )

        inf = gzip.open(infile,"r")
        headers = []
        counts = E.Counter()

        if "filename_refseq_filter" in PARAMS:
            intervals = GTF.readAndIndex( GTF.iterator( IOTools.openFile(PARAMS["filename_refseq_filter"], "r") ) )
        else:
            intervals = None
        
        for line in inf:
            data = line[:-1].split("\t")
            if line.startswith("#"):
                if not headers: headers = data[3:]
                continue

            counts.input += 1

            contig, pos, ref = data[:3]
            pos = int(pos)

            if intervals:
                if not intervals.contains( "chr%s" % contig, pos-1, pos ):
                    counts.filter += 1
                    continue

            for h, genotype in zip(headers, data[3:]):
                if genotype == "..": continue
                outfiles.write( h, 
                                "\t".join( map(str, (
                                contig,
                                pos,
                                ref,
                                Genomics.encodeGenotype( genotype ),
                                "0",
                                "0",
                                "0",
                                "0",
                                genotype,
                                "<" * len(genotype) ) ) ) + "\n" )
                
                counts.output += 1

        outfiles.close()

        E.info("%s" % str(counts) )

        E.info("starting compression and indexing" )

        # convert to bgzip and index with tabix
        outfiles = glob.glob( "mouse*.pileup" )
        
        for outfile in outfiles:
            E.info("compressing and indexing %s" % outfile )
            pysam.tabix_index( outfile, preset = "vcf" )
build samtools pileup formatted files from tabular output
dataset/ETHPy150Open CGATOxford/cgat/obsolete/pipeline_snps.py/importSNPs