生物信息--python 编程实例(1)

-回复 -浏览
楼主 2019-04-14 12:12:27
举报 只看此人 收藏本贴 楼主

                        

      

    我们用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筛选出来并保存到文件中,实现起来也不难,读者可以先试下。

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

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

我要推荐
转发到