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.
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
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?
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.