Bedtools Genomecoveragebed Usage : How To Create A Genome File?

Question

[komal.rathi]

I am using BEDTOOLS and the following command to get the coverage file:

$ ./genomeCoverageBed -ibam ~/GG_project/trim/ecoli.bam -g <genome file=""> > ~/GG_project/trim/coverage

where ecoli.bam is my sorted bam file, and coverage is my output file

From where do I get the genome file? How do I create a genome file?? Specifically I would need a ecoli.genome file.

Answer 1

[kirannbishwa01]

To make a genome file (for bed tools) using reference genome

1) Use samtools to generate fasta index samtools faidx lyrata_genome.fa - this will create a lyrata_genome.fa.fai (index file)

But, this index file won't work as genome file due to file format issue (mainly more than required number of columns).

2) take the index file, then use awk

awk -v OFS='\t' {'print $1,$2'} lyrata_genome.fa.fai > lyrata_genomeFile.txt

• this prints the 1st and 2nd column of a fai index file and separates the column by tab (OFS flag).
• use this file as genome file in bedtools.

if space desired between columns do this

awk {'print $1,"",$2'} lyrata_genome.fa.fai > lyrata_genomeFile.txt

if 'chr' needs to be added infront of the chromosome/scaffold names do this

awk {'print "chr"$1,"",$2'} lyrata_genome.fa.fai > lyrata_genomeFile.txt

Answer 2

[Ian]

If you are referring to the '-g' flag you can simply use the 'fetchChromSizes' script (link is to a 64bit Linux binary).

fetchChromSizes hg19 > hg19.chrom.sizes

'hg19' is an example, but any of the other UCSC genomes can be used, e.g. mm9, sacCer3, etc.

The '>' is a redirect to a new file.

The output is a tab-delimited file of chromosome name followed by it length.

I hope that helps.

Answer 3

[Medhat]

This is too old , but I just saw it :)

you can use this with the bam file itself if it has the header:

samtools view -H my.bam | grep -P '^@SQ' | cut -f 2,3 | awk 'BEGIN{OFS="\t"}{split($1, a, ":"); split($2, b, ":"); print a[2], b[2] }'

[finswimmer]

Or with sed:

$ samtools view -H file.bam|grep @SQ|sed 's/@SQ\tSN:\|LN://g'' > genome.txt

[jaeyoung.chun]

There is an unnecessary single quote at the end. Should be: $ samtools view -H file.bam|grep @SQ|sed 's/@SQ\tSN:\|LN://g' > genome.txt

[Medhat]

Actually I found also using genomeCoverageBed with Ibam do not require -g :D

[Carambakaracho]

Thanks, though newer versions of many bedtools seem generally more forgiving with regard to the -g option. Being stuck with 2.25 it's still required and the samtools faidx option requires htslib...

Answer 4

[Phoenix Mu]

I think this answer might be helpful although I am late. I saw this from the help page of genomeCoverageBed:

Tips:
One can use the UCSC Genome Browser's MySQL database to extract
chromosome sizes. For example, H. sapiens:

mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -e \
"select chrom, size from hg19.chromInfo" > hg19.genome

Answer 5

[swbarnes2]

samtools faidx output is pretty much what a genome file looks like. So run that on your reference genome.

[komal.rathi]

When I use this command: $ samtools faidx ref_sequence.fasta I get a ref_sequence.fasta.fai file which is only 32K large Its content are as follows:

gi|49175990|ref|NC_000913.2| 4639675 89 70 71

Just one line.

With this output my depth coverage plot is a bell shaped graph.

[Hatem Elshazly]

Hi, what tool did you use to generate the histogram (represented graphically as in x-axis and y-axis) from the output of genomeCoveragebed?

thanks

[swbarnes2]

Sounds just like what you should expect, right?