Question: Bedtools Genomecoveragebed Usage : How To Create A Genome File?
4
gravatar for komal.rathi
5.6 years ago by
komal.rathi3.3k
Children's Hospital of Philadelphia, Philadelphia, PA
komal.rathi3.3k wrote:

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.

bedtools • 19k views
ADD COMMENTlink modified 10 weeks ago by Medhat8.0k • written 5.6 years ago by komal.rathi3.3k
8
gravatar for kirannbishwa01
2.3 years ago by
United States
kirannbishwa01890 wrote:

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

ADD COMMENTlink modified 2.2 years ago • written 2.3 years ago by kirannbishwa01890
3
gravatar for Ian
5.6 years ago by
Ian5.3k
University of Manchester, UK
Ian5.3k wrote:

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.

ADD COMMENTlink written 5.6 years ago by Ian5.3k
1
gravatar for swbarnes2
5.6 years ago by
swbarnes24.5k
United States
swbarnes24.5k wrote:

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

ADD COMMENTlink written 5.6 years ago by swbarnes24.5k

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.

ADD REPLYlink modified 5.6 years ago • written 5.6 years ago by komal.rathi3.3k
1

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

ADD REPLYlink written 4.5 years ago by Hatem Elshazly60

Sounds just like what you should expect, right?

ADD REPLYlink written 5.6 years ago by swbarnes24.5k
0
gravatar for Medhat
10 weeks ago by
Medhat8.0k
Poland
Medhat8.0k wrote:

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] }'

ADD COMMENTlink written 10 weeks ago by Medhat8.0k

Or with sed:

$ samtools view -H file.bam|grep @SQ|sed 's/@SQ\tSN:\|LN://g'' > genome.txt
ADD REPLYlink written 10 weeks ago by finswimmer7.9k

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

ADD REPLYlink written 9 weeks ago by Medhat8.0k
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: 1490 users visited in the last hour