def get_bedtools_coverage_cmd(bam_filename, gff_filename,
                              output_filename,
                              require_paired=False):
    """
    Get bedtools command for getting the number of reads
    from the BAM filename that are strictly contained within
    each interval of the GFF.
    """
    args = {"bam_filename": bam_filename,
            "gff_filename": gff_filename}
    # Do not include strandedness flag since that doesn't handle
    # paired-end cases 
    intersect_cmd = "bedtools intersect -abam %(bam_filename)s " \
        "-b %(gff_filename)s -f 1 -ubam " %(args)
    coverage_cmd = "%s | bedtools coverage -abam - -b %s -counts > %s" \
        %(intersect_cmd, gff_filename, output_filename)
    return coverage_cmd