A patch fixing error that occurs when using 0:1 read mask with batch CdrBlast
A patch fixing error that occurs when using 0:1 read mask with batch CdrBlast

package com.milaboratory.migec

/**
 Copyright 2014 Mikhail Shugay (mikhail.shugay@gmail.com)

 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at

 http://www.apache.org/licenses/LICENSE-2.0

 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
 */

def scriptName = getClass().canonicalName

def cli = new CliBuilder(usage: "$scriptName [options] [checkout_dir/ or -] [assemble_dir/ or -] output_dir/\n" +
        "Either --sample-metadata or -R argument is required.")
cli.p(args: 1, "Number of threads to use")
cli._(longOpt: "sample-metadata", args: 1, argName: "file",
        "An optional tab-delimited file indicating which samples to process and containing 6 columns:\n" +
                "sample_id|species|gene|file_types|mask|quality_threshold_pair\n" +
                "Samples not included will be omitted by CdrBlast. " +
                "File types, mask and quality threshold columns could be either omitted or set to " - "")
cli._(args: 1, longOpt: "blast-path",
        "Path to blast executable.")
cli._(longOpt: "no-sort",
        "Do not sort output files, which will speed up processing. " +
                "Could be used for full pipeline as FilterCdrBlastResults will provide final clonotype table in sorted format.")
cli._(longOpt: "all-segments",
        "Use full V/D/J segment library (including pseudogens, etc).")
cli._(longOpt: "print-library",
        "Prints out allowed species-gene pairs. " +
                "To account non-functional segment data use together with --all-segments")
cli._(longOpt: "default-mask", args: 1, argName: "R1=0/1:R2=0/1",
        "Mask, default for all samples, will be ignored for unpaired and overlapped reads, " +
                "allowed values are: 1:0 (R1 assembled), 0:1 (R2 assembled) and 1:1 (both reads assembled, default)")
cli.R(longOpt: "default-gene", args: 1, argName: "gene",
        "Gene, default for all samples, allowed values: 'TRA', 'TRB', 'TRG', 'TRD', 'IGL', 'IGK' or 'IGH'. " +
                "Use --print-library for the list of allowed species-gene combinations.")
cli._(longOpt: "default-species", args: 1, argName: "species",
        "Species, default for all samples, allowed values: 'HomoSapiens'[default], 'MusMusculus', ... " +
                "Use --print-library for the list of allowed species-gene combinations.")
cli._(longOpt: "default-file-types", args: 1, argName: "type1,..",
        "Accepted file types, default for all samples, " +
                "comma-separated list of file types to be processed for a given sample, allowed values: ${Util.FILE_TYPES.join(", ")}")
cli._(longOpt: "default-quality-threshold", args: 1, argName: "Phred,CQS",
        "Quality threshold pair, default for all samples, comma-separated pair " +
                "of quality threshold values for Phred and CQS quality thresholds respectively [default = 25, 30]")

def opt = cli.parse(args)

if (opt == null) {
    println "[ERROR] Too few arguments provided"
    cli.usage()
    System.exit(-1)
}

// SEGMENTS
boolean includeNonFuncitonal = opt.'all-segments'
if (opt.'print-library') {
    println "CDR3 extraction is possible for the following data (segments include non-functional = $includeNonFuncitonal):"
    Util.listAvailableSegments(includeNonFuncitonal)
    System.exit(0)
}

if (opt.arguments().size() < 3) {
    println "[ERROR] Too few arguments provided"
    cli.usage()
    System.exit(-1)
}

def blastPath = opt.'blast-path' ?: null

// INPUT AND OUTPUT FILES
String sampleInfoFileName = opt.'sample-metadata' ?: null, checkoutDir = opt.arguments()[0],
       assembleDir = opt.arguments()[1], outputPath = opt.arguments()[2]

if (!sampleInfoFileName && !opt.R) {
    println "[ERROR] Either --sample-metadata or --default-gene options should be provided"
    System.exit(-1)
}

