Merge branch 'main' of https://github.com/TonyBolger/HelixerPost

Author: TonyBolger

LICENSE

@@ -1,27 +1,51 @@
 # HelixerPost
 ## Dependencies
+### Rust
+A recent version of Rust (see https://www.rust-lang.org/tools/install)
+
 ### hd5 libraries
 On fedora & co. you will need `hdf5-devel`, while on ubuntu & co. you will need `libhdf5-dev`.
-Currently it's also necessary to manually build the 'lzf' compression support from h5py (https://pypi.org/project/h5py/) as a shared library and install to the hdf5 plugins directory
+
+### hd5 lzf support (Skip if unsure)
+If you need the 'lzf' compression support (no longer needed on recent Helixer versions), you will need to download it from h5py (https://pypi.org/project/h5py/) and manually build it as a shared library and install to the hdf5 plugins directory. This will also require a C toolchain - the system provided GCC should be fine. 
 
 `tar -xzvf h5py-3.2.1.tar.gz`
 
 `cd h5py-3.2.1/lzf/`
 
-### Fedora
+**Fedora**
+
 `gcc -O2 -fPIC -shared -Ilzf lzf/*.c lzf_filter.c -lhdf5 -o liblzf_filter.so`
 
 `sudo mkdir -p /usr/local/hdf5/lib/plugin`
 
 `sudo cp liblzf_filter.so /usr/local/hdf5/lib/plugin`
 
-### Ubuntu
+**Ubuntu**
+
 `gcc -O2 -fPIC -shared -Ilzf -I/usr/include/hdf5/serial/ lzf/*.c lzf_filter.c -lhdf5 -L/lib/x86_64-linux-gnu/hdf5/serial -o liblzf_filter.so`
 
 `sudo mkdir /usr/lib/x86_64-linux-gnu/hdf5/plugins`
 
 `sudo cp liblzf_filter.so /usr/lib/x86_64-linux-gnu/hdf5/plugins`
 
+## Building HelixerPost
+
+`git clone https://github.com/TonyBolger/HelixerPost.git`
+
+`cd HelixerPost`
+
+`cargo build --release`
+
+The resulting binary is './targets/release/helixer_post_bin'
+
+Running `./target/release/helixer_post_bin` should show the command line parameters, giving:
+
+`HelixerPost <genome.h5> <predictions.h5> <window_size> <edge_thresh> <peak_thresh> <min_coding_length> <output.gff>`
+
+In order for Helixer to find this binary, it needs to be on the PATH. The easiest way to achieve this is to copy 
+the binary to the bin folder in the virtual environment which you previously created for Helixer 
+(e.g. `path_to_Helixer/env/bin` )
 
 ## Concept
 HelixerPost uses a sliding window assessment to determine regions of the genome which are likely gene containing.
@@ -35,8 +59,6 @@ genic score within the window exceeds the edge threshold, and continues until th
 that window. The candidate region is accepted if it also contains at least one window with a genic score above 
 the required peak threshold.
 
-
-
 ## Parameters
 
 `HelixerPost <genome.h5> <predictions.h5> <window_size> <edge_thresh> <peak_thresh> <min_coding_length> <output.gff>`
@@ -58,20 +80,13 @@ this value (e.g. 60)
 
 output.gff: The path for the output GFF file 
 
-
 ## Example Usage
 
 Using 100bp sliding window for genic detection, 0.1 edge threshold, 0.8 peak threshold (before HMM)
 
 60bp minimum CDS length per gene (after HMM)
 
-### One Step:
-
-`cargo run --release genome_data.h5 predictions.h5 100 0.1 0.8 60 output.gff`
-
-### Two Step: 
-
-`cargo build --release`
+From the HelixerPost dir:
 
 `./target/release/helixer_post_bin genome_data.h5 predictions.h5 100 0.1 0.8 60 output.gff`
 

README.md

@@ -0,0 +1,38 @@
+# Example HelixerPost data
+
+## source
+The following was created from the first 21384 base pairs
+of NC\_023163.2 in the Oryza\_brachyantha genome, version
+GCF\_000231095.2\_ObraRS2 from refseq.
+
+The genome was preprocessed into numeric form (genome\_data.h5) with `fasta2h5py`,
+and base-wise genic predictions (predictions.h5) were created with `HybridModel.h5`
+both from the Helixer project (https://github.com/weberlab-hhu/Helixer).
+
+## usage
+From this directory
+
+```
+helixer_post_bin genome_data.h5 predictions.h5 100 0.1 0.8 60 output.gff3
+```
+
+## expected output
+
+Using the provided data, the above should create the file 'output.gff3',
+and the contents there of are expected to be as follows (as of commit b674af9):
+
+```
+##gff-version 3.2.1
+##species Oryza_brachyantha
+# f4f3f7bf629d02431241bfef18cb8ade  /home/ali/.local/share/Helixer/models/land_plant/land_plant_v0.3_a_0100.h5
+##sequence-region NC_023163.2:1-21384 1 21384
+NC_023163.2:1-21384	Helixer	gene	17026	18898	.	+	.	ID=Oryza_brachyantha_NC_023163.2:1-21384_000001
+NC_023163.2:1-21384	Helixer	mRNA	17026	18898	.	+	.	ID=Oryza_brachyantha_NC_023163.2:1-21384_000001.1;Parent=Oryza_brachyantha_NC_023163.2:1-21384_000001
+NC_023163.2:1-21384	Helixer	exon	17026	17356	.	+	.	ID=Oryza_brachyantha_NC_023163.2:1-21384_000001.1.exon.1;Parent=Oryza_brachyantha_NC_023163.2:1-21384_000001.1
+NC_023163.2:1-21384	Helixer	five_prime_UTR	17026	17115	.	+	.	ID=Oryza_brachyantha_NC_023163.2:1-21384_000001.1.five_prime_UTR.1;Parent=Oryza_brachyantha_NC_023163.2:1-21384_000001.1
+NC_023163.2:1-21384	Helixer	CDS	17116	17356	.	+	0	ID=Oryza_brachyantha_NC_023163.2:1-21384_000001.1.CDS.1;Parent=Oryza_brachyantha_NC_023163.2:1-21384_000001.1
+NC_023163.2:1-21384	Helixer	exon	17707	18898	.	+	.	ID=Oryza_brachyantha_NC_023163.2:1-21384_000001.1.exon.2;Parent=Oryza_brachyantha_NC_023163.2:1-21384_000001.1
+NC_023163.2:1-21384	Helixer	CDS	17707	18305	.	+	2	ID=Oryza_brachyantha_NC_023163.2:1-21384_000001.1.CDS.2;Parent=Oryza_brachyantha_NC_023163.2:1-21384_000001.1
+NC_023163.2:1-21384	Helixer	three_prime_UTR	18306	18898	.	+	.	ID=Oryza_brachyantha_NC_023163.2:1-21384_000001.1.three_prime_UTR.1;Parent=Oryza_brachyantha_NC_023163.2:1-21384_000001.1
+```
+

example/example.md

@@ -5,5 +5,5 @@ authors = ["Tony Bolger <bolger.tony@gmail.com>"]
 edition = "2018"
 
 [dependencies]
-hdf5 = "0.7.1"
-ndarray = "0.14.0"
+hdf5 = "0.8.1"
+ndarray = "0.15.6"

example/genome_data.h5

@@ -204,12 +204,12 @@ impl<W: Write> GffWriter<W> {
     pub fn write_global_header(
         &mut self,
         species: Option<&str>,
-        helixer_model_md5sum: Option<&str>,
+        helixer_model_md5sum: Option<String>,
     ) -> std::io::Result<()> {
         const GFF_VERSION: &'static str = "3.2.1";
         write!(self.writer, "##gff-version {}\n", GFF_VERSION)?;
         if let Some(species) = species {
-            write!(self.writer, "##{}\n", species)?;
+            write!(self.writer, "##species {}\n", species)?;
         }
         if let Some(helixer_model_md5sum) = helixer_model_md5sum {
             write!(self.writer, "# {}\n", helixer_model_md5sum)?;

example/predictions.h5

@@ -3,9 +3,11 @@ use helixer_post_bin::analysis::hmm::show_hmm_config;
 use helixer_post_bin::analysis::rater::SequenceRating;
 use helixer_post_bin::analysis::Analyzer;
 use helixer_post_bin::gff::GffWriter;
+use helixer_post_bin::results::raw::RawHelixerPredictions;
 use helixer_post_bin::results::HelixerResults;
 use std::fs::File;
 use std::io::BufWriter;
+use std::path::Path;
 use std::process::exit;
 
 fn main() {
@@ -56,8 +58,11 @@ fn main() {
         1,
         "Error: Multiple Species are not allowed for GFF output."
     );
-    let model_md5sum = None; // TODO: this should fetch <hdf5>.attrs['model_md5sum']
+    let rhg = RawHelixerPredictions::new(&Path::new(predictions_path))
+        .expect("Error: Something went wrong while accessing the predictions.");
+    let model_md5sum = rhg.get_model_md5sum().ok();
     let species_name = helixer_res.get_all_species().first().map(|x| x.get_name());
+
     gff_writer
         .write_global_header(species_name, model_md5sum)
         .expect(&*format!(

helixer_post_bin/Cargo.toml

@@ -50,7 +50,7 @@ impl<'a, T: H5Type + Clone + Copy> BlockedDataset1D<'a, T> {
 
     fn get_data_for_block(&self, block_id: BlockID) -> hdf5::Result<Array1<T>> {
         let slice = s![block_id.inner(), ..];
-        self.dataset.read_slice_1d::<T, _>(&slice)
+        self.dataset.read_slice_1d(slice)
     }
 }
 
@@ -164,7 +164,7 @@ impl<'a, T: ArrayConvInto<O>, O> BlockedDataset2D<'a, T, O> {
 
     fn get_data_for_block(&self, block_id: BlockID) -> hdf5::Result<Array2<T>> {
         let slice = s![block_id.inner(), .., ..];
-        self.dataset.read_slice_2d::<T, _>(&slice)
+        self.dataset.read_slice_2d(slice)
     }
 }
 

helixer_post_bin/src/gff.rs

@@ -28,8 +28,8 @@ pub enum ErrSamples
 }
 */
 
-type SeqidsType = FixedAscii<[u8; 50]>;
-type SpeciesType = FixedAscii<[u8; 25]>;
+type SeqidsType = FixedAscii<50>;
+type SpeciesType = FixedAscii<25>;
 
 impl RawHelixerGenome {
     pub fn new(