if ( exists $bp_record->{ "SCORE" } ) {
$bed_entry[ score ] = $bp_record->{ "SCORE" };
+ } elsif ( exists $bp_record->{ "BIT_SCORE" } ) {
+ $bed_entry[ score ] = $bp_record->{ "BIT_SCORE" };
} else {
return wantarray ? @bed_entry : \@bed_entry;
}
}
+sub bed_upload_to_ucsc
+{
+ # Martin A. Hansen, September 2007.
+
+ # Upload a BED file to the UCSC database.
+
+ my ( $tmp_dir, # temporary directory
+ $file, # file to upload,
+ $options, # argument hashref
+ $append, # flag indicating table should be appended
+ ) = @_;
+
+ # Returns nothing.
+
+ my ( $args, $table, $sql_file, $fh_out, $fh_in );
+
+ if ( $append ) {
+ $args = join " ", $options->{ "database" }, $options->{ "table" }, "-tmpDir=$tmp_dir", "-oldTable", $file;
+ } else {
+ $args = join " ", $options->{ "database" }, $options->{ "table" }, "-tmpDir=$tmp_dir", $file;
+ }
+
+ if ( $options->{ "table" } =~ /rnaSecStr/ )
+ {
+ $table = $options->{ "table" };
+
+ print qq(uploading secondary structure table:"$table"\n) if $options->{ "verbose" };
+
+ $sql_file = "$tmp_dir/upload_RNA_SS.sql";
+
+ $fh_out = Maasha::Filesys::file_write_open( $sql_file );
+
+ print $fh_out qq(
+CREATE TABLE $table (
+ bin smallint not null, # Bin number for browser speedup
+ chrom varchar(255) not null, # Chromosome or FPC contig
+ chromStart int unsigned not null, # Start position in chromosome
+ chromEnd int unsigned not null, # End position in chromosome
+ name varchar(255) not null, # Name of item
+ score int unsigned not null, # Score from 0-1000
+ strand char(1) not null, # + or -
+ size int unsigned not null, # Size of element.
+ secStr longblob not null, # Parentheses and '.'s which define the secondary structure
+ conf longblob not null, # Confidence of secondary-structure annotation per position (0.0-1.0).
+ #Indices
+ INDEX(name(16)),
+ INDEX(chrom(8), bin),
+ INDEX(chrom(8), chromStart)
+);
+ );
+
+ close $fh_out;
+
+ Maasha::Common::run( "hgLoadBed", "-notItemRgb -sqlTable=$sql_file $options->{ 'database' } $options->{ 'table' } -tmpDir=$tmp_dir $file > /dev/null 2>&1" );
+
+ unlink $sql_file;
+ }
+ else
+ {
+ Maasha::Common::run( "hgLoadBed", "$args > /dev/null 2>&1" );
+ }
+}
+
+
+sub bed_analyze
+{
+ # Martin A. Hansen, March 2008.
+
+ # Given a bed record, analysis this to give information
+ # about intron/exon sizes.
+
+ my ( $entry, # BED entry
+ ) = @_;
+
+ # Returns hashref.
+
+ my ( $i, @begs, @lens, $exon_max, $exon_min, $exon_len, $exon_tot, $intron_max, $intron_min, $intron_len, $intron_tot );
+
+ $exon_max = 0;
+ $exon_min = 9999999999;
+ $intron_max = 0;
+ $intron_min = 9999999999;
+
+ $entry->{ "EXONS" } = $entry->{ "BLOCK_COUNT" };
+
+ @begs = split /,/, $entry->{ "Q_BEGS" };
+ @lens = split /,/, $entry->{ "BLOCK_LENS" };
+
+ for ( $i = 0; $i < $entry->{ "BLOCK_COUNT" }; $i++ )
+ {
+ $exon_len = $lens[ $i ];
+
+ $entry->{ "EXON_LEN_$i" } = $exon_len;
+
+ $exon_max = $exon_len if $exon_len > $exon_max;
+ $exon_min = $exon_len if $exon_len < $exon_min;
+
+ $exon_tot += $exon_len;
+ }
+
+ $entry->{ "EXON_LEN_-1" } = $exon_len;
+ $entry->{ "EXON_MAX_LEN" } = $exon_max;
+ $entry->{ "EXON_MIN_LEN" } = $exon_min;
+ $entry->{ "EXON_MEAN_LEN" } = int( $exon_tot / $entry->{ "EXONS" } );
+
+ $entry->{ "INTRONS" } = $entry->{ "BLOCK_COUNT" } - 1;
+ $entry->{ "INTRONS" } = 0 if $entry->{ "INTRONS" } < 0;
+
+ if ( $entry->{ "INTRONS" } )
+ {
+ for ( $i = 1; $i < $entry->{ "BLOCK_COUNT" }; $i++ )
+ {
+ $intron_len = $begs[ $i ] - ( $begs[ $i - 1 ] + $lens[ $i - 1 ] );
+
+ $entry->{ "INTRON_LEN_" . ( $i - 1 ) } = $intron_len;
+
+ $intron_max = $intron_len if $intron_len > $intron_max;
+ $intron_min = $intron_len if $intron_len < $intron_min;
+
+ $intron_tot += $intron_len;
+ }
+
+ $entry->{ "INTRON_LEN_-1" } = $intron_len;
+ $entry->{ "INTRON_MAX_LEN" } = $intron_max;
+ $entry->{ "INTRON_MIN_LEN" } = $intron_min;
+ $entry->{ "INTRON_MEAN_LEN" } = int( $intron_tot / $entry->{ "INTRONS" } );
+ }
+
+ return wantarray ? %{ $entry } : $entry;
+}
+
+
+sub bed_merge_entries
+{
+ # Martin A. Hansen, February 2008.
+
+ # Merge a list of given BED entries in one big entry.
+
+ my ( $entries, # list of BED entries to be merged
+ ) = @_;
+
+ # Returns hash.
+
+ my ( $i, @q_ids, @q_begs, @blocksizes, @new_q_begs, @new_blocksizes, %new_entry );
+
+ @{ $entries } = sort { $a->{ "CHR_BEG" } <=> $b->{ "CHR_BEG" } } @{ $entries };
+
+ for ( $i = 0; $i < @{ $entries }; $i++ )
+ {
+ Maasha::Common::error( qq(Attempted merge of BED entries from different chromosomes) ) if $entries->[ 0 ]->{ "CHR" } ne $entries->[ $i ]->{ "CHR" };
+ Maasha::Common::error( qq(Attempted merge of BED entries from different strands) ) if $entries->[ 0 ]->{ "STRAND" } ne $entries->[ $i ]->{ "STRAND" };
+
+ push @q_ids, $entries->[ $i ]->{ "Q_ID" } || sprintf( "ID%06d", $i );
+
+ if ( exists $entries->[ $i ]->{ "Q_BEGS" } )
+ {
+ @q_begs = split ",", $entries->[ $i ]->{ "Q_BEGS" };
+ @blocksizes = split ",", $entries->[ $i ]->{ "BLOCK_LENS" };
+ }
+ else
+ {
+ @q_begs = 0;
+ @blocksizes = $entries->[ $i ]->{ "CHR_END" } - $entries->[ $i ]->{ "CHR_BEG" } + 1;
+ }
+
+ map { $_ += $entries->[ $i ]->{ "CHR_BEG" } } @q_begs;
+
+ push @new_q_begs, @q_begs;
+ push @new_blocksizes, @blocksizes;
+ }
+
+ map { $_ -= $entries->[ 0 ]->{ "CHR_BEG" } } @new_q_begs;
+
+ %new_entry = (
+ CHR => $entries->[ 0 ]->{ "CHR" },
+ CHR_BEG => $entries->[ 0 ]->{ "CHR_BEG" },
+ CHR_END => $entries->[ -1 ]->{ "CHR_END" },
+ REC_TYPE => "BED",
+ BED_LEN => $entries->[ -1 ]->{ "CHR_END" } - $entries->[ 0 ]->{ "CHR_BEG" } + 1,
+ BED_COLS => 12,
+ Q_ID => join( ":", @q_ids ),
+ SCORE => 999,
+ STRAND => $entries->[ 0 ]->{ "STRAND" } || "+",
+ THICK_BEG => $entries->[ 0 ]->{ "THICK_BEG" } || $entries->[ 0 ]->{ "CHR_BEG" },
+ THICK_END => $entries->[ -1 ]->{ "THICK_END" } || $entries->[ -1 ]->{ "CHR_END" },
+ COLOR => 0,
+ BLOCK_COUNT => scalar @new_q_begs,
+ BLOCK_LENS => join( ",", @new_blocksizes ),
+ Q_BEGS => join( ",", @new_q_begs ),
+ );
+
+ return wantarray ? %new_entry : \%new_entry;
+}
+
+
+sub bed_split_entry
+{
+ # Martin A. Hansen, February 2008.
+
+ # Splits a given BED entry into a list of blocks,
+ # which are returned. A list of 6 column BED entry is returned.
+
+ my ( $entry, # BED entry hashref
+ ) = @_;
+
+ # Returns a list.
+
+ my ( @q_begs, @blocksizes, $block, @blocks, $i );
+
+ if ( exists $entry->{ "BLOCK_COUNT" } )
+ {
+ @q_begs = split ",", $entry->{ "Q_BEGS" };
+ @blocksizes = split ",", $entry->{ "BLOCK_LENS" };
+
+ for ( $i = 0; $i < @q_begs; $i++ )
+ {
+ undef $block;
+
+ $block->{ "CHR" } = $entry->{ "CHR" };
+ $block->{ "CHR_BEG" } = $entry->{ "CHR_BEG" } + $q_begs[ $i ];
+ $block->{ "CHR_END" } = $entry->{ "CHR_BEG" } + $q_begs[ $i ] + $blocksizes[ $i ] - 1;
+ $block->{ "Q_ID" } = $entry->{ "Q_ID" } . sprintf( "_%03d", $i );
+ $block->{ "SCORE" } = $entry->{ "SCORE" };
+ $block->{ "STRAND" } = $entry->{ "STRAND" };
+ $block->{ "BED_LEN" } = $block->{ "CHR_END" } - $block->{ "CHR_BEG" } + 1,
+ $block->{ "BED_COLS" } = 6;
+ $block->{ "REC_TYPE" } = "BED";
+
+ push @blocks, $block;
+ }
+ }
+ else
+ {
+ @blocks = @{ $entry };
+ }
+
+ return wantarray ? @blocks : \@blocks;
+}
+
+
+sub bed_overlap
+{
+ # Martin A. Hansen, February 2008.
+
+ # Checks if two BED entries overlap and
+ # return 1 if so - else 0;
+
+ my ( $entry1, # hashref
+ $entry2, # hashref
+ $no_strand, # don't check strand flag - OPTIONAL
+ ) = @_;
+
+ # Return bolean.
+
+ return 0 if $entry1->{ "CHR" } ne $entry2->{ "CHR" };
+ return 0 if $entry1->{ "STRAND" } ne $entry2->{ "STRAND" };
+
+ if ( $entry1->{ "CHR_END" } < $entry2->{ "CHR_BEG" } or $entry1->{ "CHR_BEG" } > $entry2->{ "CHR_END" } ) {
+ return 0;
+ } else {
+ return 1;
+ }
+}
+
+
+
+
# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
projects / biopieces.git / blobdiff