shell截取文件行数_截取fastq文件reads的长度

随着测序量越来越多，分析的个性化需求也不断增加，单靠别人的写好的程序已经不能满足我们的要求，万事靠自己才是王道。接下来一段时间我会不定期讲解一些用python分析数据的简单应用。今天是用python截取fastq文件中reads的特定长度，现在就开始啦！
每条reads有四行组成，第一行是物理地址，第二行是测序的序列，第三行是个加号，第四行是测序质量。下面是两条reads：

@RS-E00182:745:S3H5HTTLXX:4:1341:4340:2901 1:N:0:CGTGGTATTAAAGATTTCTCAAATATTTAATAAATGAGGCAAAACAATCATGAAAATATAAACTGCTACTATAGATTTTTCATCTTTCATTTTATGCATTTCCTCTCAACCATGTAATTTCCTGTTAAGGAATGTTGATAAATTTTTATATTTATAAT+AAAFFKKKKAFFKAKAKKKFFKKFKKAFFK@RS-E00182:745:S3H5HTTLXX:4:1341:4543:2901 1:N:0:CGTGGTATAGGCGGTCAAAAAGCCGCCCCCGGTGGCATTCAAGGTGATGTGCTTGCTACCGATAACAATACTGTAGGCATGGGTGATGCTGGTATTAAATCTGCCATTCAAGGCTCTAATGTTCCTAACCCCGATGAGGCCGTCCCTAGTTTTGTTTC+<7,

测序的reads长度是150bp，太长，map到基因组的效率会下降，我们准备截取50bp去map。

现在就开始编程啦！

#!/usr/bin/env python #建立python环境fq = open('test.fq', 'r')        #打开文件i = 0for line in fq:                  #for循环来取值    i = i + 1                    #记录reads的行数    if i%4 == 1:                 #行数除于4余数是1代表取到第一行reads的地址        print(line.strip())      #输出reads的地址, strip()函数是去除转行符    elif i%4 == 2:               #行数除于4余数是2代表取到第二行reads的序列        print(line[0:50])        #取出前50bp序列    elif i%4 == 3:               #行数除于4余数是3代表取到第三行的加号        print(line.strip())    elif i%4 == 0:              #行数除于4余数是0代表取到第四行的reads的质量        print(line[0:50])       #取出前50bp测序质量

结果如下：

@RS-E00182:745:S3H5HTTLXX:4:1341:4340:2901 1:N:0:CGTGGTATTAAAGATTTCTCAAATATTTAATAAATGAGGCAAAACAATCATGAAAATA+AAAFFKKKKAFFKAKAKKKFFKKFKKAFFK@RS-E00182:745:S3H5HTTLXX:4:1341:4543:2901 1:N:0:CGTGGTATAGGCGGTCAAAAAGCCGCCCCCGGTGGCATTCAAGGTGATGTGCTTGCTA+<

就像这样，map的速度和map率都会增加。如果你有什么要求和错误，欢迎指出！