snapATAC使用后的主观评测

2019年,一篇发表在Genome Biology上题为「Assessment of computational methods for the analysis of single-cell ATAC-seq data」的文章对目前一些常见的单细胞ATAC-seq分析工具进行了评测,其中SnapATAC是当时文章认为比较好的工具

Despite variation across methods and datasets, SnapATAC, Cusanovich2018, and cisTopic outperform other methods in separating cell populations of different coverages and noise levels in both synthetic and real datasets. Notably, SnapATAC is the only method able to analyze a large dataset (> 80,000 cells).

我曾经测试过这个工具,发现了其中的一个bug,后面就没有继续在用了。在2021年的最后几天,我又尝试使用了这个工具。如果你也在分析single cell ATAC-seq数据,想尝试snapATAC的话,希望这篇文章对你有所帮助。

  1. snapATAC有一个配套的工具, snapTools, 它能够将fastq文件处理成所需的snap文件,或者将CellRanger-ATAC输出BAM/TSV文件转成snap文件。对于CellRanger-ATAC的数据转换,我不推荐使用BAM文件,因为他需要先给bam文件里read name添加barcode信息,然后按照read name进行排序,消耗的时间非常惊人。
  2. snapATAC虽然能够处理大于80k的细胞,但是这和你选择的bin-size有关。snapATAC在分析的时候还是会将bin-by-cell的矩阵加载到内存中,因此内存还是一个比较重要的因素。我分析的时候,当bin-size比较小时,就会因为内存占用过多被kill。
  3. snapATAC的可视化基于plot3D, 而不是ggplot2, 如果想要在原先图基础上进行修改,不如ggplot2那么方便。
  4. snapATAC的聚类分析默认用的是R-igraph, leiden算法需要安装Python模块。R-igraph无法通过修改resolution来改变类群的数目,而leiden可以。不过,我好奇,为什么snapATAC不考虑使用Seurat的聚类算法?

如何绘制物理图谱和遗传图谱的对应关系

唐海宝老师开发的JCVI有一个工具,叫做ALLMAPS, 能够展示遗传图谱和物理图谱的对应关系,如下所示

alignment

但是这个图的目标是为了对ALLMAPS的scaffold结果进行可视化,并不是专门用于展示遗传图谱的标记和物理图谱的对应关系。尽管在allmaps这个组件下提供了plot函数,命令行输入只要求 input.bed 和 seqid, 但实际运行的时候还要求 allmaps path的中间文件, xxxx.lifted.bed, xxxx.agp, weight.txt等文件。

为了解决这一问题,我阅读了allmaps.py的源代码,在plot的基础上增加了一个plot2函数,只需要用户输入 input.bed, 染色体编号和染色体的长度就能够画图。

1
2
# python -m jcvi.assembly.allmaps plot2 input.bed 染色体编号 染色体长度
python -m jcvi.assembly.allmaps plot2 input.bed chr1 123140023

其中input.bed的格式要求有6列,

  • 标记所在染色体名
  • 标记所在染色体的start
  • 标记所在染色体的end, 通常就是start+1
  • 标记对应的图谱位置, 要求输入为”图谱名-连锁图谱所在组:连锁图谱的遗传距离”
  • 标记名

案例

1
2
3
4
Chr1    68185909        68185910        male-14:48.470000        Chr1:68185910
Chr1 68479621 68479622 male-14:49.380000 Chr1:68479622
Chr1 68595299 68595300 male-14:48.440000 Chr1:68595300

可以先生成如下的csv文件,然后转换成bed

input.csv

1
python -m jcvi.assembly.allmaps merge  male.csv male.bed

目前有plot2函数的代码还在我的项目下,xuzhougeng/jcvi, 待代码稳定了,再PR。

「荐书」入门深度学习的必备数学知识

以我浅薄的理解,深度学习是机器学习的一个子集。想要入门深度学习,可以先从最基础的神经网络开始。

神经网络也称之为多层感知机(multi-layer perceptron, MLP),由输入层,隐藏层,输出层构成。最简单的神经网络只需要3层,拥有数百到数千个参数,而现代复杂的神经网络则可能有成千万上百万层,拥有数万,数十万,数百万的参数,例如人工智能公司OpenAI在2020年7月发布的GPT-3拥有1750亿个参数。对于那么多层的神经网络的训练,我们就称之为深度学习(deep learning).

