Gff Convertion To Wig/Bigwig

Question

[Darked89]

I mapped 40 sets of ESTs from various species to a novel plant genome using GMAP. Importing everything as separate GFF files into GBrowse will give very cluttered view, impossible to make sense by looking at it. Instead of information about 100+ individual ESTs covering particular gene it would be nice to have plot of a frequency, how often a given small interval is covered by any EST, something that is being done for RNA-Seq data (starting with SAM files).

What is a possible route to get from GFF to BigWig?

BTW, It is easy to sort GFF files based on scaffold name (first column) then "start" and "end" of a match, but iterating over each base and checking how many times it is covered by intervals seems less obvious at this stage.

Edit: title change

Answer 1

[Pierre Lindenbaum]

From http://gmod.org/wiki/Subtrack_HOWTO "An example of GFF2WIG"

#!/usr/bin/perl -w
use strict;

my $name = shift;
my $desc = shift;

@ARGV or die "I need three args: name, desc, filename";

print qq(track type=wiggle_0 name="$name" description="$desc"\n);
while (<>) {
  chomp;
  my ($ref,$start,$end,$score) = (split)[0,3,4,5];
  $ref =~ s/CHROMOSOME_|chr//;
  $start--; # zero-based, half-open
  $score = 255 if $score !~ /^[-.Ee0-9]$/;

  print join("\t",$ref,$start,$end,$score), "\n";
}

and then use wigToBigWig as described here.

[Chris Fields]

Wiggle format is not 0-based half-open: http://genome.ucsc.edu/goldenPath/help/wiggle.html

Answer 2

[Brad Chapman]

Here's an approach using BAM as an intermediate. You'll need:

• bedtools
• samtools
• wigToBigWig
• pysam
• Python script that converts BAM to BigWig coverage

along with the FASTA file of your organism of interest. Then it's a set of commandline steps using these tools:

% samtools faidx your_genome.fasta
% bedToBam -i in.gtf -g your_genome.fasta.fai >! out.bam
% samtools sort out.bam out-sort
% samtools index out-sort.bam
% bam_to_wiggle.py out-sort.bam