rust 读取基因组文件并保存为哈希

1引言

我们尝试 使用 rust 来读取基因组文件 并 保存到哈希数据结构 里面。

2代码

这是参照 chatgpt 的代码,测试是可以看到染色体和序列分别一一对应保存到了哈希里面:

// ########################################################################
// load genome and save in hash

use std::collections::HashMap;
use std::fs::File;
use std::io::{BufRead, BufReader};
fn read_genome(genome_file: &str) -> HashMap<String, String> {
    let file = File::open(genome_file).expect("无法打开文件");
    let reader = BufReader::new(file);
    let mut seq_map: HashMap<String, String> = HashMap::new();
    let mut chr = String::new();
    let mut seq = String::new();

    for line in reader.lines() {
        let line = line.expect("读取文件失败");
        if line.starts_with('>') {  // 处理注释行
            if !chr.is_empty() {   // 保存上一个染色体的序列
                seq_map.insert(chr.clone(), seq.clone());
                seq.clear();
            }
            chr = line[1..].trim().to_string();  // 去掉注释行的 ">" 符号
        } else {
            seq.push_str(&line.trim());
        }
    }
    seq_map.insert(chr, seq);  // 保存最后一个染色体的序列
    seq_map
}

fn main() {
    let genome_file = "./test.fa";
    let seq_map = read_genome(genome_file);

    // 通过键值对的方式访问染色体序列
    for (key,val) in &seq_map {
        println!("chromosome is {}; sequnence is {}.",key,val);
    };
}
[Running] cd "c:UsersadminDesktoprust-learn" && rustc 1_readGenome.rs && "c:UsersadminDesktoprust-learn"1_readGenome
chromosome is A1BG-AS1|A1BG-AS1|NR_015380.2|1|2130|2130|NC;
sequnence is ATTTTTAGTAGAGACGGGGTTTCGTCATGTTGGGTAGACTGGTCTCGATCTCTTGACCTCATGATCCGCTGGCCT
CAGCCTCCCAAAGTGCTGGGATTATAGGCGTGAGCCACCGCACCCGGCCTCACCCTTCCATTCTTTGGGCATCCTCTGCCCAAACTCCTTTGCAACCAT.
chromosome is A1CF|A1CF|NM_001198819.2|353|2159|9481|CD;
sequnence is ATAATCAAGGAAACCTTTTCCGGGTGGGGATCTCTGAAATTACTCAGATAACAGTGCTGTGCCAAAAACCTGTGGATT
TTCTCTACAAAAATTATCATGACTGTTATCTTTTTGATGGAAAAGATGAAGTTCAGCGAAGAACAAATCTACAACTTGGCCAGAAGACTGTAGA.
chromosome is A2M|A2M|XM_006719056.4|71|4607|4953|CD;
sequnence is GGGACCAGATGGATTGTAGGGAGTAGGGTACAATACAGTCTGTTCTCCTCCAGCTCCTTCTTTCTGCAACATGGGGAAGAA
CAAACTCCTTCATCCAAGTCTGGTTCTTCTCCTCTTGGTCCTCCTGCCCACAGACGCCTCAGTCTCTGGAAAACCGCAGTATATGGTTCTGGTCCCCTC.
chromosome is A1BG|A1BG|NM_130786.4|56|1541|3382|CD;
sequnence is ATTGCTGCAGACGCTCACCCCAGACACTCACTGCACCGGAGTGAGCGCGACCATCATGTCCATGCTCGTGGTCTTTCTCTTG
CTGTGGGGTGTCACCTGGGGCCCAGTGACAGAAGCAGCCATATTTTATGAGACGCAGCCCAGCCTGTGGGCAGAGTCCGAATCACTGC.

[Done] exited with code=0 in 0.798 seconds

我们自己写一个,主要思路就是: 先建立一个空的哈希,然后遇到 > 符号开头的,就把染色体添加到哈希里作为键,如果不是 > 开头的,则把序列添加到哈希的值里。

use std::collections::HashMap;
use std::fs::File;
use std::io::{BufRead, BufReader};

// define function load genome and save in hash
fn read_genome(genome_file: &str) -> HashMap<String, String> {
    let file = File::open(genome_file).expect("can't open file!");
    let reader = BufReader::new(file);
    let mut seq_map: HashMap<String, String> = HashMap::new();
    let mut chr = String::new();

    for line in reader.lines() {
        let line = line.expect("load file failed! ");
        // 处理注释行
        if line.starts_with('>') {
            // 去掉注释行的 ">" 符号
            chr = line[1..].trim().to_string();
            seq_map.insert(chr.clone(), String::new());
        } else {
            // 获取对应键的值，并添加字符串
            if seq_map.contains_key(&chr) {
                seq_map.entry(chr.clone()).or_insert(String::new()).push_str(&line.trim());
            }
        }
    }
    seq_map
}

fn main() {
    let genome_file = "./test.fa";
    let seq_map = read_genome(genome_file);

    println!("done!");

    // 通过键值对的方式访问染色体序列
    for (key,val) in &seq_map {
        println!("chromosome is {}; sequnence is {}.",key,val);
    }
}

rust 还有关于分析生信数据的 crate,叫做 bio, 我们可以使用这个来读取基因组文件。

我们先用 cargo new bioreadgenome_1 命令建个新项目,下面是文件夹结构:

然后在 cargo.toml 文件的依赖添加:

然后在 mian.rs 里面写代码,这里我们循环打印 id 和序列:

extern crate bio;
use bio::io::fasta;
use std::fs::File;

fn main() {
    let file = File::open("C:\Users\admin\Desktop\rust-learn\test.fa").expect("can't open file!");
    let mut records = fasta::Reader::new(file).records();

    while let Some(Ok(record)) = records.next() {
        println!("id: {:?}", record.id());
        println!("seq: {:?}", String::from_utf8_lossy(record.seq()));
    }
}

cargo build 编译:

cargo run 测试运行, 运行顺利:

3结尾

有关 rust 是什么,如何安装,大家可自行网上搜索,这里不再赘述。

欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 (微信交流群需收取 20 元入群费用,一旦交费,拒不退还!(防止骗子和便于管理)) 。QQ 群可免费加入, 记得进群按格式修改备注哦。