boolean processRaw = checkoutDir != "-", processAssembled = assembleDir != "-"
if (!processRaw && !processAssembled) {
    println "[ERROR] At least one of assembled or raw (checkout) data directories should be specified"
    System.exit(-1)
}
boolean processBoth = processRaw && processAssembled

// LOAD SAMPLE ENTRIES
Map rawSamplesMap = null
if (processRaw) {
    def checkoutSampleListFile = new File(checkoutDir + "/checkout.filelist.txt")
    if (!checkoutSampleListFile.exists()) {
        println "[ERROR] Sample list (checkout.filelist.txt) absent in Checkout folder $checkoutDir"
        System.exit(-1)
    }

    rawSamplesMap = checkoutSampleListFile.readLines().findAll {
        !it.startsWith("#")   // no headers
    }.collectEntries {
        def splitLine = it.split('\t')
        def sampleName = splitLine[0], sampleType = splitLine[1]
        [("$sampleName\t$sampleType".toString()): splitLine[2..3]]
    }
}

Map assembledSampleMap = null
if (processAssembled) {
    def assembleSampleListFile = new File(assembleDir + "/assemble.log.txt")
    if (!assembleSampleListFile.exists()) {
        println "[ERROR] Assemble output (assemble.log.txt) absent in Assemble folder $assembleDir"
        System.exit(-1)
    }

    assembledSampleMap = assembleSampleListFile.readLines().findAll {
        !it.startsWith("#")   // no headers
    }.collectEntries {
        def splitLine = it.split('\t')
        def sampleName = splitLine[0], sampleType = splitLine[1], sampleKey = "$sampleName\t$sampleType".toString()

        [(sampleKey): [splitLine[2..3], splitLine[4..5]]]
    }
}

// GENERATE SAMPLE METADATA FOR CDRBLAST
def sampleInfoMap = new HashMap()
List<String> sampleInfoLines = []

if (sampleInfoFileName) {
    // Read sample metadata
    def chainsFile = new File(sampleInfoFileName)

    if (!chainsFile.exists()) {
        println "[ERROR] Sample metadata file not found"
        System.exit(-1)
    }

    sampleInfoLines = chainsFile.readLines()
} else {
    // Or generate it from default parameters
    def defaultMask = opt.'default-mask' ?: "1:1", defaultSpecies = opt.'default-species' ?: "homosapiens",
        defaultFileTypes = opt.'default-file-types' ?: Util.FILE_TYPES.join(","),
        defaultQualityThreshold = opt.'default-quality-threshold' ?: "25,30", defaultChain = opt.R

    def defaultParameters = [defaultSpecies, defaultChain,
                             defaultFileTypes, defaultMask, defaultQualityThreshold].join("\t")

    (processRaw ? rawSamplesMap : assembledSampleMap).keySet().each {
        def sampleId = it.split("\t")[0]
        sampleInfoLines.add("$sampleId\t" + defaultParameters)
    }
}

