This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to get coverage for an alignment file

hi good morning biostars I got in trouble again accepting a work that I have never performed. seems It's not that difficult so I will propose here my approximation and then you can tell me if what I'm doing is right or we can have better approximation right. I need to get coverage of whole genome from an alignment file. I really don't know how to do that possibly using bedtools. please let me know if there's an alternative to have coverage for the whole genome. thank you very much have a nice day

bam
• 5,630 views • link
updated by rlbc • rlbc • written by biotech biotech

2 answers

You can use bedtools coverage -abam <bam file> -b <bed file/gff file> -d > coverage.txt

Then you will have a coverage for every nucleotide. Because the file will be big for whole genome I suggest to write a python script for calculating of average coverage for whole genome, or coverage per gene etc. I also suggest using bed file from UCSC for whole genes.

Hope it helps,

Best,

Agata

• 0 views • link
written by agata88 agata88

Hi Agata good morning you are so helpful. In fact I just need to get coverage for each nucleotide I don't need the coverage for each gene for now. From the command you gave me seems I need a GFF file. Do you think we can skip that annotation file?

• 0 views • link
written by biotech biotech

It is not an annotation file but defined regions. I think that a coverage file for full genome would be something like 2 Gb of data?

• 0 views • link
written by Kairos Kairos

Hi thanks for your quick reply. I'm working with small bacteria is about 5 MB so things are going to be tiny. If I remember well the gff file has been produced using bedtools right. Anyway I think I should read in more detail the manual of bedtools.

• 0 views • link
written by biotech biotech

I think you can write your own bed file including coordinates start and stop of your reference. What bacteria do you have?

• 0 views • link
written by agata88 agata88

We assembled the genome in house so it's not published. Do you know what I need to generate this GFF files for bedtools? How do people generate this file?

• 0 views • link
written by biotech biotech

I think you can calculate how much nucleotides you have and then just write a bed file included:

1(start) /t XXX (end)

I am not sure it is going to work I have never done something like that but it is worth to try :) I would try that first.

• 0 views • link
written by agata88 agata88

Where are you taking this gff files from? You are so lucky you have them, maybe because you work with a model organism.

• 0 views • link
written by biotech biotech

If that file does not work cut columns 1,4,5 and that would convert it to a bed format.

This file will not give you coverage for every nucleotide for the entire genome (just for the intervals defined in your file).

If you need the entire genome I wonder if you can get away with saying 

Scaffold_1 1 End_nucleotide_#

which will give you coverage for every nucleotide in the scaffold.

• 0 views • link
written by GenoMax GenoMax

Hi genomax2. Thanks for having a look at my file. I have 31 contigs in my genome. So doing that will do the trick for all the contigs right thanks so much

• 0 views • link
written by biotech biotech

Hi Agata, do you think we are talking about gff3 or basic gff. If I remember well there are some differences.

• 0 views • link
written by biotech biotech

There is also qualimap, which gives you whole genome and per chromosome/contig coverage.

• 0 views • link
written by rlbc • rlbc

Log in to answer this question.

More posts like this