grangers 0.6.0

A rust library for working with genomic ranges and annotations.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82

use flate2::{write::GzEncoder, Compression};
use grangers::grangers_utils::IntervalType;
use grangers::options::FieldColumns;
use grangers::{options, Grangers};
use noodles::fasta::record::Record;
use polars::prelude::*;
use std::io::prelude::*;
use std::pin::Pin;

#[test]
fn test_iterator() -> anyhow::Result<()> {
    let df = df!(
        "seqname" => ["chr1", "chr1", "chr1", "chr1", "chr1", "chr1", "chr1", "chr1"],
        "feature_type" => ["gene", "transcript", "exon", "exon", "gene", "transcript", "exon", "exon"],
        "start" => [1i64, 1, 1, 21, 1, 1, 1, 41],
        "end" => [30i64, 30, 10, 30, 50, 50, 10, 50],
        "strand"=> ["+", "+", "+", "+", "-", "-", "-", "-"],
        "gene_id" => ["g1", "g1", "g1", "g1", "g2", "g2", "g2", "g2"],
        "transcript_id" => [None, Some("t1"), Some("t1"), Some("t1"), None, Some("t2"), Some("t2"), Some("t2")],
        "exon_id" => [None, None, Some("e1"), Some("e2"), None, None, Some("e3"), Some("e4")],
    )
    .unwrap();

    let mut gr = Grangers::new(
        df,
        None,
        None,
        IntervalType::Inclusive(1),
        FieldColumns::default(),
        false,
    )
    .unwrap();

    let gene_id_s = gr.get_column_name("gene_id", false)?;

    let mut e = GzEncoder::new(Vec::new(), Compression::default());
    let fa_string = b">chr1\nATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCG\n";
    e.write_all(fa_string)
        .expect("can't write FASTA string to encoder");
    let compressed_fa_string = e.finish().expect("can't encode FASTA string");

    let siter = gr.iter_sequences_from_reader(
        std::io::Cursor::new(compressed_fa_string),
        false,
        Some(&gene_id_s),
        options::OOBOption::Truncate,
    )?;

    // Is this safe?
    let mut svec: Vec<Record> = Pin::into_inner(siter).map(|x| x.1).collect();

    // could also do this
    //while let Some(r) = siter.next() {
    //svec.push(r);
    //}

    println!("total records collected by iterator = {}", svec.len());

    let grangers::GrangersSequenceCollection {
        records: seq_vec,
        signature: _,
    } = gr.get_sequences_from_read(
        std::io::Cursor::new(fa_string),
        false,
        Some(&gene_id_s),
        options::OOBOption::Truncate,
    )?;

    let mut rvec: Vec<Record> = seq_vec.into_iter().map(|x| x.1).collect();

    println!("total records collected by get_sequences = {}", rvec.len());

    // the iterator version and the eager function call should return
    // equivalent sequences modulo order. Thus, after sorting, they
    // should compare as equal.
    svec.sort_unstable_by_key(|x| x.name().to_owned());
    rvec.sort_unstable_by_key(|x| x.name().to_owned());

    assert_eq!(svec, rvec);

    Ok(())
}