Newer
Older
#!/usr/bin/env perl
# Designed to work on BED data such as that available from UCSC or 1KG:
# ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase1/analysis_results/supporting/accessible_genome_masks
#
use strict;
use warnings;
use Bio::EnsEMBL::Analysis;
use Bio::EnsEMBL::SimpleFeature;
use Bio::EnsEMBL::DBSQL::DBAdaptor;
use Bio::EnsEMBL::DBSQL::SimpleFeatureAdaptor;
use Bio::EnsEMBL::Utils::IO qw/iterate_file/;
use Getopt::Long;
my ($file,$db_name,$db_host,$db_user,$db_pass,$db_port,$db_version,$help,$species,$group);
"db_name|dbname|database=s" => \$db_name,
"db_host|dbhost|host=s" => \$db_host,
"db_user|dbuser|user|username=s" => \$db_user,
"db_pass|dbpass|pass|password=s" => \$db_pass,
"db_port|dbport|port=s" => \$db_port,
"species=s" => \$species,
'group=s' => \$group,
'logic_name=s' => \$logic_name,
'description=s' => \$description,
'display_label=s' => \$display_label,
"h!" => \$help,
"help!" => \$help,
);
if ($help) {&usage; exit 0;}
unless ($file and $db_name and $db_host) {print "Insufficient arguments\n"; &usage; exit 1;}
unless ($logic_name) { print "No logic name given\n"; usage(); exit 1; }
my $dba = Bio::EnsEMBL::DBSQL::DBAdaptor->new(
-species => $species,
-group => $group,
-dbname => $db_name,
-host => $db_host,
-user => $db_user,
-pass => $db_pass,
if(! -f $file) {
die "No file found at $file";
}
process_file($file);
return;
}
sub process_file {
my $analysis = get_Analysis();
my @features;
my $count = 0;
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
my ($line) = @_;
if($count != 0 && $count % 500 == 0) {
printf STDERR "Processed %s records\n", $count;
}
chomp $line;
my $sf = line_to_SimpleFeature($line, $analysis);
push(@features, $sf);
$count++;
});
my $sfa = $dba->get_SimpleFeatureAdaptor();
print STDERR "Storing\n";
$sfa->store(@features);
print STDERR "Done\n";
return;
}
sub line_to_SimpleFeature {
my ($line, $analysis) = @_;
my ($chr, $start, $end, $label, $score, $strand) = split(/\t/, $line);
$start++; # UCSC is 0 idx start
$score ||= 0;
$strand ||= 1;
my $slice = get_Slice($chr);
my $sf = Bio::EnsEMBL::SimpleFeature->new(
-start => $start,
-end => $end,
-score => $score,
-analysis => $analysis,
-slice => $slice,
-strand => $strand,
-display_label => $label,
);
return $sf;
}
my %slices;
sub get_Slice {
my ($original) = @_;
my $name = $original;
$name =~ s/^chr//;
return $slices{$name} if exists $slices{name};
my $slice = $dba->get_SliceAdaptor()->fetch_by_region('toplevel', $name);
if(!$slice) {
die "Could not get a Slice from the Ensembl database for the given region '$original' or '$name' and coorindate system 'toplevel'. Check your core database";
}
$slices{$name} = $slice;
return $slice;
}
sub get_Analysis {
my $aa = $dba->get_AnalysisAdaptor();
my $analysis = $aa->fetch_by_logic_name($logic_name);
if(!$analysis) {
$analysis = Bio::EnsEMBL::Analysis->new(
-displayable => 1,
-logic_name => $logic_name,
);
$analysis->description($description) if $description;
$analysis->display_label($display_label) if $display_label;
$aa->store($analysis);
}
return $analysis;
}
sub usage {
print "Launching instructions:
Run from a folder you are happy to have filled with files.
Description:
Import data from a BED file into the simple_feature table. Only supports
6 column BED files (location, name and score).
perl import_bed_simple_feature.pl -file [PATH} -db_name NAME
-directory Location to download to; defaults to current directory
-host
-dbhost
-db_user Username for the DB
-user
-username
-dbuser
-db_pass Password for the DB
-pass
-password
-dbpass
-db_port Port for the DB
-dbport
-port
-species Name of the species; defaults to human
-group Name of the DB group; defaults to core
-description Analysis description; only needed if analysis is not already in the DB
-display_label Analysis display label for the website; only needed if analysis is not already in the DB
-help
";
}