本crate提供表示和操作基因组数据序列的类型。提供了针对核苷酸和氨基酸的有效编码，并可以通过Codec trait进行扩展。

对固定长度的短序列（kmer）给予了特殊关注。

快速入门

将bio-seq添加到Cargo.toml

[dependencies]
bio-seq = "0.13"

迭代kmer以获取序列

use bio_seq::prelude::*;

// The `dna!` macro packs a static sequence with 2-bits per symbol at compile time
let seq = dna!("ATACGATCGATCGATCGATCCGT");

// iterate over the 8-mers of the reverse complement
for kmer in seq.revcomp().kmers::<8>() {
    println!("{kmer}");
}

// ACGGATCG
// CGGATCGA
// GGATCGAT
// GATCGATC
// ATCGATCG
// ...

序列类似于rust的字符串类型，遵循类似的重引用约定

use bio_seq::prelude::*;

// The `dna!` macro packs a static sequence with 2-bits per symbol at compile time.
// This is analogous to rust's string literals:
let s: &'static str = "hello!";
let seq: &'static SeqSlice<Dna> = dna!("CGCTAGCTACGATCGCAT");

// Static sequences can also be copied as kmers
let kmer: Kmer<Dna, 18> = dna!("CGCTAGCTACGATCGCAT").into();
// or with the kmer! macro:
let kmer = kmer!("CGCTAGCTACGATCGCAT");

// `Seq`s can be allocated on the heap like `String`s are:
let s: String = "hello!".into();
let seq: Seq<Dna> = dna!("CGCTAGCTACGATCGCAT").into();

// Alternatively, a `Seq` can be fallibly encoded at runtime:
let seq: Seq<Dna> = "CGCTAGCTACGATCGCAT".try_into().unwrap();

// &SeqSlices are analogous to &str, String slices:
let slice: &str = &s[1..3];
let seqslice: &SeqSlice<Dna> = &seq[2..4];

哲学

许多生物信息学crate实现自己的kmer打包逻辑。这个项目最初是为了重用kmer构建代码并使其在不同项目间兼容。很快人们意识到，没有紧密耦合到一般数据类型的kmer类型是没有意义的。crate的范围将被限制在表示序列。

有些人喜欢设计聪明的位操作技巧来反转序列互补，有些人想要快速原型化紧凑的数据结构。大多数人不想担心字节序。Rust的优势在于我们可以安全地将科学抽象化，以合作的方式同时实现这两个目标。

欢迎贡献。有很多低垂的果实可以优化，理想情况下我们只需要写一次！

内容

• Codec：生物序列符号的编码/解码方案
• Seq：编码符号的序列
• Kmer：长度为K的固定大小序列
• 可推导codec：此crate提供用于定义自己的位级编码的实用工具
• 序列之间的安全转换

序列

编码符号的字符串被打包到Seq中。切片、分块和窗口返回SeqSlice。类型Seq<<A: Codec>和&SeqSlice<<A: Codec>类似于String和&str。与标准字符串类型一样，它们存储在堆上。Kmer通常存储在栈上，实现Copy。

Kmers

kmers是长度为k的短序列，可以适应寄存器（例如usize或SIMD向量）并且通常实现Copy。这些通过const generics实现，并且k在编译时是固定的。

所有数据都以小端格式存储。这影响了序列映射到整数（“字典序”顺序）的顺序。

for i in 0..=15 {
    println!("{}: {}", i, Kmer::<Dna, 5>::from(i));
}

0: AAAAA
1: CAAAA
2: GAAAA
3: TAAAA
4: ACAAA
5: CCAAA
6: GCAAA
7: TCAAA
8: AGAAA
9: CGAAA
10: GGAAA
11: TGAAA
12: ATAAA
13: CTAAA
14: GTAAA
15: TTAAA

紧凑编码

可以通过将kmers直接作为usize来在常数时间内索引查找表。

// This example builds a histogram of kmer occurences

fn kmer_histogram<C: Codec, const K: usize>(seq: &SeqSlice<C>) -> Vec<usize> {
    // For dna::Dna our histogram will need 4^4
    // bins to count every possible 4-mer
    let mut histo = vec![0; 1 << (C::WIDTH * K as u8)];

    for kmer in seq.kmers::<K>() {
        histo[usize::from(kmer)] += 1;
    }

    histo
}

Sketching

最小化

2位表示法的核苷酸顺序为A < C < G < T。序列和kmers以小端格式存储，并且按“字典序”排序。这意味着AAAA < CAAA < GAAA < ... < AAAC < ... < TTTT。

let seq = dna!("GCTCGATCGTAAAAAATCGTATT");
let minimiser = seq.kmers::<8>().min().unwrap();

assert_eq!(minimiser, Kmer::from(dna!("GTAAAAAA")));

哈希

Hash为序列和kmer类型实现，使得这些类型的相等的值具有相同的哈希值

let seq_arr: &SeqArray<Dna, 32, 1> = dna!("AGCGCTAGTCGTACTGCCGCATCGCTAGCGCT");
let seq: Seq<Dna> = seq_arr.into();
let seq_slice: &SeqSlice<Dna> = &seq;
let kmer: Kmer<Dna, 32> = seq_arr.into();

assert_eq!(hash(seq_arr), hash(&seq));
assert_eq!(hash(&seq), hash(&seq_slice));
assert_eq!(hash(&seq_slice), hash(&kmer));

哈希最小化器

在实践中，我们希望对最小化的序列进行哈希

fn hash<T: Hash>(seq: T) -> u64 {
    let mut hasher = DefaultHasher::new();
    seq.hash(&mut hasher);
    hasher.finish()
}