// CHECK SAMPLE METADATA
sampleInfoLines.findAll { !it.startsWith("#") }.each { line ->
    def splitLine = line.split("\t")
    def sampleId = splitLine[0],
        species = splitLine[1], chain = splitLine[2],
        fileTypes = splitLine.length > 3 ? (splitLine[3] == '-' ? Util.FILE_TYPES.join(",") : splitLine[3]) : Util.FILE_TYPES.join(","),
        mask = (splitLine.length > 4 ? (splitLine[4] == '-' ? "1:1" : splitLine[4]) : "1:1"),
        qualityThreshold = splitLine.length > 5 ? (splitLine[5] == '-' ? null : splitLine[5]) : null

    qualityThreshold = qualityThreshold ? qualityThreshold.split(",") : null
    // todo: estimate q threshold by HistQ

    if (!Util.isAvailable(species, chain, includeNonFuncitonal)) {
        println "[ERROR] Sorry, no analysis could be performed for $species gene $chain " +
                "(include non-functional = $includeNonFuncitonal). " +
                "Possible variants are:\n"
        Util.listAvailableSegments(includeNonFuncitonal)
        System.exit(-1)
    }
    if (!Util.MASKS.any { mask == it }) {
        println "[ERROR] Bad mask $mask on line $line. Supported masks are ${Util.MASKS.join(", ")}"
        System.exit(-1)
    }
    mask.split(",").collect { it == "1" }
    fileTypes = fileTypes.split(",")
    fileTypes.each { fileType ->
        if (!Util.FILE_TYPES.any { fileType == it }) {
            println "[ERROR] Bad file type $fileType in line $line. Supported file types are ${Util.FILE_TYPES.join(", ")}"
            System.exit(-1)
        }
    }

    // check if sample exists
    boolean exists = false
    fileTypes.each {
        def sampleKey = "$sampleId\t$it".toString()
        if (processRaw && rawSamplesMap.containsKey(sampleKey))
            exists = true
        if (processAssembled && assembledSampleMap.containsKey(sampleKey))
            exists = true
    }
    if (!exists) {
        println "[ERROR] Sample $sampleId with file types ${fileTypes.join(",")} does not exist in checkout/assembled sample lists"
        System.exit(-1)
    }

    sampleInfoMap.put(sampleId, [species, chain, fileTypes, mask, qualityThreshold])
}

// START LOGGING
def logFile = new File(outputPath + "/cdrblast.log.txt")
if (logFile.exists())
    logFile.delete()
else
    logFile.absoluteFile.parentFile.mkdirs()

logFile.withPrintWriter { pw ->
    pw.println(Util.CDRBLAST_LOG_HEADER)

    // RUN CDRBLAST
    def baseArgs = []
    if (opt.'no-sort')
        baseArgs = [baseArgs, ["--no-sort"]]
    if (opt.p)
        baseArgs = [baseArgs, ["-p", opt.p]]
    if (blastPath)
        baseArgs = [baseArgs, ["--blast-path", blastPath]]
    if (includeNonFuncitonal)
        baseArgs = [baseArgs, ["--all-segments"]]
    String logLine
    if (processAssembled) {
        if (processBoth)
            println "[${new Date()} $scriptName] Running CdrBlast for raw and assembled data.."
        else
            println "[${new Date()} $scriptName] Running CdrBlast for assembled data.."

        sampleInfoMap.each { chainData ->
            def sampleId = chainData.key, species = chainData.value[0], chain = chainData.value[1],
                fileTypes = chainData.value[2], mask = chainData.value[3], qualityThreshold = chainData.value[4]

            def assemblyFiles = [], rawFiles = []

            fileTypes.each { fileType ->
                def sampleKey = "$sampleId\t$fileType".toString()

                if (assembledSampleMap.containsKey(sampleKey)) {
                    def correspondingFiles = assembledSampleMap[sampleKey]

                    // apply mask
                    if (fileType == 'paired') {
                        boolean zeroMask = true
                        if (mask[0]) {
                            rawFiles.add(correspondingFiles[0][0])
                            assemblyFiles.add(correspondingFiles[1][0])
                            if (correspondingFiles[1][0] != '-')
                                zeroMask = false
                        }
                        if (mask[1]) {
                            rawFiles.add(correspondingFiles[0][1])
                            assemblyFiles.add(correspondingFiles[1][1])
                            if (correspondingFiles[1][1] != '-')
                                zeroMask = false
                        }
                        if (zeroMask) {
                            println "[ERROR] Sample $sampleId with mask ${mask.collect { it ? 1 : 0 }} " +
                                    "does not have any corresponding assembled files, " +
                                    "of ${correspondingFiles[1]}"
                            System.exit(-1)
                        }
                    } else {
                        rawFiles.add(correspondingFiles[0][0])
                        rawFiles.add(correspondingFiles[0][1])
                        assemblyFiles.add(correspondingFiles[1][0])
                        assemblyFiles.add(correspondingFiles[1][1])
                    }
                }
            }

            // cleanup from masked files
            assemblyFiles.removeAll { it == '-' }
            rawFiles.removeAll { it == '-' }

            def baseArgs1 = [baseArgs, ["-R", chain], ["-S", species]]

            logLine = Util.run(new CdrBlast(), [baseArgs1, ["-a"], ["-q", qualityThreshold[1]],
                                                assemblyFiles, "$outputPath/${sampleId}.asm.cdrblast.txt"].flatten().join(" "))
            pw.println("$sampleId\t" + logLine)

            if (processBoth) {
                logLine = Util.run(new CdrBlast(), [baseArgs1, ["-q", qualityThreshold[0]],
                                                    rawFiles, "$outputPath/${sampleId}.raw.cdrblast.txt"].flatten().join(" "))
                pw.println("$sampleId\t" + logLine)
            }
        }
    } else {
        println "[${new Date()} $scriptName] Running CdrBlast for raw data.."
        sampleInfoMap.each { chainData ->
            def sampleId = chainData.key, species = chainData.value[0], chain = chainData.value[1],
                fileTypes = chainData.value[2], mask = chainData.value[3], qualityThreshold = chainData.value[4]

            def rawFiles = []

            fileTypes.each { fileType ->
                def sampleKey = "$sampleId\t$fileType".toString()

                if (rawSamplesMap.containsKey(sampleKey)) {
                    def correspondingFiles = rawSamplesMap[sampleKey]

                    // apply mask
                    if (fileType == 'paired') {
                        if (mask[0])
                            rawFiles.add(correspondingFiles[0])

                        if (mask[1])
                            rawFiles.add(correspondingFiles[1])
                    } else
                        rawFiles.add(correspondingFiles[0])
                }
            }

            logLine = Util.run(new CdrBlast(), [baseArgs, ["-R", chain], ["-S", species], ["-q", qualityThreshold[0]],
                                                rawFiles, "$outputPath/${sampleId}.raw.cdrblast.txt"].flatten().join(" "))
            pw.println("$sampleId\t" + logLine)
        }
    }
}
package com.milaboratory.migec

