Command for calculating average coverage for genes in a targeted sequencing experiment

calc_exon_cov.sh

# Input:
# - exons.bed: Exons BED file with the fourth column containing the gene name for gene coverage summarizing
# - sample.bam: BAM file using the same genome build (aka chromosome names and coordinate system) as exons.bed
# Output:
# - Tab-delimited file with two columns: gene name and average coverage

samtools bedcov exons.bed sample.bam \
    | awk ' \
        BEGIN { \
            FS=OFS="\t"} \
        { \
            exonbases=$5; \
            exonlen=$3-$2+1; \
            genebases[$4] += exonbases; \
            genelen[$4] += exonlen} \
        END { \
            for (gene in genebases) {
                genecov = genebases[gene]/genelen[gene]; \
                print gene,genecov}}'