首页 > 解决方案 > 从 fastq 文件列表中计算读取和碱基

问题描述

我使用 Trimmomatic 修剪了我的 Illumina 短读数,正向和反向。Trimmomatic 的输出是:paired_1 - unpaired_1 和paired_2 - unpaired_2.fastq.gz 文件。我想通过计算目录中每个文件的读取次数和基数来了解修剪的影响有多大。我编写了一个脚本来计算目录中每个文件的碱基数和读取数;但是,我在if __name__=='__main__'. 当我执行 for 循环时,我不知道将要运行的文件的顺序,我怎样才能让它按照我从屏幕上看到的顺序调用文件?此外,我还需要帮助来纠正脚本,因为我没有得到任何标准输出。预先感谢您的帮助。

#!/usr/bin/env python

from sys import argv 
import os 

def get_num_bases(file_content):
    total = []
    for linenumber, line in enumerate(file_content):
        mod=linenumber%4

        if mod==0: 
            ID = line.strip()[1:]
            #print(ID)

        if mod==1: 

            seq = line.strip()
            counting = 0
            counting += seq.count("T")+ seq.count("A") + seq.count("C") + seq.count("G")
            total.append(counting)
            allbases = sum(total)
            print("Number of bases are: " , allbases)

def get_num_reads(file_content):

    total_1 = []

    for line in file_content: 
        num_reads = 0 
        num_reads += content.count(line)  
        total_1.append(num_reads)
        print("Number of reads are: ", sum(total_1)/int(4))


if __name__=='__main__':
    path = os.getcwd()
    dir_files = os.listdir(path)
    list_files = []
    for file in dir_files:
        if file.endswith("fastq.gz"):
            if file not in list_files:
                file_content = open(file, "r").readlines() 
                list_files.append(file)
                print("This is the filename: ", file, get_num_bases(file_content), get_num_reads(file_content))

标签: python

解决方案


推荐阅读