/**
 Copyright 2014 Mikhail Shugay (mikhail.shugay@gmail.com)

 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at

 http://www.apache.org/licenses/LICENSE-2.0

 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
 */

def scriptName = getClass().canonicalName

def cli = new CliBuilder(usage: "$scriptName [options] [checkout_dir/ or -] [assemble_dir/ or -] output_dir/\n" +
        "Either --sample-metadata or -R argument is required.")
cli.p(args: 1, "Number of threads to use")
cli._(longOpt: "sample-metadata", args: 1, argName: "file",
        "An optional tab-delimited file indicating which samples to process and containing 6 columns:\n" +
                "sample_id|species|gene|file_types|mask|quality_threshold_pair\n" +
                "Samples not included will be omitted by CdrBlast. " +
                "File types, mask and quality threshold columns could be either omitted or set to " - "")
cli._(args: 1, longOpt: "blast-path",
        "Path to blast executable.")
cli._(longOpt: "no-sort",
        "Do not sort output files, which will speed up processing. " +
                "Could be used for full pipeline as FilterCdrBlastResults will provide final clonotype table in sorted format.")
cli._(longOpt: "all-segments",
        "Use full V/D/J segment library (including pseudogens, etc).")
cli._(longOpt: "print-library",
        "Prints out allowed species-gene pairs. " +
                "To account non-functional segment data use together with --all-segments")
cli._(longOpt: "default-mask", args: 1, argName: "R1=0/1:R2=0/1",
        "Mask, default for all samples, will be ignored for unpaired and overlapped reads, " +
                "allowed values are: 1:0 (R1 assembled), 0:1 (R2 assembled) and 1:1 (both reads assembled, default)")
