Gene 12 Column BED file

[obifro]

Hello,
I'd like a gene level 12 column bed file for Ensembl v65 Mus musculus, however I cannot find such a file (those in the ensembl archive site are 12 column transcript bed files).

This file should contain one 12 column bed entry per gene, with all exons from all transcripts merged if necessary (eg boundaries extended if the exons overlap).

Does anyone know if such a file exists?

Alternatively, I can make one, but would like to avoid having to code it myself. At my disposal I have a GTF file, and can make a 12 column transcript BED file. BEDtools has the mergeBed tool, however this merges by exon coordinates, not gene ID.

Has anyone coded this?

Cheers!

[obifro]

I ended up coding it up. Only briefly tested but seems to do the trick, hope it's useful if anyone else needs such a script. Requires Bedtools.

#! /usr/bin/perl -w

# merge the exons from all transcripts for a given gene into a single 12 col bed line
# for each individual gene.

# Requires BEDTools to be installed and in $PATH

# ARGV[0] is a transcript BED file (eg as downloaded from ensembl)
# ARGV[1] is a tab delimited text file with geneID<tab>transcriptID for every transcript
# ARGV[2] is a working directory for tmp files.

# prints new BED file to STDOUT.

use strict;
use warnings;
use Data::Dumper;

my $usage = "usage: perl $0 <transcriptBedFile> <transcript2GeneFile> <workingdir>\n";

my $txBedFile = shift or die $usage;
my $tx2GeneFile = shift or die $usage;
my $workingdir = shift or die $usage;


# load all the transcripts and gene info
my %gene2tx = ();
open(TX2GENE, $tx2GeneFile) or die;
while (<TX2GENE>) {
	chomp;
	my @d = split/\t/;
	push(@{$gene2tx{$d[0]}}, $d[1]);
}
close TX2GENE;

# load all the bed lines:
my %bed = ();
open(TXBED, $txBedFile) or die;
while (<TXBED>) {
	my @d = split/\t/;
	$bed{$d[3]} = $_;
}
close TXBED;


# for each gene write all 12 col transcript lines, transform to a 6 col bed, sort, then merge exons, load, and write new
# bed file.
foreach my $gene (keys %gene2tx) {
	my @transcripts = @{$gene2tx{$gene}};
	my $geneTx12Bed = "$workingdir/$gene.tx.12.bed";
	my $geneTx6Bed = "$workingdir/$gene.tx.6.bed";
	my $geneExonsBed = "$workingdir/$gene.exons.merged.bed";
	
	open(TX12, ">$geneTx12Bed") or die "Failed to write to $geneTx12Bed: $!\n";
	foreach my $tx (@transcripts) { print TX12 $bed{$tx}; }
	close TX12;
	
	system "bed12ToBed6 -i $geneTx12Bed | sort -k1,1 -k2,2n > $geneTx6Bed";
	system "mergeBed -i $geneTx6Bed | sort -k1,1 -k2,2n > $geneExonsBed";

	my $strand = `cut -f 6 $geneTx6Bed | sort | uniq`;
	chomp $strand;

	open(EXONS, "$geneExonsBed") or die "Faield to read from $geneExonsBed: $!\n";
	my $geneChr = undef;
	my $geneChrStart = undef;
	my $geneChrEnd = undef;
	my @exonStarts = ();
	my @exonLengths = ();
	while (<EXONS>) {
		chomp;
		my ($chr, $start, $end) = split /\t/;
		unless(defined($geneChrStart)) {
			$geneChrStart = $start;
			$geneChr = $chr;
		}

		push(@exonStarts, $start - $geneChrStart);
		push(@exonLengths, ($end - $start));
		$geneChrEnd = $end;
	}
	
	my @bedLine = (
		$geneChr,
		$geneChrStart,
		$geneChrEnd,
		$gene,
		0,
		$strand,
		$geneChrStart,
		$geneChrEnd,
		0,
		scalar (@exonStarts),
		join(",",@exonLengths).",",
		join(",",@exonStarts).",",
		);
	
	print join("\t", @bedLine),"\n";
	
	system "rm $geneTx12Bed $geneTx6Bed $geneExonsBed";
}

[swaraj]

Hey Obifro,

Cool job!!! I was about to code for a script when I came across your post. It saved me some time. Thanks.