Solexa Eland Export format to BED file
Here is a perl script which converts from Eland export format to BED format.
#!/usr/bin/perl
# Program to convert eland export format to BED format
# for running MACS
# Chris Seidel, June 2009
#
# Requires tab delim file of chromosome or contig names
# (eland fa match files) in the format:
# UCSC_chr_name chr_length eland_name
# corrects for alignments that go off the ends of the chrs
# negative bases are trimmed to 1,
# bases > chr_length are set to chr_length
# (I know the former exist, I don't know if the latter exist)
# results are not sorted, but can be sorted in linux by:
# sort -o infile.bed -k 1,1 -k 2,2n infile.bed
# (sort in place, first column, then by second column numeric)
die("usage: $0 chrmap.txt eland_export.txt") unless(scalar(@ARGV) == 2);
# create output filename
$outfile = $ARGV[1];
$outfile =~ s/\.txt$/\.bed/;
open(FOUT, ">$outfile") || die("can't open output file: $outfile");
# get info on chromosomes
open(cmap, $ARGV[0]) || die("no chromosome name mapping file!");
%chrmap = {};
while($line = <cmap>){
chomp($line);
($newval, $size, $oldval) = split(/\t/, $line);
$chrmap{$oldval} = $newval;
$chrsize{$oldval} = $size;
}
# open input file
open(fp, $ARGV[1]) || die("can't open eland file");
$lines = 0;
while(<fp>){
chop;
++$lines;
@bits = split(/\t/);
# skip reads that didn't pass filtering
next if($bits[21] eq "N");
# get match name
$seqname = $bits[10];
# skip No Matches or QC failures
# next if($seqname =~ /NM|QC/);
# skip repeat matches
# next if($seqname =~ /\d+:\d+:\d+/);
# we're only interested in sequences that match our chrs
next unless(exists($chrmap{$seqname}));
$seqlen = length($bits[8]);
$start = $bits[12];
$end = $start + $seqlen - 1;
$strand = $bits[13];
# parse match descriptor
$n = ($bits[14] =~ tr/[ACGTN]/[ACGTN]/);
# skip reads beyond a certain threshold
next if($n > 2);
$read_code = "U".$n;
# correct for alignments off the chromosome ends
if( $start <= 0 ){
print STDERR "start less than or equal to 0: ", $start, "\n";
print STDERR join("\t", @bits), "\n";
$start = 1;
}
if($end > $chrsize{$seqname}){
print STDERR "end greater than chr end $chrsize{$seqname}: $end, diff: ", $end - $chrsize{$seqname}, "\n";
print STDERR join("\t", @bits), "\n";
$end = $chrsize{$seqname};
}
if($strand eq "F"){
$strand = "+";
$color = "0,0,255";
}
else{
$strand = "-";
$color = "255,0,0";
}
$score = 0;
print FOUT join("\t", $chrmap{$seqname}, $start, $end, $read_code, $score, $strand, $start, $end, $color), "\n";
# give some feedback
print STDERR "$lines processed\n" if(!($lines % 100000));
}
close(FOUT);
print STDERR "output file: $outfile\n";
it takes a companion file of chromosome descriptions. This file is used by the program to map the chromosome names used for the alignment - which are often big ugly names to nice short chromosome names. It also has the length of each chromosome. It's a tab delimited text file that you create for your organism to reflect the alignment values.
chr2L 23011544 Drosophila_melanogaster.BDGP5.4.54.dna.chromosome.2L.fa chr2LHet 368872 Drosophila_melanogaster.BDGP5.4.54.dna.chromosome.2LHet.fa chr2R 21146708 Drosophila_melanogaster.BDGP5.4.54.dna.chromosome.2R.fa chr2RHet 3288761 Drosophila_melanogaster.BDGP5.4.54.dna.chromosome.2RHet.fa chr3L 24543557 Drosophila_melanogaster.BDGP5.4.54.dna.chromosome.3L.fa chr3LHet 2555491 Drosophila_melanogaster.BDGP5.4.54.dna.chromosome.3LHet.fa chr3R 27905053 Drosophila_melanogaster.BDGP5.4.54.dna.chromosome.3R.fa chr3RHet 2517507 Drosophila_melanogaster.BDGP5.4.54.dna.chromosome.3RHet.fa chr4 1351857 Drosophila_melanogaster.BDGP5.4.54.dna.chromosome.4.fa chrX 22422827 Drosophila_melanogaster.BDGP5.4.54.dna.chromosome.X.fa chrXHet 204112 Drosophila_melanogaster.BDGP5.4.54.dna.chromosome.XHet.fa chrYHet 347038 Drosophila_melanogaster.BDGP5.4.54.dna.chromosome.YHet.fa