cli.R(longOpt: "default-gene", args: 1, argName: "gene",
        "Gene, default for all samples, allowed values: 'TRA', 'TRB', 'TRG', 'TRD', 'IGL', 'IGK' or 'IGH'. " +
                "Use --print-library for the list of allowed species-gene combinations.")
cli._(longOpt: "default-species", args: 1, argName: "species",
        "Species, default for all samples, allowed values: 'HomoSapiens'[default], 'MusMusculus', ... " +
                "Use --print-library for the list of allowed species-gene combinations.")
cli._(longOpt: "default-file-types", args: 1, argName: "type1,..",
        "Accepted file types, default for all samples, " +
                "comma-separated list of file types to be processed for a given sample, allowed values: ${Util.FILE_TYPES.join(", ")}")
cli._(longOpt: "default-quality-threshold", args: 1, argName: "Phred,CQS",
        "Quality threshold pair, default for all samples, comma-separated pair " +
                "of quality threshold values for Phred and CQS quality thresholds respectively [default = 25, 30]")

def opt = cli.parse(args)

if (opt == null) {
    println "[ERROR] Too few arguments provided"
    cli.usage()
    System.exit(-1)
}

// SEGMENTS
boolean includeNonFuncitonal = opt.'all-segments'
if (opt.'print-library') {
    println "CDR3 extraction is possible for the following data (segments include non-functional = $includeNonFuncitonal):"
    Util.listAvailableSegments(includeNonFuncitonal)
    System.exit(0)
}

if (opt.arguments().size() < 3) {
    println "[ERROR] Too few arguments provided"
    cli.usage()
    System.exit(-1)
}

def blastPath = opt.'blast-path' ?: null

// INPUT AND OUTPUT FILES
String sampleInfoFileName = opt.'sample-metadata' ?: null, checkoutDir = opt.arguments()[0],
       assembleDir = opt.arguments()[1], outputPath = opt.arguments()[2]

if (!sampleInfoFileName && !opt.R) {
    println "[ERROR] Either --sample-metadata or --default-gene options should be provided"
    System.exit(-1)
}

boolean processRaw = checkoutDir != "-", processAssembled = assembleDir != "-"
if (!processRaw && !processAssembled) {
    println "[ERROR] At least one of assembled or raw (checkout) data directories should be specified"
    System.exit(-1)
}
boolean processBoth = processRaw && processAssembled

// LOAD SAMPLE ENTRIES
Map rawSamplesMap = null
if (processRaw) {
    def checkoutSampleListFile = new File(checkoutDir + "/checkout.filelist.txt")
    if (!checkoutSampleListFile.exists()) {
        println "[ERROR] Sample list (checkout.filelist.txt) absent in Checkout folder $checkoutDir"
        System.exit(-1)
    }

    rawSamplesMap = checkoutSampleListFile.readLines().findAll {
        !it.startsWith("#")   // no headers
    }.collectEntries {
        def splitLine = it.split('\t')
        def sampleName = splitLine[0], sampleType = splitLine[1]
        [("$sampleName\t$sampleType".toString()): splitLine[2..3]]
    }
}

Map assembledSampleMap = null
if (processAssembled) {
    def assembleSampleListFile = new File(assembleDir + "/assemble.log.txt")
    if (!assembleSampleListFile.exists()) {
        println "[ERROR] Assemble output (assemble.log.txt) absent in Assemble folder $assembleDir"
        System.exit(-1)
    }

    assembledSampleMap = assembleSampleListFile.readLines().findAll {
        !it.startsWith("#")   // no headers
    }.collectEntries {
        def splitLine = it.split('\t')
        def sampleName = splitLine[0], sampleType = splitLine[1], sampleKey = "$sampleName\t$sampleType".toString()

        [(sampleKey): [splitLine[2..3], splitLine[4..5]]]
    }
}

// GENERATE SAMPLE METADATA FOR CDRBLAST
def sampleInfoMap = new HashMap()
List<String> sampleInfoLines = []

