生物信息-python 编程实例（1）

我们用python处理小的文本文件时一般使用文件对象的read()、readline()和readlines()方法就可以了，但是，NGS测序质量文件一般都比较大，一个文件可能几G，几十G，上面的方法是把文件内容一次性全部读入内存，对内存的消耗是很大的，甚至会造成内存溢出，那样处理就不适合了。

我们可以把文件当作迭代器，用可迭代对象的next()方法，减轻对内存的消耗。

当迭代完最后⼀项数据之后，再次调⽤next()函数会抛出StopIteration的异常，说明所有数据都已迭代完，不能再执⾏next()方法了，我们捕捉异常，从而完成操作，解决了读取和处理大文件的问题。

FASTQ文件格式简要说明：

每四行一个单位，

1.第一行以’@’开头，表示序列标识以及相关的描述信息；

2.第二行是序列信息；

3.第三行以‘+’开头，后面是序列标示符、描述信息，或者留空；

4.第四行，是质量评分，和第二行的序列相对应，每一个序列都有一个质量评分。字符数与第二行相同，每个字符代表对应第二行碱基的质量。

如果我们要解析双端测序文件，可以根据其特点进行构造解析函数：

代码要点注释：

第2行定义函数, 第一和第二个参数分别为read1、read2的两个文件对象，

5-8行表示每次循环将read1和read2的每4行内容分别放入列表第0项和第1项，

第9行表示对数据进行的处理方式，我们以打印出来为例，实际环境中可以根据需要进行各种操作，如统计测序信息、质控、去N和去接头等。

第10和第11行表示，遇到异常终止循环 ！

这样设计，就解决了大文件对内存消耗过大的问题 ！

我们可以对以上代码稍加改动，就可以成为真正的NGS测序数据质控软件啦！ 请看下面代码示例：

仅仅对第一个例子的代码增加了简单的判断，就可以打印出符合条件的read的位置、read1和read2中分别N的个数。。

我测试的数据是 SRA ERR1138631， read1和read2 分别取前5000条数据 ！

实际生产环境中，仅仅打印N的个数实在没啥大的用处，我们要把符合条件的read筛选出来并保存到文件中，实现起来也不难，读者可以先试下。

到今天为止，咱也算会写质控软件啦 ！哈哈，虽然有些简陋！

加油 ！ 后面的内容更精彩 ！