A quick question: I want to count the total aligned bases in a BAM file. I try to apply pysam (samtools) and pybedtools (bedtools) to do that. How can I make it?
Thanks.

**EDIT**:
I have solved this question with `pybedtools` and `pandas`:

     genomecov_base = bd.BedTool.genome_coverage(bam, dz=True) ## Depth per base on genome    
     
     gencov_perbase = pan.read_table(genomecov_base.fn, names=['chrom', 'position', 'depth'])
     gencov_allbase = (gencov_perbase.iloc[:, 2]).sum()    ## Total bases mapped on genome
Hi, if it is available you could do a quick bedtools and awk calculation:

    bedtools genomecov -ibam input.bam -bga > output_aln.bam
    awk '{s+=$4}END{print s}' output_aln.bam
    
 
this will output the total of reads at each base from your reference genome in your bam file, and then awk will just some the number of reads at each base (giving you the total aligned bases).

I hope this helps!