目前关于深度学习的框架很多,包括但不限于MXNET, PyTorch, TensorFlow等,因此利用这些框架训练一个手写数字识别,猫和狗分类的神经网络并不算特别困难。但是如果我们有更高的追求,就得了解下它背后的数学基础。

目前阅读到大部分课程,虽然都会提到深度学习所需要的数学知识,例如线性代数,微积分,概率论和信息论,但是要么受限于篇幅没有详细介绍,要么是因为“知识诅咒”,会直接使用高阶的抽象的数学符号介绍原理,使得我迷失各种符号定义中。

之前推荐的Python神经网络编程虽然也涉及到数学,但比较浅显,在此基础上直接去阅读动手学习深度学习这类书,还是会有很大压力,于是我找到了另一本缓冲书籍,深度学习的数学

深度学习的数学

我个人觉得这本书最大的亮点在于,它创建了一个Excel就能训练的神经网络数据集。这个数据集一共64张图片,记录数字0和1。神经网络一共三层,输入层是12个节点,对应输入图像的12个像素点,隐藏层只有一层,3个节点,输出层是2个节点对应0和1,总计 12 x 3 + 3 + 3 x 2 + 2 = 47个参数。

由于参数少,所以作者就可以展开所有的公式,即 $z = w_1x_1 + w_2x2+ b$,而不必使用线性代数中的矩阵和向量乘法, 即 $z = w \cdot x + b$

虽然前者看起来很繁琐,不如后者简洁。但,就我而言,我更喜欢前者的直观,因为我可以自己在脑中将前者转换成后者,建立两者的联系,那么,以后看到矩阵乘法的形式,就能自动展开。

同样由于参数少,我们既可以使用Excel以线性模型的方式对代价函数优化,也可以使用Excel采取(随机)梯度下降的方式对代价函数优化。相对于使用Python根据原理手动搭建简单的神经网络,使用Excel让我们对数据有了更多的掌控感。

这本书的另一个亮点是,它涉及到数学知识刚好是我能够得上的水平。例如,在我学习「动手学习深度学习」的4.7节时,我搞不懂为什么损失函数对权重参数求导式子是下面这个样子

$$
\frac{\partial J}{\partial \mathbf{W}^{(2)}}= \text{prod}\left(\frac{\partial J}{\partial \mathbf{o}}, \frac{\partial \mathbf{o}}{\partial \mathbf{W}^{(2)}}\right) + \text{prod}\left(\frac{\partial J}{\partial s}, \frac{\partial s}{\partial \mathbf{W}^{(2)}}\right)= \frac{\partial J}{\partial \mathbf{o}} \mathbf{h}^\top + \lambda \mathbf{W}^{(2)}.
$$

看了「深度学习的数学」的2.8节才知道这是多变量函数的链式法则。

学习就相当于探险,一本合适的书会告诉我们需要了解哪些知识,这就相当于在探险过程中拿到了地图,虽然我们还没有真正抵达目的地,但起码我们知道该怎么走了,心里就有底了。

「荐书」Python神经网络编程,入门神经网络的儿童读本

之前计划学习「深度学习」相关内容,于是我在JD上买了很多和这个主题有关的书,其中就包括这篇文章推荐的「Python神经网络编程」。

为什么我在标题中,称其为”入门神经网络的儿童读本”呢? 因为同时期我还买了一本书,深度学习(DEEP LEARNING),由于封面是一个花园,所以它一般也被称为”花书”。花书大约500页,每页都是密密麻麻的字,和一堆数学公式,如果你计划直接阅读花书来入门深度学习的话,就相当于让一个刚识字的孩子直接去看四大名著,可能字都认识,连起来不明白它到底想表达什么意思。

而这本书200多页,一共22多万字,分为3个主要章节和2个附录,由于文字浅显易懂,读起来不会有太多障碍,因此一个周末的时间就能读完。看这本书,相当于让孩子在直接读四大名著前,先看注音白话版(含插画),先增加自信心,提起后续继续的阅读的兴趣。