if (sampleInfoFileName) {
    // Read sample metadata
    def chainsFile = new File(sampleInfoFileName)

    if (!chainsFile.exists()) {
        println "[ERROR] Sample metadata file not found"
        System.exit(-1)
    }

    sampleInfoLines = chainsFile.readLines()
} else {
    // Or generate it from default parameters
    def defaultMask = opt.'default-mask' ?: "1:1", defaultSpecies = opt.'default-species' ?: "homosapiens",
        defaultFileTypes = opt.'default-file-types' ?: Util.FILE_TYPES.join(","),
        defaultQualityThreshold = opt.'default-quality-threshold' ?: "25,30", defaultChain = opt.R

    def defaultParameters = [defaultSpecies, defaultChain,
                             defaultFileTypes, defaultMask, defaultQualityThreshold].join("\t")

    (processRaw ? rawSamplesMap : assembledSampleMap).keySet().each {
        def sampleId = it.split("\t")[0]
        sampleInfoLines.add("$sampleId\t" + defaultParameters)
    }
}

// CHECK SAMPLE METADATA
sampleInfoLines.findAll { !it.startsWith("#") }.each { line ->
    def splitLine = line.split("\t")
    def sampleId = splitLine[0],
        species = splitLine[1], chain = splitLine[2],
        fileTypes = splitLine.length > 3 ? (splitLine[3] == '-' ? Util.FILE_TYPES.join(",") : splitLine[3]) : Util.FILE_TYPES.join(","),
        mask = (splitLine.length > 4 ? (splitLine[4] == '-' ? "1:1" : splitLine[4]) : "1:1"),
        qualityThreshold = splitLine.length > 5 ? (splitLine[5] == '-' ? null : splitLine[5]) : null

    qualityThreshold = qualityThreshold ? qualityThreshold.split(",") : null
    // todo: estimate q threshold by HistQ

    if (!Util.isAvailable(species, chain, includeNonFuncitonal)) {
        println "[ERROR] Sorry, no analysis could be performed for $species gene $chain " +
                "(include non-functional = $includeNonFuncitonal). " +
                "Possible variants are:\n"
        Util.listAvailableSegments(includeNonFuncitonal)
        System.exit(-1)
    }
    if (!Util.MASKS.any { mask == it }) {
        println "[ERROR] Bad mask $mask on line $line. Supported masks are ${Util.MASKS.join(", ")}"
        System.exit(-1)
    }
    mask.split(",").collect { it == "1" }
    fileTypes = fileTypes.split(",")
    fileTypes.each { fileType ->
        if (!Util.FILE_TYPES.any { fileType == it }) {
            println "[ERROR] Bad file type $fileType in line $line. Supported file types are ${Util.FILE_TYPES.join(", ")}"
            System.exit(-1)
        }
    }

    // check if sample exists
    boolean exists = false
    fileTypes.each {
        def sampleKey = "$sampleId\t$it".toString()
        if (processRaw && rawSamplesMap.containsKey(sampleKey))
            exists = true
        if (processAssembled && assembledSampleMap.containsKey(sampleKey))
            exists = true
    }
    if (!exists) {
        println "[ERROR] Sample $sampleId with file types ${fileTypes.join(",")} does not exist in checkout/assembled sample lists"
        System.exit(-1)
    }

    sampleInfoMap.put(sampleId, [species, chain, fileTypes, mask, qualityThreshold])
}

// START LOGGING
def logFile = new File(outputPath + "/cdrblast.log.txt")
if (logFile.exists())
    logFile.delete()
else
    logFile.absoluteFile.parentFile.mkdirs()

