python 随机抽取Fastq文件

参考 http://pythonforbiologists.com/index.php/randomly-sampling-reads-from-a-fastq-file/

最近要做一个二代测序的模拟，所以网上找了个小脚本，做了些注释，希望能够帮助大家。

from __future__ import division
import random

number_to_sample = 3000000
number_of_replicates = 10

#计算行数 
with open("test.fastq") as input:
    num_lines = sum([1 for line in input])

# 全部reads数目
total_records = int(num_lines / 4)
print("sampling " + str(number_to_sample) + " out of " + str(total_records) + " records")

# 可以设置多个输出文件
output_files = []
output_sequence_sets = []

#打开10个文件
for i in range(number_of_replicates):
    output_files.append(open("sample.fastq." + str(i), "w"))

# 随机抽取    output_sequence_sets.append(set(random.sample(xrange(total_records + 1), number_to_sample)))

# 一次读取4行
record_number = 0
with open("test.fastq") as input:
        for line1 in input:
            line2 = input.next()
            line3 = input.next()
            line4 = input.next()
            for i, output in enumerate(output_files):
                if record_number in output_sequence_sets[i]:
                        output.write(line1)
                        output.write(line2)
                        output.write(line3)
                        output.write(line4)
                record_number += 1

for output in output_files:
    output.close()