Adjust MuTect output to match other callers.

    - Rename allelic fraction field in mutect output from FA to FREQ to standarize with other tools
    - Remove extra 'none' samples introduced when calling tumor-only samples
326473
def _fix_mutect_output(orig_file, config, out_file, is_paired):
    
    out_file_noc = out_file.replace(".vcf.gz", ".vcf")
    none_index = -1
    with file_transaction(config, out_file_noc) as tx_out_file:
        with open_gzipsafe(orig_file) as in_handle:
            with open(tx_out_file, 'w') as out_handle:
                for line in in_handle:
                    if not is_paired and line.startswith("#CHROM"):
                        parts = line.rstrip().split("\t")
                        none_index = parts.index("none")
                        del parts[none_index]
                        line = "\t".join(parts) + "\n"
                    elif line.startswith("##FORMAT=<ID=FA"):
                        line = line.replace("=FA", "=FREQ")
                    elif not line.startswith("#"):
                        if none_index > 0:
                            parts = line.rstrip().split("\t")
                            del parts[none_index]
                            line = "\t".join(parts) + "\n"
                        line = line.replace("FA", "FREQ")
                    out_handle.write(line)
    return bgzip_and_index(out_file_noc, config)
1036069
def reporter(self):
        
        logging.info('Preparing reports')
        # Populate self.plusdict in order to reuse parsing code from an assembly-based method
        for sample in self.runmetadata.samples:
            self.plusdict[sample.name] = dict()
            self.matchdict[sample.name] = dict()
            if sample.general.bestassemblyfile != 'NA':
                for gene in sample[self.analysistype].allelenames:
                    self.plusdict[sample.name][gene] = dict()
                    for allele, percentidentity in sample[self.analysistype].results.items():

                        if gene in allele:
                            # Split the allele number from the gene name using the appropriate delimiter
                            if '_' in allele:
                                splitter = '_'
                            elif '-' in allele:
                                splitter = '-'
                            else:
                                splitter = ''
                            self.matchdict[sample.name].update({gene: allele.split(splitter)[-1]})
                            # Create the plusdict dictionary as in the assembly-based (r)MLST method. Allows all the
                            # parsing and sequence typing code to be reused.
                            try:
                                self.plusdict[sample.name][gene][allele.split(splitter)[-1]][percentidentity] \
                                    = sample[self.analysistype].avgdepth[allele]
                            except KeyError:
                                self.plusdict[sample.name][gene][allele.split(splitter)[-1]] = dict()
                                self.plusdict[sample.name][gene][allele.split(splitter)[-1]][percentidentity] \
                                    = sample[self.analysistype].avgdepth[allele]
                    if gene not in self.matchdict[sample.name]:
                        self.matchdict[sample.name].update({gene: 'N'})
        self.profiler()
        self.sequencetyper()
        self.mlstreporter()
312260
def model_name(allele, num):
        
        random_string = hashlib.sha1(
            str(time.time()).encode()).hexdigest()[:16]
        return "%s-%d-%s" % (allele.upper(), num, random_string)
706806
def _clean_name(fname, data):
    
    for to_remove in dd.get_batches(data) + [dd.get_sample_name(data), data["sv"]["variantcaller"]]:
        for ext in ("-", "_"):
            if fname.startswith("%s%s" % (to_remove, ext)):
                fname = fname[len(to_remove) + len(ext):]
        if fname.startswith(to_remove):
            fname = fname[len(to_remove):]
    return fname
749762
def _run_scalpel_caller(align_bams, items, ref_file, assoc_files,
                          region=None, out_file=None):
    
    config = items[0]["config"]
    if out_file is None:
        out_file = "%s-variants.vcf.gz" % os.path.splitext(align_bams[0])[0]
    if not utils.file_exists(out_file):
        with file_transaction(config, out_file) as tx_out_file:
            if len(align_bams) > 1:
                message = ("Scalpel does not currently support batch calling!")
                raise ValueError(message)
            input_bams = " ".join("%s" % x for x in align_bams)
            tmp_path = "%s-scalpel-work" % utils.splitext_plus(out_file)[0]
            tx_tmp_path = "%s-scalpel-work" % utils.splitext_plus(tx_out_file)[0]
            if os.path.exists(tmp_path):
                utils.remove_safe(tmp_path)
            opts = " ".join(_scalpel_options_from_config(items, config, out_file, region, tmp_path))
            opts += " --dir %s" % tx_tmp_path
            min_cov = "3"  # minimum coverage
            opts += " --mincov %s" % min_cov
            perl_exports = utils.get_perl_exports(os.path.dirname(tx_out_file))
            cmd = ("{perl_exports} && "
                   "scalpel-discovery --single {opts} --ref {ref_file} --bam {input_bams} ")
            do.run(cmd.format(**locals()), "Genotyping with Scalpel", {})
            shutil.move(tx_tmp_path, tmp_path)
            # parse produced variant file further
            scalpel_tmp_file = bgzip_and_index(os.path.join(tmp_path, "variants.indel.vcf"), config)
            compress_cmd = "| bgzip -c" if out_file.endswith("gz") else ""
            bcftools_cmd_chi2 = get_scalpel_bcftools_filter_expression("chi2", config)
            sample_name_str = items[0]["name"][1]
            fix_ambig = vcfutils.fix_ambiguous_cl()
            add_contig = vcfutils.add_contig_to_header_cl(dd.get_ref_file(items[0]), tx_out_file)
            cl2 = ("{bcftools_cmd_chi2} {scalpel_tmp_file} | "
                   r"sed 's/FORMAT\tsample\(_name\)\{{0,1\}}/FORMAT\t{sample_name_str}/g' "
                   "| {fix_ambig} | vcfallelicprimitives -t DECOMPOSED --keep-geno | vcffixup - | vcfstreamsort "
                   "| {add_contig} {compress_cmd} > {tx_out_file}")
            do.run(cl2.format(**locals()), "Finalising Scalpel variants", {})
    return out_file