举个例子,在「深度学习」和「Python神经网络编程」都有一个关于XOR(异或)的案例,用于说明线性分类器(线性模型: y=ax+b)的局限性。在读「深度学习」的时候,我只是知道XOR不行,但是不知道为什么作者要举XOR这个例子。而「Python神经网络编程」的作者解决了我这个疑惑。

下图中,我们可以发现AND和OR能够使用一条直线,根据输入x1和x2判断输出是0还是1.

AND和OR

但是对于XOR,你无法借助于一条直线将图中的红色和绿色点进行区分,所以为什么不去试试神奇的”神经网络”呢?

XOR

另外一个例子,就是我理解了为什么线性代数在神经网络中那么重要。在「Python神经网络编程」的1.7节中,作者先带着我们一步一步在神经网络中追踪信号的传递, 先计算了第二层的第一个节点,也就是x=(第一个节点的输出 * 链接权重 ) + (第二个节点的输出 * 链接权重),然后再按照相同的思路计算了其他节点。在1.7节的最后,作者说到

…好在,即使是面对一个具有很多层,众多节点的神经网络,数学可以帮助我们以非常简洁的方式写下计算出所有输出值的命令。由于这种简洁性,指令变得非常短,执行起来也更有效率,因此这种简洁性不仅仅是对人类读者有益处,对计算机而言,也一样大有裨益。这一简洁方式就是使用矩阵…

在1.8节中,作者引入了矩阵的计算方式,然后将之前的繁琐的计算过程化简成 $X = W \cdot I$, 使得我知其然并知其所以然。

最后,我花了3天晚上把这本书看完了,让我从”花书”的打击中走了出来。这一幕就跟当年不知好歹,想通过「算法导论」 入门算法一样,最后还是靠着「图解算法」活了过去。

如何向GitHub上传超过100MB的大文件

Git是一个版本控制工具,Github是一个项目托管网站。一般来说,我们会对代码进行版本控制,而记录代码的文件通常也很小,所以,我一般也不会用这种上传大文件的需求。并且,往GitHub上传超过50MB的文件时,它会警告你,认为这会影响性能,超过100Mb,直接报错。综上,我觉得往GitHub上传大文件是我不需要掌握的技能。

直到最近,我发现我上传的一个项目里面居然有一个400Mb的数据时,我就不得不去处理之前我觉得没有必要的蠢问题。(总不可能把数据放在百度网盘让人下载吧,这不太体面)

往Github上传大文件,我们需要使用GitHub的LFS(Large File Storage)服务,其中免费版允许2GB(超过2GB就需要付流量费)。

首先,点击Github仓库的Setting中,

Setting

在其中的Archives中勾选Git LFS对应选项

Git LFS

第二步,在服务器上安装Git LFS

1
2
3
4
# centos
sudo yum install git-lfs
# ubuntu
sudo apt install git-lfs

第三步: 设置LFS

1
2
3
git lfs install
# 提示信息如下
Git LFS initialized.

第四步: 在你项目下设置你需要跟踪的大文件,例如我需要跟踪R语言保存的二进制文件,一般是以Rds, Rds, Rdata, Rda结尾的文件

1
2
3
4
5
6
git lfs track "*.rds"
git lfs track "*.Rds"
git lfs track "*.Rdata"
git lfs track "*.Rda"
# 让git跟踪记录lfs信息的文件
git add .gitattributes

第五步:正常使用Git管理你的项目,上传到GitHub

