#fasta #fastq #生物信息学 #基因组学 #数据流

bio-streams

流式生物信息学数据类型

2 个版本

0.3.1 2024年5月14日
0.3.0 2024年5月12日

#107 in 生物学

Download history 134/week @ 2024-05-07 221/week @ 2024-05-14 25/week @ 2024-05-21 1/week @ 2024-06-11

每月下载量62次

MIT 协议

25KB
458

Docs.rs CI status

bio-steams

流式基因组数据类型和数据结构

该crate处于早期开发阶段。欢迎贡献。

Webassembly示例: 从流式fastq中移除非M. TB读取基于扩增子的SARS-CoV-2组装

特性

通过 FastqFasta 流共享 Record 类型

pub struct Record<T: for<'a> TryFrom<&'a [u8]> = Vec<u8>> {
    pub fields: Vec<u8>,
    pub seq: T,
    pub quality: Option<Vec<Phred>>, // fasta records set quality to `None`
}

可以将记录读取到自定义类型:pub struct Fastq<R: BufRead, T = Seq<Dna>>

示例

流式传输一对fastq并检查其名称字段上的某些条件

// Open a pair of gzipped fastq files as streams of `Record`s with `Seq<Dna>` sequences

let fq1: Fastq<BufReader<MultiGzDecoder<File>>> = Fastq::new(BufReader::new(
    MultiGzDecoder::new(File::open(&file1).unwrap()),
));

let fq2: Fastq<BufReader<MultiGzDecoder<File>>> = Fastq::new(BufReader::new(
    MultiGzDecoder::new(File::open(&file2).unwrap()),
));

for zipped in fq1.zip(fq2) {
    match zipped {
        (Ok(r1), Ok(r2)) => {
            // check that the last characters of the name strings are 1 and 2
            if r1.fields[r1.fields.len() - 1] != b'1' || r2.fields[r2.fields.len() - 1] != b'2'
            {
                eprintln!("paired records do not end in 1/2");
            }

            // check that the description fields are equal up to the last character
            if r1.fields[..r1.fields.len() - 1] != r2.fields[..r2.fields.len() - 1] {
                eprintln!("reads do not have the same names");
            }
        }
        _ => {
            eprintln!("Parse error in fastq files");
        }
    }
}

使用读取文件 r1.fq.gzf2.fq.gz 运行 fqcheck 示例程序

$ cargo build --example fqcheck --release
$ target/release/examples/fqcheck r1.fq.gz r2.fq.gz

计算氨基酸k-mer

// this opens a gzipped data stream and parses it into `Records` with `Seq<Amino>` sequence fields
let faa: Fasta<BufReader<File>, Seq<Amino>> =
    Fasta::new(BufReader::new(File::open(&faa_file).unwrap()));

// we can convert amino acid k-mers directly into usizes and use them to index into a table
let mut histogram = Box::new([0u64; 1 << (K * Amino::BITS as usize)]);

for contig in faa {
    // here "contig" is a fasta record
    for kmer in contig.unwrap().seq.kmers::<K>() {
        histogram[usize::from(kmer)] += 1;
    }
}

使用fasta文件 proteins.faa 运行 aminokmers 示例程序

$ cargo build --example fqcheck --release
$ target/release/examples/aminokmers proteins.faa

路线图

输入流

  • fastq
  • fasta
  • 待办事项 sam/bam
  • 待办事项 gfa

待办事项

  • 质量分数特质,Phredu8 的别称
  • 为异步使用futures::streams
  • GAT 借用迭代器
  • 基准测试
  • 示例

依赖关系

~2–2.7MB
~59K SLoC