1365005
def run_rnaseq_variant_calling(data):
    
    variantcaller = dd.get_variantcaller(data)
    if isinstance(variantcaller, list) and len(variantcaller) > 1:
        logger.error("Only one variantcaller can be run for RNA-seq at "
                     "this time. Post an issue here "
                     "(https://github.com/bcbio/bcbio-nextgen/issues) "
                     "if this is something you need to do.")
        sys.exit(1)

    if variantcaller:
        if "gatk-haplotype" in variantcaller:
            data = variation.rnaseq_gatk_variant_calling(data)
        if vardict.get_vardict_command(data):
            data = variation.rnaseq_vardict_variant_calling(data)
        vrn_file = dd.get_vrn_file(data)
    return [[data]]
441367
def make_separate_indels_and_one_alt_with_all_snps_no_combinations(self, ref_seq):
        ''''''
        final_start_position = min([x.POS for x in self.vcf_records])
        final_end_position = max([x.ref_end_pos() for x in self.vcf_records])
        snps = []
        new_vcf_records = []

        for record in self.vcf_records:
            if record.is_snp():
                snps.append(copy.copy(record))
            else:
                new_record = copy.copy(record)
                new_record.add_flanking_seqs(ref_seq, final_start_position, final_end_position)
                new_vcf_records.append(new_record)

        if len(snps):
            new_record = copy.copy(snps[0])
            for snp in snps[1:]:
                merged = new_record.merge(snp, ref_seq)
                if merged is not None:
                    new_record = merged
            new_record.add_flanking_seqs(ref_seq, final_start_position, final_end_position)
            new_vcf_records.append(new_record)

        alts = ','.join(sorted(list(set([x.ALT[0] for x in new_vcf_records]))))
        new_record = vcf_record.VcfRecord('\t'.join([self.vcf_records[0].CHROM, str(final_start_position + 1), '.', new_vcf_records[0].REF, alts, '.', 'PASS', '.']))

        return new_record
259234
def run(vrn_info, calls_by_name, somatic_info, do_plots=True, handle_failures=True):
    
    if "seq2c" in calls_by_name:
        cnv_info = calls_by_name["seq2c"]
    elif "cnvkit" in calls_by_name:
        cnv_info = calls_by_name["cnvkit"]
    else:
        raise ValueError("BubbleTree only currently support CNVkit and Seq2c: %s" % ", ".join(calls_by_name.keys()))
    work_dir = _cur_workdir(somatic_info.tumor_data)
    class OutWriter:
        def __init__(self, out_handle):
            self.writer = csv.writer(out_handle)
        def write_header(self):
            self.writer.writerow(["chrom", "start", "end", "freq"])
        def write_row(self, rec, stats):
            self.writer.writerow([_to_ucsc_style(rec.chrom), rec.start, rec.stop, stats["tumor"]["freq"]])
    vcf_csv = prep_vrn_file(vrn_info["vrn_file"], vrn_info["variantcaller"],
                            work_dir, somatic_info, OutWriter, cnv_info["cns"])
    cnv_csv = _prep_cnv_file(cnv_info["cns"], cnv_info["variantcaller"], work_dir,
                             somatic_info.tumor_data)
    wide_lrr = cnv_info["variantcaller"] == "cnvkit" and somatic_info.normal_bam is None
    return _run_bubbletree(vcf_csv, cnv_csv, somatic_info.tumor_data, wide_lrr, do_plots,
                           handle_failures)