1
2
3
4
5
6
7
8
9
git add data/*.rds
git commit -m "add rds"
git push origin master
``

上传的时候,提示信息如下

```nbash
Uploading LFS objects: 0% (0/1), 3.4 MB | 129 KB/s

整体上,并不复杂,相当于在基础款的Git上安装一个插件而已。

参考资料:

对几个月前看的动漫【无职转生】的碎碎念

今年年初,B站上线了一个番剧【无职转生】,我看完了前2集后,就发现它因为一些原因下线了,考虑到动画周期太长,所以我就直接下载了已经完结的轻小说看了。

小说故事很简单,讲的就是一个因为校园霸凌事件而躲在房间里肥宅,在双亲去世那天也不出场,因此被亲戚赶出门,在路上冲去救一对可能会被车撞的情侣而被送去异世界,最后决定在异世界努力重新度过一生的故事。

更加具体一点,就是这个小说以一个传记形式讲述了主角鲁迪乌斯的一生,从幼年的无忧无欲,到青少年时在魔大陆的努力求生,青年时在学校学习,到后来为了守护自己的家庭而开始谋划各种事件。由于一生过于漫长,所以作者在鲁迪乌斯经历人生最重要的boss战后,就直接写了他最后的时光,而整个故事最触动我的一点就是他最后的时光了。

很多故事都在主角打败最终boss后结束,或者赋予主角永恒的寿命,但这本小说作者没有让鲁迪乌斯发现永生的秘密,让他打破人类寿命的上限。在他最后时光里,他做了一个梦,梦到了他的红发妻子,但醒来却发现他的妻子早已先他离世,无比的悲伤,两位寿命远长于他的妻子发现他醒来后,立马叫来了他的后代和好友,在众人的围绕下长逝。最后的最后,他还是醒了一次,以前世的躯壳见到了’人神’被他后代击败的场景。

我看到最后的时光非常的悲伤,因为他描写了我们人类终将面对的问题,’死亡’。在小说中,作者如果愿意,主角可以从创世活到世界的结束,但是在这个世界中,除非在你心中坚定人死后可以去另一个世界,否

Rust第四课:通过一个案例认识自己的不足

为了避免「一学就会,一写就废」,因此我给自己设置了一个练习题,就是读取一个只有两列的CSV文件,用第一列当做key,第一列当做value,保存到HashMap(哈希表)中。

用Python完成这个任务非常的简单,代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import csv
import sys

csv_file = sys.argv[1]

hash_map = {}
with open(csv_file, "r") as f:
for line in csv.reader(f):
key = line[0]
value = line[1]
if key in hash_map:
hash_map[key].append(value)
else:
hash_map[key] = [ value ]

根据Python的代码,我查找了rust中csv和HashMap的相关文档并学习,但是我对Rust本身的语言特性并不是特别熟悉,导致花了一天才处理完报错并理解背后的逻辑。

Rust的变量默认不可变

和我之前学的C/C++, Python, R等编程语言不同,Rust的变量默认是不可改变的,也就是下面的代码是错误

1
2
let hm = HashMap::new();
hm.insert("key","value");

你必须显式用关键词 mut声明,Rust才允许你对一个变量进行改变。

1
2
3
let mut hm = HashMap::new();
hm.insert("key","value");

一开始觉得很不习惯,但是后来想想也挺合理的。如果一个变量,你没有改变它的想法,但是你编写的某些代码却能改动它,就可能引起bug。同时,如果一个值,你认为它需要改变,但是实际上它没有变动,那就意味着你应该去掉这个mut声明。

Rust要求你时刻考虑错误处理。

Rust通过类型系统来处理异常情况,因此许多函数并不会直接返回结果,而是返回Option或Result。这就导致我们从函数里得到的值无法直接使用,需要多写一点代码从中提取目标结果。

下面代码中csv_reader的类型是 Result<Reader<File>, Error>, 因此无法使用for循环进行遍历。

1
2
3
4
5
use csv::Reader;
let mut csv_reader = Reader::from_path(csv_file);
for result in csv_reader.records() {
....
}

处理方法有两种,一种是使用?操作符将错误往后传播处理

1
2
3
for result in csv_reader.records()? {
....
}

另一种则是调用.unwrap()直接退出程序,

1
2
3
for result in csv_reader.records().unwrap() {
....
}

同样的,遍历得到result也不能直接使用,需要通过模式匹配,从中提取出目标结果

1
2
3
4
if let Ok(record) = result{
println!("{}", record);
...
}

Rust要求你注意值和借用的生命周期

在Rust中,每个值都有其生命周期,对值的引用(借用)不能超过(outlive)值的生存期。

在下面代码中,我的目标是遍历每一行,然后将x进行split构建成一个vector,然后提取k,v,进行插入

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
fn main() -> io::Result<()>{
let args: Vec<String> = env::args().collect();
let csv_file = &args[1];

let mut hm = HashMap::new();

for line in io::BufReader::new(fs::File::open(csv_file)?).lines(){
if let Ok(x) = line {
let s: Vec<&str> = x.split(',').collect();
//直接用 hm.insert(s[0], s[1])也是会报错的
let k = s[0];
let v = s[1];
hm.insert(k, v);
}
}
Ok(())
}

但是代码会报错,出现如下提示

1
2
3
4
5
6
7
error[E0597]: `x` does not live long enough     
let s: Vec<&str> = x.split(',').collect();
^ borrowed value does not live long enough
hm.insert(s[0],s[1]);
-- borrow later used here
}
- `x` dropped here while still borrowed

从我的角度看,虽然x是被借用的数据,但是通过split和collect得到的s,里面记录的 &str指向的字符串应该是新分配的,跟x没有关系才对,所以我觉得就算是报错,也应该出在s的生命周期上。但实际上,通过对文档的阅读,我才知道vector里存放的值实际是x的切片,每一次循环结束后,x的生命周期就结束了,切片数据也同样不复存在。

An iterator over substrings of this string slice, separated by characters matched by a pattern.

同时由于&str只是字符串的借用,因此无法转移所有权, 需要将&str转成String才行。

1
2
3
4
5
6
7
8
9
10
11
12
let k = s[0].to_string();
let v = s[1].to_string();

/* to_string方法的源代码如下
// 其中to_owned函数通过clone从借用数据创建具有所有权的数据
impl ToString for String {
#[inline]
fn to_string(&self) -> String {
self.to_owned()
}
}
*/

