code: 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()
docs: %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