Question: Gene coverage from Bam files
0
gravatar for David
4 weeks ago by
David140
David140 wrote:

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 • 101 views
ADD COMMENTlink modified 4 weeks ago by finswimmer9.9k • written 4 weeks ago by David140
1
gravatar for finswimmer
4 weeks ago by
finswimmer9.9k
Germany
finswimmer9.9k wrote:

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

ADD COMMENTlink written 4 weeks ago by finswimmer9.9k

Thanks, Indeed the following command did work:

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

ADD REPLYlink written 4 weeks ago by David140
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1143 users visited in the last hour