如果用的csv库进行解析,也同样需要考虑这个问题。下面代码会因为record的生命周期出错,也需要将其转成String。

1
2
3
4
5
6
7
8
9
10
11
let mut rdr = ReaderBuilder::new().from_path(csv_file)?;

for result in rdr.records(){
if let Ok(record) = result{
let k = record.get(0).unwrap();
let v= record.get(1).unwrap();
//println!("{}:{}", cluster, cell_name);
let vec = hm.entry(k).or_insert(Vec::new());
vec.push(v);
}
}

除开上面两个主要问题外,还有一些小问题,也简单记录一下

  • Rust并不总能推导出变量的类型,比如说由于有很多集合类型都实现了collect方法,就导致Rust不知道你调用collect到底要什么

  • Rust的vector添加单个元素是push, append是合并另外一个vector

  • HashMap的get函数只能返回key所对应的value的引用,改引用无法修改,get_mut返回才是可变引用

  • HashMap除了使用get_mut函数获取可变引用进行修改外,还可以使用entry函数。

最终代码如下:

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
use std::error::Error;
use std::env;
use std::collections::HashMap;
use csv::Reader;

fn main() -> Result<(), Box<dyn Error>>{
let args: Vec<String> = env::args().collect();
let csv_file = &args[1];

let mut hm= HashMap::new();

let mut rdr = Reader::from_path(csv_file)?;
for result in rdr.records(){
if let Ok(record) = result{
let k = record.get(0).unwrap();
let v= record.get(1).unwrap();
//println!("{}:{}", cluster, cell_name);
let vec = hm.entry(k).or_insert(Vec::new());
vec.push(v);
}
}

for (key, value) in &hm{
println!("{}: {:?}", key, value);
}


Ok(())
}

Rust第三课:模式匹配match语法

我目前阅读Rust代码的一个不适应点就是Rust中的模式匹配match语法, 以及两个语法糖if letwhile let,为了强化自己的记忆,于是便有了第三课。

所谓的模式匹配,就是对于不同情况采取不同的处理方法,

让我们想象一个应用场景,将一段DNA序列进行互补,即A->T, T->A, C->G, G->。

对于Python而言,我们可以通过字典(dict)进行翻译

1
2
3
4
5
6
7
dna = "ATCG"
trans_dict = {'A' : 'T', 'T': 'A', 'C':'G', 'G':'C'}
translated = ""
for base in dna:
translated+=(trans_dict[base])

print(translated)

如果是C语言的话,我会用switch语法,如下

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
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

