Entering edit mode
5.0 years ago
Mehmet
▴
820
Dear All,
I have mapped reads of several samples to genome for SNP calling and now I would like to calculate coverage of each gene in each sample (reads) output (bam files).
I have checked bedtools, samtools and GATK with -T DepthOfCoverage option by providing a set of gene names, but none of them provided what I want.
Any tool you can suggest?
Thank you.
but none of them provided what I want
So what exactly do you want? Please provide expected output examples.OK. I want coverage of each gene in each bam file (each sample ) that were produced by mapping of DNA reads to the reference genome of the same species.
OK. Sorry I forgot. Below an example of read count of each gene in one sample generated by htseq. For coverage I want to have something like that. Coverage per gene in a sample.
the first question - what is a gene for you? How is it defined? Is it canonical gene coordinates, or all the exons/transcripts, or just region between gene's start/end?
starting and ending region.
so you extract starting + ending coordinates of your genes in BED format (chrom, start, end) and then run bedtools/samtools/ngs-bits. may be sounds a bit old-school, but it will work.
You can have the average depth or the amount of reads.
It is possible that not every basepair in that gene has the same depth. You can get the depth per basepair with
samtools depth
after that you need to make your own script that calculates the average depth per region. O a tool already exist but dont know that.Also see this post: samtools bedcov vs. bedtools coverage