301357
def _dup_samples_by_variantcaller(samples, require_bam=True):
    
    samples = [utils.to_single_data(x) for x in samples]
    samples = germline.split_somatic(samples)
    to_process = []
    extras = []
    for data in samples:
        added = False
        for i, add in enumerate(handle_multiple_callers(data, "variantcaller", require_bam=require_bam)):
            added = True
            add = dd.set_variantcaller_order(add, i)
            to_process.append([add])
        if not added:
            data = _handle_precalled(data)
            data = dd.set_variantcaller_order(data, 0)
            extras.append([data])
    return to_process, extras
321060
def _check_freqs(parts, indexes):
    
    thresh_ratio = 2.7
    try:  # FreeBayes
        ao_index = parts[8].split(":").index("AO")
        ro_index = parts[8].split(":").index("RO")
    except ValueError:
        ao_index, ro_index = None, None
    try:  # VarDict
        af_index = parts[8].split(":").index("AF")
    except ValueError:
        af_index = None
    if af_index is None and ao_index is None:
        # okay to skip if a gVCF record
        if parts[4].find("<*>") == -1:
            raise NotImplementedError("Unexpected format annotations: %s" % parts[8])
    def _calc_freq(item):
        try:
            if ao_index is not None and ro_index is not None:
                ao = sum([int(x) for x in item.split(":")[ao_index].split(",")])
                ro = int(item.split(":")[ro_index])
                freq = ao / float(ao + ro)
            elif af_index is not None:
                freq = float(item.split(":")[af_index])
            else:
                freq = 0.0
        except (IndexError, ValueError, ZeroDivisionError):
            freq = 0.0
        return freq
    tumor_freq, normal_freq = _calc_freq(parts[indexes["tumor"]]), _calc_freq(parts[indexes["normal"]])
    return normal_freq <= 0.001 or normal_freq <= tumor_freq / thresh_ratio
228589
def _filter_seqs(fn):
    """"""
    out_file = op.splitext(fn)[0] + "_unique.fa"
    idx = 0
    if not file_exists(out_file):
        with open(out_file, 'w') as out_handle:
            with open(fn) as in_handle:
                for line in in_handle:
                    if line.startswith("@") or line.startswith(">"):
                        fixed_name = _make_unique(line.strip(), idx)
                        seq = in_handle.next().strip()
                        counts = _get_freq(fixed_name)
                        if len(seq) < 26 and (counts > 1 or counts == 0):
                            idx += 1
                            print(fixed_name, file=out_handle, end="\n")
                            print(seq, file=out_handle, end="\n")
                        if line.startswith("@"):
                            in_handle.next()
                            in_handle.next()
    return out_file
1192181
def rnaseq_variant_calling(samples, run_parallel):
    
    samples = run_parallel("run_rnaseq_variant_calling", samples)
    variantcaller = dd.get_variantcaller(to_single_data(samples[0]))
    if variantcaller and ("gatk-haplotype" in variantcaller):
        out = []
        for d in joint.square_off(samples, run_parallel):
            out.extend([[to_single_data(xs)] for xs in multi.split_variants_by_sample(to_single_data(d))])
        samples = out
    if variantcaller:
        samples = run_parallel("run_rnaseq_ann_filter", samples)
    if variantcaller and ("gatk-haplotype" in variantcaller):
        out = []
        for data in (to_single_data(xs) for xs in samples):
            if "variants" not in data:
                data["variants"] = []
            data["variants"].append({"variantcaller": "gatk-haplotype", "vcf": data["vrn_file_orig"],
                                     "population": {"vcf": data["vrn_file"]}})
            data["vrn_file"] = data.pop("vrn_file_orig")
            out.append([data])
        samples = out
    return samples
5364
def _organize_variants(samples, batch_id):
    
    caller_names = [x["variantcaller"] for x in samples[0]["variants"]]
    calls = collections.defaultdict(list)
    for data in samples:
        for vrn in data["variants"]:
            calls[vrn["variantcaller"]].append(vrn["vrn_file"])
    data = samples[0]
    vrn_files = []
    for caller in caller_names:
        fnames = calls[caller]
        if len(fnames) == 1:
            vrn_files.append(fnames[0])
        else:
            vrn_files.append(population.get_multisample_vcf(fnames, batch_id, caller, data))
    return caller_names, vrn_files
673965
def _covert_to_diploid(in_file, data):
    
    sample = dd.get_sample_name(data)
    out_file = "%s-diploid.vcf" % utils.splitext_plus(in_file)[0]
    in_vcf = pysam.VariantFile(in_file)
    out_vcf = pysam.VariantFile(out_file, 'w', header=in_vcf.header)
    for record in in_vcf:
        gt = list(record.samples[sample]['GT'])
        if 'SOMATIC' in record.info:
            for allele in set(gt):
                if allele != gt[0]:
                    record.samples[sample]['GT'] = gt[0], allele
                    out_vcf.write(record)
        else:
            if len(gt) == 1:
                record.samples[sample]['GT'] = gt
            else:
                record.samples[sample]['GT'] = gt[0], gt[1]
            out_vcf.write(record)
    in_vcf.close()
    out_vcf.close()
    return vcfutils.bgzip_and_index(out_file, data["config"])
