Python读Fastq新姿势

分析二代测序数据时免不了要对fastq文件进行操作，然而也会出现想进行的操作并没有现成的工具可以实现的情况，这时候就需要自己写一个小脚本来读取fastq文件并进行操作。

其实用Python读取fastq文件的逻辑也很简单，根据fastq每四行为一个read的特点，边读边记行数，行数除4余2的行即位序列所在的行。这种方法看起来很naive，不过好像也没有更优雅的方式，直到我发现了mappy。

mappy简介

mappy是minimap2的Python版本，同样也是出自Li Heng，你可以使用它在Python程序直接实现和minimap2同样性能的比对工具，而不在需要从命令行调用。不要问为什么这么强，问就是Cython。

而这里我们需要用到的只是mappy的一个小小的功能：

可以看到mappy读取fastq只需要一行：mp.fastx_read('./test.fq',read_comment=True)，然后返回一个迭代器，内容是一个元组，read ID，序列包括质量等信息都在里面。

mappy性能对比

这里通过一个小任务来比较mappy和我以前使用的方法的区别（其实也是Cython和Python的区别）

任务是将一个有1千万条read的fastq文件截取前50bp写进新fastq，read ID不变。

相比之下，不仅代码少了很多，而且时间缩短了约25%。

换成gzip压缩文件后：

mappy的速度优势更明显了。

今天也是羡慕C的速度的一天呢😐。