This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Gene coverage from Bam files

Hi,

I would like to extract gene coverage from a bacterial genome that i have assembled. I have done the following steps:

I have created a bed file with all genes i´m interested in from my gff file (corresponding to my assembled genome. Only 1 full chromosome). My bed file is genome.map.bed. Sequence1 is the name of the bacterial chromosome. Last column represents gene_names.

sequence1       1625    2738    Prodigal_2
sequence1       2749    2956    Prodigal_3
sequence1       2961    4041    Prodigal_4
sequence1       4086    5988    Prodigal_5
sequence1       5996    8459    Prodigal_6
sequence1       8451    8754    Prodigal_7
sequence1       8757    9090    Prodigal_8
sequence1       9092    9881    Prodigal_9
sequence1       10482   10872   Prodigal_10

I have mapped paired-end reads back to my genome assembly using bwt and got a bam file. The output file is 1.mapped.bam.sorted

Next, i´m trying to extract coverage information from bam file as follows:

bedtools coverage -hist -abam 1.mapped.bam.sorted -b genome.map.bed > Sample.hist

Below the header lines of my Sample.hist file.

sequence1 0 100
J00173:25:HJ3VYBBXX:4:1101:12418:6765/2 60 + 0 100
0,0,0 1 100, 0, 0 100 100 1.0000000 sequence1 0 41
J00173:25:HJ3VYBBXX:4:1105:11617:35092/1 60 + 0
41 0,0,0 1 41, 0, 0 41 41
1.0000000 sequence1 0 112 J00173:25:HJ3VYBBXX:4:1105:24870:19654/1 60 + 0
112 0,0,0 1 112, 0, 0 112 112
1.0000000 sequence1 0 129 J00173:25:HJ3VYBBXX:4:1106:32309:31875/2 60 + 0
129 0,0,0 1 129, 0, 0 129 129
1.0000000 sequence1 0 62 J00173:25:HJ3VYBBXX:4:1107:7273:15891/1 60 + 0 62
0,0,0 1 62, 0, 0 62 62 1.0000000

I was expecting to get coverage for each gene by using both the bed file and bam file but each row represents a read instrad of a gene. What´s wrong ??

Thanks,

bedtools
• 2,789 views • link
updated by finswimmer finswimmer • written by David David

1 answer

Hello,

you have to use the bed file for the -a parameter and the bam file for the -b parameter to get the expected result.

fin swimmer

• 1 views • link
written by finswimmer finswimmer

Thanks, Indeed the following command did work:

bedtools coverage -hist -b 1.mapped.bam.sorted -a genome.map.bed > Sample.hist

• 1 views • link
written by David David

Log in to answer this question.

More posts like this