773930
def run_rnaseq_ann_filter(data):
    
    data = to_single_data(data)
    if dd.get_vrn_file(data):
        eff_file = effects.add_to_vcf(dd.get_vrn_file(data), data)[0]
        if eff_file:
            data = dd.set_vrn_file(data, eff_file)
        ann_file = population.run_vcfanno(dd.get_vrn_file(data), data)
        if ann_file:
            data = dd.set_vrn_file(data, ann_file)
    variantcaller = dd.get_variantcaller(data)
    if variantcaller and ("gatk-haplotype" in variantcaller):
        filter_file = variation.gatk_filter_rnaseq(dd.get_vrn_file(data), data)
        data = dd.set_vrn_file(data, filter_file)
    # remove variants close to splice junctions
    vrn_file = dd.get_vrn_file(data)
    vrn_file = variation.filter_junction_variants(vrn_file, data)
    data = dd.set_vrn_file(data, vrn_file)
    return [[data]]
1274910
def _filter_vcf(out_file):
    
    in_file = out_file.replace(".vcf", "-ori.vcf")
    FILTER_line = ('##FILTER=<ID=SBIAS,Description="Due to bias">\n'
                   '##FILTER=<ID=5BP,Description="Due to 5BP">\n'
                   '##FILTER=<ID=REJECT,Description="Not somatic due to qSNP filters">\n')
    SOMATIC_line = '##INFO=<ID=SOMATIC,Number=0,Type=Flag,Description="somatic event">\n'
    if not utils.file_exists(in_file):
        shutil.move(out_file, in_file)
    with file_transaction(out_file) as tx_out_file:
        with open(in_file) as in_handle, open(tx_out_file, "w") as out_handle:
            for line in in_handle:
                if line.startswith("##normalSample="):
                    normal_name = line.strip().split("=")[1]
                if line.startswith("##patient_id="):
                    tumor_name = line.strip().split("=")[1]
                if line.startswith("#CHROM"):
                    line = line.replace("Normal", normal_name)
                    line = line.replace("Tumour", tumor_name)
                if line.startswith("##INFO=<ID=FS"):
                    line = line.replace("ID=FS", "ID=RNT")
                if line.find("FS=") > -1:
                    line = line.replace("FS=", "RNT=")
                if "5BP" in line:
                    line = sub("5BP[0-9]+", "5BP", line)
                if line.find("PASS") == -1:
                    line = _set_reject(line)
                if line.find("PASS") > - 1 and line.find("SOMATIC") == -1:
                    line = _set_reject(line)
                if not _has_ambiguous_ref_allele(line):
                    out_handle.write(line)
                if line.startswith("##FILTER") and FILTER_line:
                    out_handle.write("%s" % FILTER_line)
                    FILTER_line = ""
                if line.startswith("##INFO") and SOMATIC_line:
                    out_handle.write("%s" % SOMATIC_line)
                    SOMATIC_line = ""
    return out_file
356599
def demultiplex_cells(fastq, out_dir, readnumber, prefix, cb_histogram, 
                      cb_cutoff):
    ''' 
    '''
    annotations = detect_fastq_annotations(fastq)
    re_string = construct_transformed_regex(annotations)
    parser_re = re.compile(re_string)
    readstring = "" if not readnumber else "_R{}".format(readnumber)
    filestring = "{prefix}{sample}{readstring}.fq"
    cb_set = set()
    if cb_histogram:
        cb_set = get_cb_depth_set(cb_histogram, cb_cutoff)
    sample_set = set()
    batch = collections.defaultdict(list)
    parsed = 0
    safe_makedir(out_dir)
    for read in read_fastq(fastq):
        parsed += 1
        match = parser_re.search(read).groupdict()
        sample = match['CB']
        if cb_set and sample not in cb_set:
            continue
        sample_set.add(sample)
        batch[sample].append(read)
        # write in batches to avoid opening up file handles repeatedly
        if not parsed % 10000000:
            for sample, reads in batch.items():
                out_file = os.path.join(out_dir, filestring.format(**locals()))
                with open(out_file, "a") as out_handle:
                    for read in reads:
                        out_handle.write(read)
            batch = collections.defaultdict(list)

    for sample, reads in batch.items():
        out_file = os.path.join(out_dir, filestring.format(**locals()))
        with open(out_file, "a") as out_handle:
            for read in reads:
                out_handle.write(read)