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
2 answers
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!
No problem! Glad to help. The -a flag means output all positions (even those with 0 coverage) which may or may not be useful downstream. You don't actually have to include it if you just want the aligned bases coverage, but sometimes it is useful if you want regions with 0 coverage too, in case you have a targeted region and you are trying to see how well it is covered.
This was a mistake on my part! it should be:
bedtools genomecov -ibam input.bam -bga > output_aln.bam
the -a flag is for bedtools coverage, you are correct! The bga flag means bed graph all so it reports positions with 0 coverage as well. If you do not want 0 coverage positions you can use the flag -bg instead. I edited the answer to reflect the change. Thanks for catching this!
This will sum the reads per base, while the question was about the number of aligned bases (i.e. each covered base should be counted once, regardless of how many reads aligns over it). In this case, just sum the length of the aligned contig (also, using just -bg should speed up the process):
bedtools genomecov -ibam input.bam -bg > output_aln.bg
awk '{ s+=($3-$2) } END { print s }' output_aln.bg
I am worried that this will significantly underestimate the number of aligned bases.
With -bga, each row has one coverage value for a range of bases, not each individual base. Illustrated in the manual here: https://bedtools.readthedocs.io/en/latest/content/tools/genomecov.html You would need to multiply the column 4 integer by the number of basepairs spanned by columns 2-3 to get total number of bp.
I think you might need to edit the first command to use -d instead of -bga to get per-base coverage:
bedtools genomecov -ibam input.bam -d > output_aln.bam
awk '{s+=$4}END{print s}' output_aln.bam
Hello vw,
beside other metrics you can get this value with CollectAlignmentSummaryMetrics from picard tools.
fin swimmer
Log in to answer this question.
-
advice for wattersons theta calculation using vcf file
written by ekirsch •
-
Run metrics of Oxford Nanopore reads and visualisation of alignment
written by giulia.trauzzi •
-
Pysam pileup and Rsamtools pileup output discrepancy
written by Xinqi •
-
Calculating coverage/depth for nanopore dRNAseq reads with reference transcriptome/genome
written by kenneditodd •
-
How to obtain zero-based coordinates read depth using bedtools coverage for a specific region?
written by vladimiralan •
-
samtools depth and sambamba depth produce different results
written by Johannes •
-
bedtools genomecov max function not working correctly
written by RBright21 •
-
Single Cell Sequencing Quality
written by hsbinf •
-
How to compute target aligned reads/total aligned reads (read enrichment) with Python?
written by vw •
-
What is the definition of "read depth" vs "coverage"? (again...)
written by ariel.balter
Hi, this sounds like something I need as well. Does this count the actual number of reads that contain a particular base or just how many reads aligned at that position? For example at position 1 if I know I have 100 reads aligned there, would it tell me I have 50 A's, 25 a, 25 C? Something like that? Or does it just say that there are 100 bases aligned there? Cheers :)
Hi DNAngel, it does not return the per base output, but one option I would recommend for such a calculation is sambamba depth. This will return the per base output at each, although the files become quite large.