int main(){
char *dna = "ATCG";
int dna_len = strlen(dna);
char *translated = (char*)malloc( sizeof(char) * dna_len);

for (int i = 0; i < dna_len; i++){
char base = dna[i];
switch (base) {
case 'A':
translated[i] = 'T';
break;
case 'T':
translated[i] = 'A';
break;
case 'C':
translated[i] = 'G';
break;
case 'G':
translated[i] = 'C';
break;
default:
translated[i] = 'N';
break;
}
}
printf("%s\n", translated);
return 0;
}

Rust提供了match用于模式匹配。以一个碱基为例,我们输入’A’, 希望翻译成’T’, 下面是非常粗糙的代码

1
2
3
4
5
6
7
8
9
10
fn main(){
let base = 'A';
let trans = match base {
'A' => 'T',
'T' => 'A',
'C' => 'G',
'G' => 'C',
_ => 'N',
};
println!("{}", trans);

其中match语法描述如下

1
2
3
4
5
6
match 变量名 {
匹配类型1 => 代码1,
匹配类型2 => 代码2,
....
_ => 其他
}

match要求我们考虑所有的情况, 而一个utf-8字符肯定不只是ATCG这四个,因此其他情况需要用_指代,否则会报错。

对于一段序列, 我们可以写一个循环,直接在里面嵌套一个match语法。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
fn main(){
let dna = "ATCG".to_string();
let mut translated = String::new();

for i in dna.chars() {
translated.push(
match i {
'A' => 'T',
'T' => 'A',
'C' => 'G',
'G' => 'C',
_ => 'N',
}
)
}
println!("{}", translated );

}

或者考虑将match匹配的部分封装到一个函数中

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
fn translate(base: char ) -> char {

match base {
'A' => 'T',
'T' => 'A',
'C' => 'G',
'G' => 'C',
_ => 'N',
}
}

fn main(){

let dna = "ATCG".to_string();
let mut translated = String::new();

for base in dna.chars() {
translated.push(translate(base));

}
println!("{}", translated );
}

但是match不仅

Rust第二课:为什么我的Rust比Python慢!

我的Rust第一课, 我写了一个程序对fasta中的ATCG进行计数。后面,我就想到一个非常常见的需求,对文件进行读取,统计行数,类似于 wc -l

下面是我写的第一个版本的代码, 我命名为myRead.rs

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
use std::io::BufReader;
use std::fs::File;
use std::env;
use std::io::BufRead;

fn main() -> std::io::Result<()> {

let args: Vec<String> = env::args().collect();
let filename = &args[1];

let f = File::open(filename)?;
let reader = BufReader::new(f);
let mut line_num = 0;

for _line in reader.lines() {
line_num += 1

}
println!("{}", line_num);
Ok(())
}

然后用rustc进行编译

1
rustc myRead.rs

接着我用一个记录cdna的fasta文件(约300MB),进行测试,耗时约5.6秒

1
2
time ./myRead Homo_sapiens.GRCh38.cdna.all.fa
# 5.50s user 0.13s system 99% cpu 5.647 total

同样,我写了5行python脚本进行比较

1
2
3
4
5
6
import sys
count = 0
for line in open(sys.argv[1]):
count += 1
print(f"line number {count}")

python代码不到1秒就完成了任务

1
2
 time python ./read_file.py Homo_sapiens.GRCh38.cdna.all.fa
# 0.75s user 0.16s system 99% cpu 0.904 total

看到这个结果我直接震惊. Rust的运行速度居然比Python慢了6倍左右。经过高强度的检索,终于被我找到了靠谱的答案, BufReader 100x slower than Python — am I doing something wrong?

总结下原因就是

  1. 默认的优化不行,对于rustc 需要设置 -C opt-level=2 或者等价的 -O, 对于cargo则是设置 --release
  2. .lines() 会为每一行都重新分配内存,因此不仅仅是处理UTF-8的问题。

根据第一个建议, 我重新用rustc编译了代码, 速度直接超过了Python

1
2
3
rustc -O ./myRead.rs
time ./myRead Homo_sapiens.GRCh38.cdna.all.fa
# 0.51s user 0.12s system 72% cpu 0.874 total

根据第二个建议,重新写了如下代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
use std::fs::File;
use std::env;
use std::io::Read;

fn main() {
let args: Vec<String> = env::args().collect();
let filename = &args[1];
let mut file = File::open(filename).unwrap();
let mut lines = 0;
let mut buf = [0u8; 4096*32];
while let Ok(num_bytes) = file.read(&mut buf) {
if num_bytes == 0 { break; }
lines += buf[..num_bytes].iter().filter(|&&byte| byte == b'\n').count();
}
println!("line number {}", lines);

}

使用 rustc -O编译后,运行速度又提升了2倍。

实际上,对于我这个Rust初学者,只需要记住Rust程序在编译的时候要设置优化参数, 进一步的优化代码,一时半会我还看不懂。

bwa索引文件amb/ann/pac格式说明

高通量数据比对讲究的就是一个快和准,因此大部分软件都是用C语言实现。BWA是目前基因组序列比对最常用的工具,由于自我感觉已经入门C语言,为了提高自己的水平,因此开始从源码角度学习李恒大神开发的BWA工具。

使用 bwa index 建立索引后,会得到以下后缀的文件, .amb, .ann, .pac, .bwt, .sa. 我们通常也不在乎这些文件是什么,除非你像我一样,想搞懂bwa建立索引的具体过程。这次我们先来介绍.amb, .ann和.pac这三个文件是什么。

先说结论,amb是ambiguous的缩写,也就是模棱两可(兼并碱基)的意思,也就是除了ATCG/atcg以外的字符. amb和ann用来记录基因组中除了ATCG以外碱基的信息。而pac文件则是碱基信息高度压缩。

碱基表示方法

接下来的介绍,主要涉及到两个文件,bntseq.h和bntseq.c, 这里面的代码就是用于将fasta转成amb, ann和pac这三个文件。

bntseq.h的开头,就定义了对应输出文件的三种数据结构, bntann1_t, bntamb1_tbntseq_t

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
typedef struct {
int64_t offset; //染色体偏移量
int32_t len; //染色体长度
int32_t n_ambs; //多少个模糊碱基
uint32_t gi;
int32_t is_alt;
char *name, *anno; //染色体名名和注释
} bntann1_t;

typedef struct {
int64_t offset; //偏移量
int32_t len; //长度
char amb; //模式碱基的字符
} bntamb1_t;

typedef struct {
int64_t l_pac;
int32_t n_seqs; //基因组序列数
uint32_t seed; //随机数种子
bntann1_t *anns; // n_seqs elements
int32_t n_holes; //染色体有多少个空缺
bntamb1_t *ambs; // n_holes elements
FILE *fp_pac;
} bntseq_t;

前两种数据结构能够保存成文本,我们以拟南芥基因组为例。ann文件前几行如下,它记录每个染色体信息,以及里面amb数目

1
2
3
4
5
6
239335500 7 11 # bns->l_pac, bns->n_seqs, bns->seed)
0 Chr1 (null) # p->gi, p->name, p->anno
0 30427671 363 # p->offset, p->len, p->n_ambs
0 Chr2 (null)
30427671 19698289 58
...

amb文件前几行如下,它 记录确切的amb位置

1
2
3
4
239335500 7 563 # bns->l_pac, bns->n_seqs, bns->n_holes
12192623 1 Y # p->offset, p->len, p->amb
13201252 100 N # p->offset, p->len, p->amb
...

第三种数据结构,bntseq_t,则是用于打包将前两者的信息进行了整合,同时记录pac文件的文件指针地址。

对于pac文件,他通过哈希映射的方法(哈希函数如下),把一个比较大的位置映射到比较小的位置中,保证了最终序列文件大小低于原始文件。

1
2
#define _set_pac(pac, l, c) ((pac)[(l)>>2] |= (c)<<((~(l)&3)<<1))
#define _get_pac(pac, l) ((pac)[(l)>>2]>>((~(l)&3)<<1)&3)

总结一下,amb,ann,pac这三个文件主要是用于压缩原始的参考基因组序列,将碱基信息分为两部分,ATCG和其他兼并碱基的序列。amb记录兼并碱基的具体位置,ann记录参考基因组中染色体信息,每条染色体的偏移信息和长度信息。我们原本可以去掉染色体名,将所有的序列保存到一个文件中,但是这会导致文件过大,不利于后续的读取操作。因此,李恒采用了一种计算方式,使得只用四分之一的体质保存原来的信息。