paramsweep.snk

import os

#Load Config File
configfile: "paramsweep.yaml"

runs = config['runs']

#For each run we detect the barcodes it contains (individual samples)
barcodes = {}

#It is possible to use a user defined subset here as well
if config['useSubsetOfBarcodes']:
    barcodes = config['barcodes']
    print("Using the following subset of barcodes: {}".format(barcodes))
else:
    for run in runs:
        #We use the medaka output here to determine the existing barcodes, this is arbritrary, any other file could be chosen for this purpose
        barcodes[run] = glob_wildcards('data/input/'+run+'/barcode{barcode}.medaka.'+config['vcf_suffix']).barcode


### REMOVE LATER (THIS SECTION IS ONLY USED FOR EVALUATION OF THE PIPELINE AND IRRELEVANT FOR PRODUCTIVE USE)
gisaidMapping = {}
gisaidMappingInverse = {}

with open('data/input/mappingRunsGisaid.csv', 'r') as infile:
    lines = infile.read().splitlines()
    for l in lines:
        data = l.split()
        run = data[0]
        barcode = data[1]
        if len(data) != 2: #entry in the table
            gisaidID = data[2]
            if not run in gisaidMapping:
                gisaidMapping[run] = {}
            gisaidMapping[run][barcode] = gisaidID
            gisaidMappingInverse[gisaidID] = run+'_'+barcode

def getGisaidFile(run,barcode):
    if run in gisaidMapping:
        if barcode in gisaidMapping[run]:
            return gisaidMapping[run][barcode]
    return None

for run in runs:
    for barcode in barcodes[run]:
        gisaidid = getGisaidFile(run,barcode)

        if gisaidid:
            gisaidpath =  gisaidid  + ".fasta"

            #print(gisaidpath)
            if os.path.isfile(os.path.join('data/input/gisaidseqs',gisaidpath)):
                continue
            else:
                #print('No sequence for gisaid id: {} (run: {} bc: {})'.format(gisaidpath,run,barcode))
                del gisaidMapping[run][barcode]
        else:
            #print('No mapping for run: {} bc: {}'.format(run,barcode))
            pass

illuminaMapping = {}
for run in runs:
    illuminaMapping[run] = set()
    for barcode in barcodes[run]:
        if os.path.isdir(os.path.join('data/input/illumina',(run+'_'+barcode))):
            illuminaMapping[run].add(barcode)

### REMOVE LATER END


# This function assembles all the required output files that serve as a target for the snakemake pipeline
def getInput(wildcards):

    inputList = []

    for run in runs:
        if config['VarAnnotSnpEff']:
            inputList += expand('data/auxiliary/pangenome_vc/'+run+'/{barcode}/filter.annoted.vcf', barcode=barcodes[run])

        #VCF Based
        inputList += ['data/output/evaluation/vcfbased/comparison_pancov.tsv']
        inputList += ['data/output/evaluation/vcfbased/comparison_medaka.tsv']
        inputList += ['data/output/evaluation/vcfbased/comparison_nanopolish.tsv']

        for barcode in barcodes[run]:
            if run in illuminaMapping and barcode in illuminaMapping[run]:
                pass
            else:
                #print('skipping run {} barcode {} for illumina comparison as we have no seq yet ...'.format(run,barcode))
                continue

            inputList += [
                'data/auxiliary/illuminaVarCalls/'+run+'_'+barcode+'/ivar.vcf.tsv'
            ]

    return inputList

#Main rule that aggregates all the targets and is used when they are not specified, see function above for output files that are being created
rule all:
    input:
        getInput

include: 'rules/shared.snk'
include: 'rules/errorCorrection.snk'
include: 'rules/readAnalysis.snk'
include: 'rules/variantAnalysis.snk'
include: 'rules/consensus.snk'
include: 'rules/pangenome.snk'
include: 'rules/discovery.snk'
include: 'rules/pangenome_variant_call.snk'
### REMOVE LATER (THE INCLUDES BELOW ARE ONLY FOR EVALUATION)
include: 'rules/pangenome_eval.snk'
include: 'rules/illumina.snk'