Bed File Of Reads To Wiggle File Of Numbers Of Features At Each Position

Question

[KCC]

This is probably really simple, but how do you go from a BED file of genomic features to a wiggle file which contains the number of features at a given position on the chromosome? So, every position in the wiggle file gives the total number of features that overlap with that position in the chromosome.

Answer 1

[brentp]

In additions to @David's solution, if your data is large, you can get directly to bigwig format by using

bedtools genomecov (genomeCoverageBed) and

then bedGraphToBigWig

[David Langenberger]

By the way... is there also a possibility to distinguish between the strands? With my script it's possible, but absolutely not memory efficient for a whole genome, when working with NGS data... but with genomeCoverageBed it's not doable, is it?

[brentp]

David, genomeCoverageBed has a -strand option so you can specify to check reads only from a single strand. Even if it didn't, one could just use awk to filter to a specific strand, then pipe to bedtools.

[David Langenberger]

ah ok thx, didn't know that. But anyways, I always use a script, since I also want to normalize for multiple hits... I work with short RNAseq and there it happens a lot. And that is not possible with such tools.

Answer 2

[David Langenberger]

I think the simplest way is to write a small perl-script...

Copy this to a file called something like bed2wig.pl:

use strict;
use warnings;

my $file = $ARGV[0];
my %hash = ();

open(FILE, "<$file") || die "cannot open $file\n";
while(<FILE>){
    chomp;
    my ($chr, $start, $end, $id, $score, $strand) = split(/\s+/,$_);
    for(my $i=$start; $i<=$end; $i++){
        $hash{$chr}{$i}++;
    }
}
close(FILE);

foreach my $chr (keys %hash){
    print "variableStep chrom=$chr\n";
    foreach my $i (sort {$a<=>$b} keys %{$hash{$chr}}){
        print "$i\t$hash{$chr}{$i}\n";
    }
}

Now call it like ./bed2wig.pl yourFile.bed > yourFile.wig

It's not very memory efficient for big files (because of the hash), but it should give you an idea of how to do it.