def contamination(args):
    """
    %prog contamination Ecoli.fasta genome.fasta read.fastq

    Check read contamination on a folder of paired reads. Use bowtie2 to compare
    the reads against:
    1. Ecoli.fsata - this will tell us the lower bound of contamination
    2. genome.fasta - this will tell us the upper bound of contamination
    """
    from jcvi.apps.bowtie import BowtieLogFile, align

    p = OptionParser(contamination.__doc__)
    p.set_firstN()
    opts, args = p.parse_args(args)

    if len(args) != 3:
        sys.exit(not p.print_help())

    ecoli, genome, fq = args
    firstN_opt = "--firstN={0}".format(opts.firstN)
    samfile, logfile = align([ecoli, fq, firstN_opt])
    bl = BowtieLogFile(logfile)
    lowerbound = bl.rate
    samfile, logfile = align([genome, fq, firstN_opt])
    bl = BowtieLogFile(logfile)
    upperbound = 100 - bl.rate

    median = (lowerbound + upperbound) / 2

    clogfile = fq + ".Ecoli"
    fw = open(clogfile, "w")
    lowerbound = "{0:.1f}".format(lowerbound)
    upperbound = "{0:.1f}".format(upperbound)
    median = "{0:.1f}".format(median)

    print("\t".join((fq, lowerbound, median, upperbound)), file=fw)
    print("{0}: Ecoli contamination rate {1}-{2}".\
                        format(fq, lowerbound, upperbound), file=sys.stderr)
    fw.close()
def contamination ( args ) : from jcvi . apps . bowtie import BowtieLogFile , align p = OptionParser ( contamination . __doc__ ) p . set_firstN ( ) opts , args = p . parse_args ( args ) if len ( args ) != 3 : sys . exit ( not p . print_help ( ) ) ecoli , genome , fq = args firstN_opt = "--firstN={0}" . format ( opts . firstN ) samfile , logfile = align ( [ ecoli , fq , firstN_opt ] ) bl = BowtieLogFile ( logfile ) lowerbound = bl . rate samfile , logfile = align ( [ genome , fq , firstN_opt ] ) bl = BowtieLogFile ( logfile ) upperbound = 100 - bl . rate median = ( lowerbound + upperbound ) / 2 clogfile = fq + ".Ecoli" fw = open ( clogfile , "w" ) lowerbound = "{0:.1f}" . format ( lowerbound ) upperbound = "{0:.1f}" . format ( upperbound ) median = "{0:.1f}" . format ( median ) print ( "\t" . join ( ( fq , lowerbound , median , upperbound ) ) , file = fw ) print ( "{0}: Ecoli contamination rate {1}-{2}" . format ( fq , lowerbound , upperbound ) , file = sys . stderr ) fw . close ( )
%prog contamination Ecoli.fasta genome.fasta read.fastq

    Check read contamination on a folder of paired reads. Use bowtie2 to compare
    the reads against:
    1. Ecoli.fsata - this will tell us the lower bound of contamination
    2. genome.fasta - this will tell us the upper bound of contamination