@transform( "run*.bed.gz",
               suffix(".bed.gz"),
               ".bam" )
    def buildBAM( infile, outfile ):
        '''build bam formatted files from bed files.
        The resultant files are sorted and indexed.
        The BAM files are filled with dummy sequences
        and quality scores.
        '''

        tmpfilename1 = P.getTempFilename()
        tmpfilename2 = P.getTempFilename()

        outf = open( tmpfilename1, "w" )
        fasta = IndexedFasta.IndexedFasta( PARAMS["genome"] )
        
        # convert to SAM
        for x,b in enumerate( Bed.iterator( gzip.open( infile ) )):
            l = b.end - b.start
            strand = b.fields[2]
            if strand == "+": flag = "0"
            elif strand == "-": flag = "16"
            else: raise ValueError("unknown strand in %s" % str(b))
            outf.write( "\t".join( 
                ( "read%i" % x,
                  flag,
                  b.contig,
                  str(b.start+1),
                  "10",
                  "%iM" % l,
                  "*", "0", "0",
                  "A" * l,
                  "A" * l,
                  l % "MD:Z:%i" )) + "\n" )
        
        outf.close()
            
        statement = '''
        samtools import <( cut -f 1,4 %(genome)s.idx ) %(tmpfilename1)s %(tmpfilename2)s >& %(outfile)s.log
        '''
        P.run( **dict( locals().items() + PARAMS.items() ) )

        prefix = outfile[:-len(".bam")]
        pysam.sort( tmpfilename2, prefix ) 
        pysam.index( outfile )

        os.unlink( tmpfilename1 )
        os.unlink( tmpfilename2 )

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