0

我有两条路

path1 = "/home/x/nearline"
path2 = "/home/x/sge_jobs_output"

在 path1 中,我有一堆 fastq 文件:

ERR001268_1.recal.fastq.gz
ERR001268_2.recal.fastq.gz
ERR001269_1.recal.fastq.gz
ERR001269_2.recal.fastq.gz
.............

在path2中,我有很多.txt对应于path1中的fastq文件:

ERR001268_1.txt
ERR001268_2.txt
ERR001269_1.txt
ERR001269_2.txt
.............

现在我已经制作了脚本来从 path1 中的 fastq 文件计算 fastq_seq_num,见下文:

for file in os.listdir(path1):
  if re.match('.*\.recal.fastq.gz', file):
    fullpath1 = os.path.join(path1, file)
#To calculate the sequence number in fastq.gz files  
    result = commands.getoutput('zcat ' + fullpath1 + ' |wc -l')
    fastq_seq_num = int(result)/4.0
  print file,fastq_seq_num 

并且还从 path2 中的 .txt 文件计算 num_seq_processed_sai,见下文:

for file in os.listdir(path2):
  if re.match('.*\.txt', file):
      fullpath2 = os.path.join(path2, file)
#To calculate how many sequences have been processed in .sai file
      linelist = open (fullpath2,'r').readlines
      lastline = linelist[len(linelist)-1]
      num_seq_processed_sai = lastline.split(']')[1].split()[0]
  print file,num_seq_processed_sai

好的,现在我的问题是:我想创建一个循环,在其中计算 path1 中第一个 fastq 文件的 fastq_seq_num;然后计算path2中FIRST txt文件的num_seq_processed;然后比较这两个数字;然后结束循环。然后第二个循环开始......我怎样才能设计一些循环来实现这一点?谢谢!!!

4

2 回答 2

2

两个目录中的文件数是否保证相同?如果是这样,您可以使用该zip函数来完成此操作:

for fastqFile, txtFile in zip(glob.glob(path1+'/*.recal.fastq.gz'), glob.glob(path2+'/*.txt')):
    result = commands.getoutput('zcat ' + fastqFile + ' |wc -l')
    fastq_seq_num = int(result)/4.0

    lastline = linelist[-1]
    num_seq_processed_sai = lastline.split(']')[1].split()[0]


    print fastqFile, fastq_seq_num 
    print txtFile, num_seq_processed_sai

编辑:作为旁注,使用glob.glob()几乎总是比手动过滤输出更可取,os.listdir()并且该单词file是 Python 中的内置类型,您永远不应该将其用作变量名。此外,要到达 a 的最后一项list,您只需使用 访问它listName[-1]。用它来访问它listName[len(listName)-1]是不合常理的。

于 2011-06-26T23:42:29.367 回答
2

迭代 fastq 文件,检查对应的 .txt 文件,如果存在,运行你的处理并比较输出

import commands
import glob
from os import path

dir1 = '/home/x/nearline'
dir2 = '/home/x/sge_jobs_output'

for filepath in glob.glob(path.join(dir1, '*.recal.fastq.gz')):
    filename = path.basename(filepath)
    job_id = filename.split('.', 1)[0]

    ## Look for corresponding .txt file
    txt_filepath = path.join(dir2, '%s.txt' % job_id)
    ## Fail early if corresponding .txt file is missing
    if not path.isfile(txt_filepath):
        print('Missing %s for %s' % (txt_filepath, filepath))
        continue

    ## Both exist, process each
    ## This is from your code snippet
    result = commands.getoutput('zcat ' + fullpath1 + ' |wc -l')
    fastq_seq_num = int(result)/4.0

    linelist = open(txt_filepath).readlines()
    lastline = linelist[len(linelist)-1]
    num_seq_processed_sai = lastline.split(']')[1].split()[0]

    if fastq_seq_num == num_seq_processed_sai:
        print "Sequence numbers match (%d : %d)" % (fastq_seq_num, num_seq_processed_sai)
    else:
        print "Sequence numbers do not match (%d : %d)" % (fastq_seq_num, num_seq_processed_sai)

我建议使用glob.glob()来列出你的文件(就像我在这个片段中所做的那样)。

我还建议替换commandssubprocess,它是较新的,我认为命令已被弃用。

此外,读取 .txt 文件的所有行以获取最后一行是低效的,如果文件很大,可能会导致内存问题。考虑调用tail以获取最后一行(我认为只有 *nix)。

于 2011-06-26T23:51:10.787 回答