logFile.withPrintWriter { pw ->
    pw.println(Util.CDRBLAST_LOG_HEADER)

    // RUN CDRBLAST
    def baseArgs = []
    if (opt.'no-sort')
        baseArgs = [baseArgs, ["--no-sort"]]
    if (opt.p)
        baseArgs = [baseArgs, ["-p", opt.p]]
    if (blastPath)
        baseArgs = [baseArgs, ["--blast-path", blastPath]]
    if (includeNonFuncitonal)
        baseArgs = [baseArgs, ["--all-segments"]]
    String logLine
    if (processAssembled) {
        if (processBoth)
            println "[${new Date()} $scriptName] Running CdrBlast for raw and assembled data.."
        else
            println "[${new Date()} $scriptName] Running CdrBlast for assembled data.."

        sampleInfoMap.each { chainData ->
            def sampleId = chainData.key, species = chainData.value[0], chain = chainData.value[1],
                fileTypes = chainData.value[2], mask = chainData.value[3], qualityThreshold = chainData.value[4]

            def assemblyFiles = [], rawFiles = []

            fileTypes.each { fileType ->
                def sampleKey = "$sampleId\t$fileType".toString()

                if (assembledSampleMap.containsKey(sampleKey)) {
                    def correspondingFiles = assembledSampleMap[sampleKey]

                    // apply mask
                    if (fileType == 'paired') {
                        boolean zeroMask = true
                        if (mask[0]) {
                            rawFiles.add(correspondingFiles[0][0])
                            assemblyFiles.add(correspondingFiles[1][0])
                            if (correspondingFiles[1][0] != '-')
                                zeroMask = false
                        }
                        if (mask[1]) {
                            rawFiles.add(correspondingFiles[0][1])
                            assemblyFiles.add(correspondingFiles[1][1])
                            if (correspondingFiles[1][1] != '-')
                                zeroMask = false
                        }
                        if (zeroMask) {
                            println "[ERROR] Sample $sampleId with mask ${mask.collect { it ? 1 : 0 }} " +
                                    "does not have any corresponding assembled files, " +
                                    "of ${correspondingFiles[1]}"
                            System.exit(-1)
                        }
                        assemblyFiles.removeAll { it == '-' } // cleanup from files masked on assembly stage
                    } else {
                        rawFiles.add(correspondingFiles[0][0])
                        assemblyFiles.add(correspondingFiles[1][0])
                    }
                }
            }

            def baseArgs1 = [baseArgs, ["-R", chain], ["-S", species]]

            logLine = Util.run(new CdrBlast(), [baseArgs1, ["-a"], ["-q", qualityThreshold[1]],
                                                assemblyFiles, "$outputPath/${sampleId}.asm.cdrblast.txt"].flatten().join(" "))
            pw.println("$sampleId\t" + logLine)

            if (processBoth) {
                logLine = Util.run(new CdrBlast(), [baseArgs1, ["-q", qualityThreshold[0]],
                                                    rawFiles, "$outputPath/${sampleId}.raw.cdrblast.txt"].flatten().join(" "))
                pw.println("$sampleId\t" + logLine)
            }
        }
    } else {
        println "[${new Date()} $scriptName] Running CdrBlast for raw data.."
        sampleInfoMap.each { chainData ->
            def sampleId = chainData.key, species = chainData.value[0], chain = chainData.value[1],
                fileTypes = chainData.value[2], mask = chainData.value[3], qualityThreshold = chainData.value[4]

            def rawFiles = []

            fileTypes.each { fileType ->
                def sampleKey = "$sampleId\t$fileType".toString()

                if (rawSamplesMap.containsKey(sampleKey)) {
                    def correspondingFiles = rawSamplesMap[sampleKey]

                    // apply mask
                    if (fileType == 'paired') {
                        if (mask[0])
                            rawFiles.add(correspondingFiles[0])

                        if (mask[1])
                            rawFiles.add(correspondingFiles[1])
                    } else
                        rawFiles.add(correspondingFiles[0])
                }
            }

            logLine = Util.run(new CdrBlast(), [baseArgs, ["-R", chain], ["-S", species], ["-q", qualityThreshold[0]],
                                                rawFiles, "$outputPath/${sampleId}.raw.cdrblast.txt"].flatten().join(" "))
            pw.println("$sampleId\t" + logLine)
        }
    }
}