0

我是生物信息学领域的新手，也是 Python 初学者。

我正在使用一个名为 Ruffus 的框架来运行一个管道，使文件通过工具导航以获得有趣的输出。

在我的管道的一部分（第三步）中，我需要做一些解决方法来让 Ruffus 考虑我的文件：使用的工具几乎可以任意更改输出名称。

我的问题如下：我必须第一次运行管道以使其运行到第三步，它说每个文件都已完成。然后当我重新运行时，它从第 3 次到最后。我认为问题是我使用了与特定模式有关的 glob，该 glob 函数应该在我的脚本开始时运行，但我正在寻找的文件 (.*avinput) 尚未创建。然后，一旦我重新运行它，它们就可以了，然后该过程可以继续进行。但这是不现实的。

我如何“强制”这样的 var 等到特定时刻运行，或者在达到某个步骤时重新运行它？

最好的，

from ruffus import *
import sys
import subprocess
import os
import re
import glob
import time
import termcolor



#====================================
#
#Input area
#
#
#====================================

#Create a method to get all the files from the folder
input_files = (["/genetics/ongoing_work/pipeline/78_14_TR_p_BWA_mem_PIC_conv_ord_mdup_group_recal_rawvariants_Joint_INDEL_recalibrated_extracted_filter_appl.vcf",
               "/genetics/ongoing_work/pipeline/78_14_TR_p_BWA_mem_PIC_conv_ord_mdup_group_recal_rawvariants_Joint_SNP_recalibrated_extracted_filter_appl.vcf"])
# DIRECTORIES
outdir = '/genetics/ongoing_work/pipeline'
annovar_db='/genetics/apps/annovar/humandb/'

#SOFTWARES
gatk = '/genetics/apps/gatk3.8/GenomeAnalysisTK.jar'
convert2annovar='/genetics/apps/annovar/convert2annovar.pl'
annotate_variation='/genetics/apps/annovar/annotate_variation.pl'
table_annovar='/genetics/apps/annovar/table_annovar.pl'

#FILES
reference_genome = '/genetics/reference/ucsc.hg19.fasta'


# ---------- 1st step -----------
#  Filter Low Quality Genotypes
#  Filter Low Quality Genotypes
#  GATK
#
# ---------------------------

@transform(input_files, formatter("([^/]+).vcf$"), os.path.join(outdir, "{1[0]}_flqg.vcf"))

def LowQualityGenotype(input_file, output_file):
    print(termcolor.colored(' ---------- 1st step ----------- \n  \
  Filter Low Quality Genotypes \n \
  Filter Low Quality Genotypes \n  \
  GATK \n  \
 \n  \
 --------------------------- \n','red'))

    subprocess.run("java -jar {tool} \
    -T VariantFiltration \
    -R {ref} \
    -V {i_file} \
    -G_filter 'GQ < 20.0' \
    -G_filterName 'LowQG' \
    -o {o_file}".format(
        o_file=output_file,
        i_file=input_file,
        tool=gatk,
        ref=reference_genome
    ), shell=True)


# ---------- 2nd step -----------
#  Convert to Avinput
#  ANNOVAR
# UPDATE TO FIX
# ---------------------------


@subdivide(LowQualityGenotype, formatter("([^/]+).vcf$"), os.path.join(outdir, "{1[0]}.avinput"))

def ConvertAvinput(input_file, output_files):
    subprocess.run("perl {tool} \
    -format vcf4 -allsample\
    -withzyg \
    -includeinfo \
    {i_file} \
    -outfile {o_file}".format(
        o_file=output_files,
        i_file=input_file,
        tool=convert2annovar,
    ), shell=True)
    #modifying outputs because 2 gives 4 with the samples names in it

    file_list_avinput = glob.glob(outdir+'/*.avinput')
    print('Impression du glob')
    print(file_list_avinput)
    #use of glob to get all the files created
    for i in file_list_avinput:
        if i:
            os.rename(i, re.sub('(.avinput)(?:\.)(.{5})(?:.*)',"_sample_\\2\\1",i))

    time.sleep(3)





# ---------- 3rd step -----------
#  Filter 1000G
#  ANNOVAR
#
# ---------------------------
pre1000g_files = glob.glob(outdir + '/*.avinput')
print('LIGNES POUR PRE1000G')
print(pre1000g_files)


@follows(ConvertAvinput)

@subdivide(pre1000g_files, formatter("([^/]+).avinput"),
           os.path.join(outdir, "{1[0]}_1000G_filtered"))

def Filter1000G(input_file, output_files):
    print('LIGNES POUR PRE1000G')
    print(pre1000g_files)
    subprocess.run("perl {tool} \
    -filter \
    -dbtype 1000g2015aug_eur \
    -buildver hg19 \
    -maf {maf_freq} \
    -out {o_file} \
    {i_file} \
    {annodb} \
    -outfile {o_file}".format(
        o_file=output_files,
        i_file=input_file,
        maf_freq=0.03,
        annodb=annovar_db,
        tool=annotate_variation
    ), shell=True)

    os.rename(output_files+'.hg19_EUR.sites.2015_08_filtered',output_files[0:-8]+'filtered')
    os.rename(output_files + '.hg19_EUR.sites.2015_08_dropped', output_files[0:-7]+'_dropped')


# ---------- 4th step -----------
#  Filter InHouse
#  ANNOVAR
#
# ---------------------------


@transform(Filter1000G,  formatter("([^/]+)_filtered"),
                                        os.path.join(outdir, "{1[0]}_Norge_filtered"))
def FilterInHouse(input_file, output_file):

    subprocess.run("perl {tool} \
    -filter \
    -dbtype generic \
    -genericdbfile nimblegen_v3_and_v2_combined_jul2015.singleallelic.avinput.AC_gt_5 \
    -buildver hg19 \
    -out {o_file} \
    {i_file} \
    {annodb} \
    -outfile {o_file}".format(
        o_file=output_file,
        i_file=input_file,
        maf_freq=0.005,
        annodb=annovar_db,
        tool=annotate_variation
    ), shell=True)
    os.rename(output_file+'.hg19_generic_filtered',output_file[0:-8]+'filtered')
    os.rename(output_file + '.hg19_generic_dropped', output_file[0:-7]+'_dropped')
pipeline_run(FilterInHouse)
4

0 回答 0