python 从fastq文件中挑选出序列长度在规定范围的序列

要求：

FASTQ文件长度过滤：编写一个程序，读取FASTQ文件，过滤掉碱基序列长度在60-80之外的序列，将长度在60-80之内的序列输出到结果文件中。

FASTQ格式文件如下：每四行表示一个测序序列，第二行是碱基序列。

@M2,HWI-7001455:326:HGTFVADXX:1:1101:2666:105721:N:0:GCGCTA

TTTAGTTTTGTAGTAATTGTTTGTAGTAAAATTTGTATTAGTTTTTTTATTTGTA

+

CCCFFFDDHHFHHIJJJJJHHIJJJJIIGEIJIJJIIJJJJHIJJJJJIJIJIHI

编程要求如下:

1)      程序命名为fq_filter.py

2)      程序采用optparse从命令行输入，参数共有3个，-infile，-outfile，-h；

3)      其中-infile用于接收输入的FASTQ文件名；

4)      -outfile用于给出符合条件的结果文件名，输出结果格式要与输入的FASTQ格式相同。

5)      -h用于给出程序的使用说明；

fq_filter.py

#!/usr/bin/python
from optparse import OptionParser


parser = OptionParser()
parser.add_option("--infile", dest="infile", help="give a fasta file to me", metavar="FILE")
parser.add_option("--outfile",  dest="outfile", help="the name of oufput file [fasta]", metavar="FILE")
(options, args) = parser.parse_args()


fqfile=file(options.infile,'r')
fqfilter=file(options.outfile,'w')


i=0
for line in fqfile:
        #print i,i%4,line
        if i%4==0:
                seqID=line.strip("\n")
        elif i%4==1:
                sequence=line.strip("\n")
        elif i%4==2:
                mark=line.strip("\n")
        elif i%4==3:
                if len(sequence)>=60 and len(sequence)<=80:
                        qual=line.strip("\n")
                        fqfilter.write(seqID+"\n"+sequence+"\n"+mark+"\n"+qual+"\n")
        i+=1


fqfile.close()
fqfilter.close()

运行命令：

py fq_filter.py --infile test.fq --outfile out.fa