导读

做最基本的fastq序列统计。这里的输入数据是一个上游文件，一个下游文件。

一、新建文件

vi summary_fastq.py

二、写函数

• 函数1：查询文件

思路：
1 for os.listdir 遍历文件
2 if in判断是否含关键字符串

• 函数2：统计fastq序列数、碱基数、GC%

思路：
1 for in 逐行遍历fastq文件
2 行数%4==2是序列行
3 累加序列数、碱基数、GC数【format格式化GB单位】，记录序列长度，求GC%和MaxLength、MinLength
4 加到列表return

#!/usr/bin/env python3
import os, re, sys

def search(path, key):
    result = []
    for filename in os.listdir(path):
        file = os.path.join(path, filename)
        if os.path.isfile(file) and key in filename:
            result.append(file)
    return(result)

def reads(file):
    with open(file) as f:
        result = []
        lengths = []
        line_num = 0
        gc_num = 0
        base_num = 0
        read_num = 0
        for line in f:
            line_num += 1
            if line_num%4 ==2:
                read_num += 1
                gc_num += line.count("G") + line.count("g") + line.count("C") + line.count("c")
                base_num += len(line)
                lengths.append(len(line))
                
        print("\033[32m---------------read_num, base_num, gc_percent, max_length, min_length---------------\033[0m")
        gc_percent = format(gc_num/base_num*100, '0.2f')
        
        result.append(read_num)
        result.append(base_num)
        result.append(gc_percent)
        result.append(max(lengths))
        result.append(min(lengths))
        return(result)

三、主函数

思路：
1 合并raw clean数据r1 r2
2 reads函数处理
3 格式化输出
4 判断并删除中间文件

if __name__ == '__main__':
    print("\033[32m---------------merge rawdata---------------\033[0m")
    os.system("cat " + search("rawdata", "R")[0] + " " + search("rawdata", "R")[1] + " > merge_raw.fq")
    print("\033[32m---------------merge cleandata---------------\033[0m")
    os.system("cat " + search("Result/kneaddata", "trimmed.1")[0] + " " + search("Result/kneaddata", "trimmed.2")[0] + " > merge_clean.fq")
    
    with open("Result/kneaddata/fastq_summary.txt", 'w') as o:
        o.write("\tInsertSize(bp)\tSeqStrategy\tReads(#)\tBase(GB)\tGC(%)\tMaxLength\tMinLength\n")
        if os.path.isfile("./merge_raw.fq") and os.path.isfile("./merge_raw.fq"):
            f1 = "./merge_raw.fq"
            reads_info = reads(f1)
            read_num = reads_info[0]
            base_num = format(reads_info[1]/10**9, '0.2f')
            gc_percent = reads_info[2]
            max_length = reads_info[3]
            min_length = reads_info[4]
            #print(read_num, base_num, gc_percent, max_length, min_length)
            o.write("RawData\t350\t(150:150)\t{}\t{}\t{}\t{}\t{}\n".format(read_num, base_num, gc_percent, max_length, min_length))
            
            f2 = "./merge_clean.fq"
            reads_info = reads(f2)
            read_num = reads_info[0]
            base_num = format(reads_info[1]/10**9, '0.2f')
            gc_percent = reads_info[2]
            max_length = reads_info[3]
            min_length = reads_info[4]
            #print(read_num, base_num, gc_percent, max_length, min_length)
            o.write("CleanData\t350\t(150:150)\t{}\t{}\t{}\t{}\t{}".format(read_num, base_num, gc_percent, max_length, min_length))
        if os.path.isfile("Result/kneaddata/fastq_summary.txt"):
            os.remove("./merge_raw.fq")
            os.remove("./merge_clean.fq")
        else:
            print("\033[31m no fastq_summary.txt! please check!\033[0m")

四、结果

python3 summary_fastq.py

    InsertSize(bp)  SeqStrategy Reads(#)    Base(GB)    GC(%)   MaxLength   MinLength
RawData 350 (150:150)   16493814    2.49    32.90   151 151
CleanData   350 (150:150)   14872620    2.14    32.82   151 51