samtools coverage usage
0
0
Entering edit mode
19 months ago

Hello,

I want to see the coverage of every bases and want the results in percentage, for e.g. 100% coverage covered, 99% coverage, etc. But I got the following results:

#rname  startpos    endpos  numreads covbases coverage meandepth meanbaseq  meanmapq
NC_000001.10    1   249250621   10793257    50016111    20.0666 6.34772 36.4    55.7
NT_113878.1 1   106433  2390    21426   20.131  3.27799 36.3    11.6
NT_167207.1 1   547496  21064   123882  22.627  5.60021 36.5    3.43


Here coverage is 20, 20, 22. What does it mean? Is it 20% or 20X. I couldnt understand. How can I get the results as 100% coverage?

NGS • 2.1k views
1
Entering edit mode

In case of samtools coverage it looks like it is a percentage, so 20% (http://www.htslib.org/doc/samtools-coverage.html). If you want the coverage (depth) per base I think samtools depth is the best option. An alternative can be samtools mpileup. Keep an eye on the settings of depth you may need to change something because otherwise it uses some default settings.

0
Entering edit mode

I used normal samtools coverage input.bam and samtools depth input.bam. Can you please tell me what parameters should I use? if 20% is the coverage, then it means poor coverage, right?

0
Entering edit mode

And I checked this link. If 9.7 is the coverage, that doesnt mean it is 9%. There should be some concept behind it.

0
Entering edit mode

If you want to see the finer details of your results you may need to use IGV (https://software.broadinstitute.org/software/igv/). Based on this part:

I want to see the coverage of every bases and want the results in percentage

I would say that you need samtools depth or mpileup. I may do not understand your question but based on what other number you want to calculate a percentage? In other words, if position 1 has a coverage/depth of 20 reads how do you want to calculate the percentage?

0
Entering edit mode

Thats what I am asking. How do we calculate the coverage in percentage? the coverage covered by variant should be 90-100%, so I am asking that only. How to calculate that?

0
Entering edit mode

you can't .

there is no point in putting the per-base coverage in percentage. As others have already asked: what would be for instance 100% for you?

what you can put in % is the fraction of a region that is covered by at least one read. Is it that you're after?

0
Entering edit mode

Yeah I just want to see my variants should cover 100% coverage or atleast 99% coverage

0
Entering edit mode

Please explain what you mean with "100% coverage" ?

0
Entering edit mode

I meant gene coverage. For example EGFR is the gene found in patient's data and coverage is 100%. So I want to know how to calculate gene coverage

0
Entering edit mode

Ah! You want the percentage of variants? In this case you need a vcf file. If you dont have that yet you need to do a step back. And do variant calling first. You could use freebayes for example.

EDIT: OP gave new info again, still leave it here just in case.

0
Entering edit mode

No not new info, I was being just clear. I need gene coverage. Can you please tell me how to do that?

0
Entering edit mode

IF you want gene coverage you need to dig further. I usually use is deepTools and bedtools coverage :

First you can use bamCoverageand then you can get the counts four our genes in the bam file using bedtools coverage with - counts option.

0
Entering edit mode

what is -a file and -b file here? which file I have to use?

1
Entering edit mode

these things are very clearly explained in the manual. So please do make use of it!

0
Entering edit mode

0
Entering edit mode

why did you remove your command? It would be helpful for me.

1
Entering edit mode

For two reasons: 1) I posted them by mistake 2) I want to write the tools and not my pipeline. I don't have any problem in share it but I give you a lot of information and my pipeline may not suits your data

0
Entering edit mode

This is why I asked first the options you used for it... Anyway, the best option to get the coverage for every single base in a bed file, is to use the option depth eg. samtools depth -aa -b your_regions.bed your_input.bam

0
Entering edit mode

I already mentioned the command. I didnt use any option.

0
Entering edit mode

Have you tried depth command? have a look into this post, there are so many options

0
Entering edit mode

yeah I have used

0
Entering edit mode

Which options did you use? When you ask for help with a specific results, is always better to write the code or pipeline you use. If not, is very complicated to understand the output