@split( PARAMS["filename_vcf"], "*.pileup.gz" )
    def buildPileups( infile, outfile ):
        '''build samtools pileup formatted files from vcf formatted files.

        The column to strain mapping are determined dynamically.

        #CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  129P2   129S1   129S5   AKR     A_J     BALB    C3H     C57BL   CAST    CBA     DBA     LP_J    NOD     NZO     PWK     SPRET   WSB
        #CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  129P2   129S5   AKR     C3H     CAST    DBA     NOD     PWK     WSB     129S1   A_J     BALB    C57BL   CBA     LP_J    NZO     SPRET

        The parser is incomplete (only looks at unphased data, etc.)

        IT ALSO IGNORES HETEROZYGOUS CALLS.

        Both vcf and pileup employ 1-based coordinate systems. Adds "chr" prefix.

        This method applies two potential filters:

        1. filename_refseq_filter: 
               remove all variants not in refseq regions
        2. filename_snp_filter: 
               remove all SNPs in a blacklist. Note that the code
               assumes that the blacklist is not too large.

        '''

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

        if "filename_snp_filter" in PARAMS:
            def _defdict(): return collections.defaultdict( list )
            filter_snps = collections.defaultdict( _defdict )
            x = 0
            for line in IOTools.openFile(PARAMS["filename_snp_filter"], "r" ):
                if line.startswith("track"):continue
                if line.startswith("#"):continue
                
                data = line[:-1].split("\t")
                track, contig, pos = data[0], data[1], int(data[2])
                track = track[len("mouse"):]
                filter_snps[track][contig].append( pos )
                x += 1

            E.info("removing %i false positive SNPs" % x )
        else:
            filter_snps = None

        if "filename_refseq_filter" in PARAMS:
            E.info( "reading segment filter")
            intervals = GTF.readAndIndex( GTF.iterator( IOTools.openFile(PARAMS["filename_refseq_filter"], "r") ) )
            E.info( "read segment filter")
        else:
            intervals = None

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

        for line in inf:
            data = line[:-1].split("\t")
            if line.startswith("#CHROM"):
                if not headers: headers = data[9:]
                continue
            elif line.startswith("#"):
                continue

            contig, pos, ref = data[0], data[1], data[3]

            pos = int(pos)
            variants = [ref]
            variants.extend( data[4].split(",") )
            counts.input += 1

            contig = "chr%s" % contig

            if intervals:
                if not intervals.contains( contig, pos-1, pos ):
                    counts.filter += 1
                    continue
                
            for h, genotype_info in zip(headers, data[9:]):

                # no variant for this strain - skip
                if genotype_info == "." or genotype_info.startswith("./."): continue

                # determine the genotype base - this is a hard-coded order
                # revise if input file formats change.
                consensus_quality, genotype_quality, read_depth = "0", "0", "0"

                dd = genotype_info.split(":")

                if len(dd) == 5:
                    genotype, mapping_quality, hcg, genotype_quality, read_depth = dd
                    if hcg != "1": continue
                elif len(dd) == 4:
                    genotype, mapping_quality, genotype_quality, read_depth = dd
                elif len(dd) == 2:
                    genotype, genotype_quality = dd
                elif len(dd) == 1:
                    genotype = dd[0]
                else:
                    raise ValueError( "parsing error for %s: line=%s" % (genotype_info, line) )

                genotype = genotype.split("/")

                if filter_snps:
                    if pos-1 in filter_snps[h][contig]:
                        counts.filtered_snps += 1
                        continue
                
                # ignore heterozygous calls
                if len(set(genotype)) != 1: continue

                genotype = [ variants[int(x)] for x in genotype ]
                lengths = [len(x) for x in genotype] + [len(ref)]
                is_snp = len( set(lengths) ) == 1 and lengths[0] == 1

                # skip genotypes for which no call can be made
                if "." in genotype: continue

                if is_snp:
                    genotype = "".join(genotype)

                    # skip wild type 
                    if genotype == "%s%s" % (ref,ref):
                        continue
                    
                    outfiles.write( h, 
                                "\t".join( map(str, (
                                    contig,
                                    pos,
                                    ref,
                                    Genomics.encodeGenotype( genotype ),
                                    consensus_quality,
                                    genotype_quality,
                                    mapping_quality,
                                    read_depth,
                                    genotype,
                                    "<" * len(genotype) ) ) ) + "\n" )
                else:

                    def getPrefix( s1, s2 ):
                        '''get common prefix of strings s1 and s2.'''
                        n = min( len( s1), len( s2 ) )
                        predix = []
                        for x in range( n ):
                            if s1[x] != s2[x]: return s1[:x]
                        return s1[:n]

                    def getSuffix( s1, s2 ):
                        '''get common sufix of strings s1 and s2.'''
                        n = min( len( s1), len( s2 ) )
                        predix = []
                        if s1[-1] != s2[-1]: return ""
                        for x in range( -2, -n - 1, -1 ):
                            if s1[x] != s2[x]: return s1[x+1:]
                        return s1[-n:]

                    def getGenotype( variant, ref ):
                        if variant == ref: return "*", 0

                        if len(ref) > len(variant):
                            # is a deletion
                            if ref.startswith(variant):
                                return "-%s" % ref[len(variant):], len(variant) - 1
                            elif ref.endswith( variant ):
                                return "-%s" % ref[:-len(variant)], -1
                            else:
                                prefix = getPrefix( ref, variant )
                                suffix = getSuffix( ref, variant )
                                shared = len(prefix) + len(suffix) - len(variant) 
                                # print "-", prefix, suffix, ref, variant, shared, len(prefix), len(suffix), len(ref)
                                if shared < 0:
                                    raise ValueError()
                                return "-%s" % ref[len(prefix):-(len(suffix)-shared)], len(prefix) - 1

                        elif len(ref) < len(variant):
                            # is an insertion
                            if variant.startswith(ref):
                                return "+%s" % variant[len(ref):], len(ref) - 1
                            elif variant.endswith(ref):
                                return "+%s" % variant[:len(ref)], 0
                            else:
                                prefix = getPrefix( ref, variant )
                                suffix = getSuffix( ref, variant )
                                shared = len(prefix) + len(suffix) - len(ref) 
                                if shared < 0:
                                    raise ValueError()

                                return "+%s" % variant[len(prefix):-(len(suffix)-shared)], len(prefix)
                        else:
                            assert 0, "snp?"

                    # in pileup, the position refers to the base
                    # after the coordinate, hence subtract 1
                            #pos -= 1

                    genotypes, offsets = [], []
                    is_error = True
                    for variant in genotype:
                        try:
                            g, offset = getGenotype( variant, ref ) 
                        except __HOLE__:
                            break

                        assert len(g) > 1, "incomplete genotype %s" % g
                        genotypes.append( g )
                        offsets.append( offset )
                    else: 
                        is_error = False
                    if is_error: 
                        print line,
                        counts.errors += 1
                        continue

                    assert len(set(offsets )) == 1
                    offset = offsets[0]

                    genotypes = "/".join( genotypes )
                    outfiles.write( h, 
                                "\t".join( map(str, (
                                contig,
                                pos + offset,
                                "*",
                                genotypes,
                                "0",
                                "0",
                                "0",
                                "0",
                                genotypes,
                                "<" * len(genotype),
                                "0", 
                                "0",
                                "0") ) ) + "\n" )

            counts.output += 1
        outfiles.close()

        E.info("%s" % str(counts))
dataset/ETHPy150Open CGATOxford/cgat/obsolete/pipeline_snps.py/buildPileups