let (minimiser, min_hash) = seq
    .kmers::<16>()
    .map(|kmer| (kmer, hash(&kmer)))
    .min_by_key(|&(_, hash)| hash)
    .unwrap();

规范kmers

在最小化时考虑kmers的前向和反向互补

let (canonical_minimiser, canonical_hash) = seq
    .kmers::<16>()
    .map(|kmer| {
        let canonical_hash = min(hash(&kmer), hash(&kmer.revcomp()));
        (kmer, canonical_hash)
    })
    .min_by_key(|&(_, hash)| hash)
    .unwrap();

尽管分别最小化seq和seq.revcomp()可能更有效率。

编码器

Codec特质描述了生物序列符号的编码/解码过程。这个特质可以通过程序性推导。有四个内置的编码器

codec::Dna

使用字典序排列的2位表示法

codec::Iupac

IUPAC核苷酸模糊代码用4位表示。这支持通过位操作进行成员资格解析。逻辑or是并集

assert_eq!(iupac!("AS-GYTNA") | iupac!("ANTGCAT-"), iupac!("ANTGYWNA"));

逻辑and是两个iupac序列的交集

assert_eq!(iupac!("ACGTSWKM") & iupac!("WKMSTNNA"), iupac!("A----WKA"));

codec::Text

可以从常见的纯文本文件格式中直接读取的utf-8字符串可以被视为序列。可以定义额外的逻辑以确保'a'=='A'以及处理'N'。

codec::Amino

氨基酸序列用6位表示。氨基酸的表示设计成容易从2位编码的DNA序列中转换过来。

定义新的编解码器

可以通过实现Codec特质来定义自定义编解码器。

编解码器也可以从枚举类型的变体名称和区分符派生而来

use bio_seq_derive::Codec;
use bio_seq::codec::Codec;

#[derive(Clone, Copy, Debug, PartialEq, Codec)]
#[repr(u8)]
pub enum Dna {
    A = 0b00,
    C = 0b01,
    G = 0b10,
    T = 0b11,
}

注意，您需要在枚举中显式提供“区分符”（例如，0b00）。

#[width(n)]属性指定每个符号所需的位数。支持的最大值是8。如果没有指定此属性，则将选择最优宽度。

#[alt(...,)]和#[display('x')]属性可以用来定义替代表示或使用特殊字符显示项目。以下是codec::Amino中终止密码子的定义。

pub enum Amino {
    #[display('*')] // print the stop codon as a '*'
    #[alt(0b001011, 0b100011)] // TGA, TAG
    X = 0b000011, // TAA (stop)

序列转换

翻译表特质

翻译表提供了将密码子转换为氨基酸的方法。

在Cargo.toml中启用翻译功能

[dependencies]
bio-seq = { version="0.13", features=["translation"] }

pub trait TranslationTable<A: Codec, B: Codec> {
    fn to_amino(&self, codon: &SeqSlice<A>) -> B;
    fn to_codon(&self, amino: B) -> Result<Seq<A>, TranslationError>;
}

/// A partial translation table where not all triples of characters map to amino acids
pub trait PartialTranslationTable<A: Codec, B: Codec> {
    fn try_to_amino(&self, codon: &SeqSlice<A>) -> Result<B, TranslationError>;
    fn try_to_codon(&self, amino: B) -> Result<Seq<A>, TranslationError>;
}

标准遗传密码以translation::STANDARD常量提供

use crate::prelude::*;
use crate::translation::STANDARD;
use crate::translation::TranslationTable;

let seq = dna!("AATTTGTGGGTTCGTCTGCGGCTCCGCCCTTAGTACTATGAGGACGATCAGCACCATAAGAACAAA");

let aminos: Seq<Amino> = seq
    .windows(3)
    .map(|codon| STANDARD.to_amino(&codon))
    .collect::<Seq<Amino>>();

assert_eq!(
    aminos,
    Seq<Amino>::try_from("NIFLCVWGGVFSRVSLCARGALSPRAPPLL*SVYTLYM*ERGDTRDISQSAHTPHI*KRENTQK").unwrap()
);

自定义翻译表

从实现Into<HashMap<Seq<A>, B>>的类型的实例化翻译表

let codon_mapping: [(Seq<Dna>, Amino); 6] = [
    (dna!("AAA"), Amino::A),
    (dna!("ATG"), Amino::A),
    (dna!("CCC"), Amino::C),
    (dna!("GGG"), Amino::E),
    (dna!("TTT"), Amino::D),
    (dna!("TTA"), Amino::F),
];

let table = CodonTable::from_map(codon_mapping);

let seq: Seq<Dna> = dna!("AAACCCGGGTTTTTATTAATG");
let mut amino_seq: Seq<Amino> = Seq::new();

for codon in seq.chunks(3) {
    amino_seq.push(table.try_to_amino(codon).unwrap());
}

直接实现TranslationTable特质

struct Mitochondria;

impl TranslationTable<Dna, Amino> for Mitochondria {
    fn to_amino(&self, codon: &SeqSlice<Dna>) -> Amino {
        if *codon == dna!("AGA") {
            Amino::X
        } else if *codon == dna!("AGG") {
            Amino::X
        } else if *codon == dna!("ATA") {
            Amino::M
        } else if *codon == dna!("TGA") {
            Amino::W
        } else {
            Amino::unsafe_from_bits(Into::<u8>::into(codon))
        }
    }

    fn to_codon(&self, _amino: Amino) -> Result<Seq<Dna>, TranslationError> {
        unimplemented!()
    }
}

贡献

欢迎贡献和建议。开始贡献的最简单方法是通过github问题。