Puzzle: manage DNA variant resources.
510334
def base(ctx, verbose, config):
    """"""
    # configure root logger to print to STDERR
    loglevel = LEVELS.get(min(verbose, 3))
    configure_stream(level=loglevel)
    ctx.obj = {}
    if config and os.path.exists(config):
        ctx.obj = yaml.load(open(config, 'r')) or {}
        ctx.obj['config_path'] = config
    # launch the command line interface
    logger.debug('Booting up command line interface')
1032462
def cli(variant_file, vep, split, outfile, verbose, silent, check_info,
        allele_symbol, logfile, loglevel):
    
    from vcf_parser import logger, init_log

    if not loglevel:
        if verbose:
            loglevel = 'INFO'

    init_log(logger, logfile, loglevel)
    nr_of_variants = 0
    start = datetime.now()
    
    # with open(variant_file, 'r', encoding="utf-8") as f:
    #     for line in f:
    #         if not line.startswith('#'):
    #             nr_of_variants += 1
    
    if variant_file == '-':
        logger.info("Start parsing variants from stdin")
        my_parser = VCFParser(
            fsock=sys.stdin, 
            split_variants=split,
            check_info=check_info, 
            allele_symbol=allele_symbol
        )
    else:
        logger.info("Start parsing variants from file {0}".format(variant_file))
        my_parser = VCFParser(
            infile = variant_file,
            split_variants=split, 
            check_info=check_info, 
            allele_symbol=allele_symbol
        )

    if outfile:
        f = open(outfile, 'w', encoding='utf-8')
        logger.info("Printing vcf to file {0}".format(outfile))
    
    if not silent:
        logger.info("Printing vcf to stdout")
    else:
        logger.info("Skip printing since silent is active")
    
    for line in my_parser.metadata.print_header():
        if outfile:
            f.write(line+'\n')
        else:
            if not silent:
                print(line)
    try:
        for variant in my_parser:
            variant_line = '\t'.join([variant[head] for head in my_parser.header])
            if outfile:
                f.write(variant_line + '\n')
            else:
                if not silent:
                    print(variant_line)
            nr_of_variants += 1
    except SyntaxError as e:
        print(e)

    logger.info('Number of variants: {0}'.format(nr_of_variants))
    logger.info('Time to parse file: {0}'.format(str(datetime.now() - start)))
458110
def get_genomes(cfg):
    
    if cfg['host']=='homo_sapiens':
        contigs=['1',
                 '2',
                 '3',
                 '4',
                 '5',
                 '6',
                 '7',
                 '8',
                 '9',
                 '10',
                 '11',
                 '12',
                 '13',
                 '14',
                 '15',
                 '16',
                 '17',
                 '18',
                 '19',
                 '20',
                 '21',
                 '22',
                 'X','Y']
    else:
        runbashcmd(f"pyensembl install --reference-name {cfg['genomeassembly']} --release {cfg['genomerelease']} --species {cfg['host']}")

        import pyensembl
        ensembl = pyensembl.EnsemblRelease(species=pyensembl.species.Species.register(
                                            latin_name=cfg['host'],
                                            synonyms=[cfg['host']],
                                            reference_assemblies={
                                                cfg['genomeassembly']: (cfg['genomerelease'], cfg['genomerelease']),
                                            }),release=cfg['genomerelease'])
        contig_mito=['MTDNA','MITO','MT']
        contigs=[c for c in ensembl.contigs() if ((not '.' in c) and (len(c)<5) and (c not in contig_mito))]    
    if len(contigs)==0:
        contigs=[c for c in ensembl.contigs()]
        # logging.error([c for c in ensembl.contigs()])
        # logging.error('no contigs identified by pyensembl; aborting')
        # sys.exit(0)
    logging.info(f"{len(contigs)} contigs/chromosomes in the genome")
    logging.info(contigs)
    # raw genome next
    if 'human' in cfg['host'].lower():
        cfg['host']='homo_sapiens'
    if 'yeast' in cfg['host'].lower():
        cfg['host']='saccharomyces_cerevisiae'
    host_="_".join(s for s in cfg['host'].split('_')).capitalize()
    if 'GRCh37' in cfg['genomeassembly']:    
        ensembl_fastad=f"pub/grch37/update/fasta/{cfg['host']}/dna/"
    else:
        ensembl_fastad=f"pub/release-{cfg['genomerelease']}/fasta/{cfg['host']}/dna/"        
    genome_fastad=f"{cfg['genomed']}/{ensembl_fastad}"
    cfg['genomep']=f'{genome_fastad}/genome.fa'.format()
    if not exists(cfg['genomep']):
        logging.error(f"not found: {cfg['genomep']}")
        ifdlref = input(f"Genome file are not there at {genome_fastad}.\n Download?[Y/n]: ")
        if ifdlref=='Y':
        # #FIXME download contigs and cat and get index, sizes
            for contig in contigs:
                fn=f"{cfg['host'].capitalize()}.{cfg['genomeassembly']}.dna_sm.chromosome.{contig}.fa.gz"
                fp=f'{ensembl_fastad}/{fn}'
                if not exists(f"{cfg['genomed']}{fp.replace('.gz','')}"):
                    if not exists(f"{cfg['genomed']}{fp}"):
                        logging.info(f'downloading {fp}')
                        cmd=f"wget -q -x -nH ftp://ftp.ensembl.org/{fp} -P {cfg['genomed']}"
                        try:
                            runbashcmd(cmd,test=cfg['test'])
                        except:
                            fn=f"{cfg['host'].capitalize()}.{cfg['genomeassembly']}.dna_sm.toplevel.fa.gz"
                            fp='{}/{}'.format(ensembl_fastad,fn)                        
                            if not exists(fp):
                                cmd=f"wget -q -x -nH ftp://ftp.ensembl.org/{fp} -P {cfg['genomed']}"
                                # print(cmd)
                                runbashcmd(cmd,test=cfg['test'])
                                break
                #                 break
            # make the fa ready
            if not exists(cfg['genomep']):
                cmd=f"gunzip {genome_fastad}*.fa.gz;cat {genome_fastad}/*.fa > {genome_fastad}/genome.fa;"
                runbashcmd(cmd,test=cfg['test'])
        else:
            logging.error('abort')
            sys.exit(1)
    if not exists(cfg['genomep']+'.bwt'):
        cmd=f"{cfg['bwa']} index {cfg['genomep']}"
        runbashcmd(cmd,test=cfg['test'])
    else:        
        logging.info('bwa index is present')
    if not exists(cfg['genomep']+'.fai'):
        cmd=f"{cfg['samtools']} faidx {cfg['genomep']}"
        runbashcmd(cmd,test=cfg['test'])
    else:
        logging.info('samtools index is present')
    if not exists(cfg['genomep']+'.sizes'):
        cmd=f"cut -f1,2 {cfg['genomep']}.fai > {cfg['genomep']}.sizes"            
        runbashcmd(cmd,test=cfg['test'])
    else:
        logging.info('sizes of contigs are present')

    if 'GRCh37' in cfg['genomeassembly']:    
#     ftp://ftp.ensembl.org/pub/grch37/update/gff3/homo_sapiens/Homo_sapiens.GRCh37.87.gff3.gz
        ensembl_gff3d=f"pub/grch37/update/gff3/{cfg['host']}/"    
    else:
        ensembl_gff3d=f"pub/release-{cfg['genomerelease']}/gff3/{cfg['host']}/"    
    
    genome_gff3d=f"{cfg['genomed']}/{ensembl_gff3d}"
    cfg['genomegffp']=f'{genome_gff3d}/genome.gff3'
    if not exists(cfg['genomegffp']):
        logging.error('not found: {}'.format(cfg['genomegffp']))
        ifdlref = input("Download genome annotations at {}?[Y/n]: ".format(genome_gff3d))
        if ifdlref=='Y':
        # #FIXME download contigs and cat and get index, sizes
            fn=f"{cfg['host'].capitalize()}.{cfg['genomeassembly']}.{cfg['genomerelease']}.gff3.gz"
            fp=f"{ensembl_gff3d}/{fn}"
            if not exists(fp):
                cmd=f"wget -x -nH ftp://ftp.ensembl.org/{fp} -P {cfg['genomed']}"
                runbashcmd(cmd,test=cfg['test'])
            # move to genome.gff3
                cmd=f"cp {genome_gff3d}/{fn} {cfg['genomegffp']}"
                runbashcmd(cmd,test=cfg['test'])

        else:
            logging.error('abort')
            sys.exit(1)
    logging.info('genomes are installed!')
    return cfg
273825
def _square_batch_bcbio_variation(data, region, bam_files, vrn_files, out_file,
                                  todo="square"):
    
    ref_file = tz.get_in(("reference", "fasta", "base"), data)
    cores = tz.get_in(("config", "algorithm", "num_cores"), data, 1)
    resources = config_utils.get_resources("bcbio-variation-recall", data["config"])
    # adjust memory by cores but leave room for run program memory
    memcores = int(math.ceil(float(cores) / 5.0))
    jvm_opts = config_utils.adjust_opts(resources.get("jvm_opts", ["-Xms250m", "-Xmx2g"]),
                                        {"algorithm": {"memory_adjust": {"direction": "increase",
                                                                         "magnitude": memcores}}})
    # Write unique VCFs and BAMs to input file
    input_file = "%s-inputs.txt" % os.path.splitext(out_file)[0]
    with open(input_file, "w") as out_handle:
        out_handle.write("\n".join(sorted(list(set(vrn_files)))) + "\n")
        if todo == "square":
            out_handle.write("\n".join(sorted(list(set(bam_files)))) + "\n")
    variantcaller = tz.get_in(("config", "algorithm", "jointcaller"), data).replace("-joint", "")
    cmd = ["bcbio-variation-recall", todo] + jvm_opts + broad.get_default_jvm_opts() + \
          ["-c", cores, "-r", bamprep.region_to_gatk(region)]
    if todo == "square":
        cmd += ["--caller", variantcaller]
    cmd += [out_file, ref_file, input_file]
    bcbio_env = utils.get_bcbio_env()
    cmd = " ".join(str(x) for x in cmd)
    do.run(cmd, "%s in region: %s" % (cmd, bamprep.region_to_gatk(region)), env=bcbio_env)
    return out_file
52131
def _variant_filtration(in_file, ref_file, vrn_files, data, filter_type,
                        hard_filter_fn):
    
    # Algorithms multiplied by number of input files to check for large enough sample sizes
    algs = [data["config"]["algorithm"]] * len(data.get("vrn_files", [1]))
    if (not config_utils.use_vqsr(algs, in_file) or
          _already_cutoff_filtered(in_file, filter_type)):
        logger.info("Skipping VQSR, using cutoff-based filers: we don't have whole genome input data")
        return hard_filter_fn(in_file, data)
    elif not _have_training_data(vrn_files):
        logger.info("Skipping VQSR, using cutoff-based filers: genome build does not have sufficient training data")
        return hard_filter_fn(in_file, data)
    else:
        sensitivities = {"INDEL": "98.0", "SNP": "99.97"}
        recal_file, tranches_file = _run_vqsr(in_file, ref_file, vrn_files,
                                              sensitivities[filter_type], filter_type, data)
        if recal_file is None:  # VQSR failed
            logger.info("VQSR failed due to lack of training data. Using cutoff-based soft filtering.")
            return hard_filter_fn(in_file, data)
        else:
            return _apply_vqsr(in_file, ref_file, recal_file, tranches_file,
                               sensitivities[filter_type], filter_type, data)
751401
def _run_genotype_gvcfs_genomicsdb(genomics_db, region, out_file, data):
    
    if not utils.file_exists(out_file):
        with file_transaction(data, out_file) as tx_out_file:
            broad_runner = broad.runner_from_config(data["config"])
            params = ["-T", "GenotypeGVCFs",
                      "--variant", "gendb://%s" % genomics_db,
                      "-R", dd.get_ref_file(data),
                      "--output", tx_out_file,
                      "-L", bamprep.region_to_gatk(region)]
            params += ["-ploidy", str(ploidy.get_ploidy([data], region))]
            # Avoid slow genotyping runtimes with improved quality score calculation in GATK4
            # https://gatkforums.broadinstitute.org/gatk/discussion/11471/performance-troubleshooting-tips-for-genotypegvcfs/p1
            params += ["--use-new-qual-calculator"]
            resources = config_utils.get_resources("gatk", data["config"])
            params += [str(x) for x in resources.get("options", [])]
            cores = dd.get_cores(data)
            memscale = {"magnitude": 0.9 * cores, "direction": "increase"} if cores > 1 else None
            broad_runner.run_gatk(params, memscale=memscale)
    return vcfutils.bgzip_and_index(out_file, data["config"])
849701
def _variant_checkpoints(samples):
    
    checkpoints = {}
    checkpoints["vc"] = any([dd.get_variantcaller(d) or d.get("vrn_file") for d in samples])
    checkpoints["sv"] = any([dd.get_svcaller(d) for d in samples])
    checkpoints["jointvc"] = any([(dd.get_jointcaller(d) or "gvcf" in dd.get_tools_on(d))
                                  for d in samples])
    checkpoints["hla"] = any([dd.get_hlacaller(d) for d in samples])
    checkpoints["align"] = any([(dd.get_aligner(d) or dd.get_bam_clean(d)) for d in samples])
    checkpoints["align_split"] = not all([(dd.get_align_split_size(d) is False or
                                           not dd.get_aligner(d))
                                          for d in samples])
    checkpoints["archive"] = any([dd.get_archive(d) for d in samples])
    checkpoints["umi"] = any([dd.get_umi_consensus(d) for d in samples])
    checkpoints["ensemble"] = any([dd.get_ensemble(d) for d in samples])
    checkpoints["cancer"] = any(dd.get_phenotype(d) in ["tumor"] for d in samples)
    return checkpoints
542809
def parallel_combine_variants(orig_files, out_file, ref_file, config, run_parallel):
    
    file_key = "vcf_files"
    def split_by_region(data):
        base, ext = utils.splitext_plus(os.path.basename(out_file))
        args = []
        for region in [x.name for x in ref.file_contigs(ref_file, config)]:
            region_out = os.path.join(os.path.dirname(out_file), "%s-regions" % base,
                                      "%s-%s%s" % (base, region, ext))
            utils.safe_makedir(os.path.dirname(region_out))
            args.append((region_out, ref_file, config, region))
        return out_file, args
    config = copy.deepcopy(config)
    config["file_key"] = file_key
    prep_files = run_multicore(p_bgzip_and_index, [[x, config] for x in orig_files], config)
    items = [[{file_key: prep_files}]]
    parallel_split_combine(items, split_by_region, run_parallel,
                           "merge_variant_files", "concat_variant_files",
                           file_key, ["region", "sam_ref", "config"], split_outfile_i=0)
    return out_file
603006
def _balanced_strand_gen(base_aligns, limit):
    ''''''
    predicate=lambda align: align.is_reverse
    gen1, gen2 = itertools.tee((predicate(align), align) for align in base_aligns)
    for_strand_gen = (align for pred, align in gen1 if not pred)
    rev_strand_gen = (align for pred, align in gen2 if pred)
    combined_gen = iter_zip(for_strand_gen, rev_strand_gen)
    interleaved_gen = (align for for_rev_aligns in combined_gen for align in for_rev_aligns)
    return itertools.islice(interleaved_gen, limit)
750060
def run(align_bams, items, ref_file, assoc_files, region=None, out_file=None):
    
    paired = vcfutils.get_paired_bams(align_bams, items)
    assert paired and not paired.normal_bam, ("Pisces supports tumor-only variant calling: %s" %
                                              (",".join([dd.get_sample_name(d) for d in items])))
    vrs = bedutils.population_variant_regions(items)
    target = shared.subset_variant_regions(vrs, region,
                                            out_file, items=items, do_merge=True)
    min_af = float(dd.get_min_allele_fraction(paired.tumor_data)) / 100.0
    if not utils.file_exists(out_file):
        base_out_name = utils.splitext_plus(os.path.basename(paired.tumor_bam))[0]
        raw_file = "%s.vcf" % utils.splitext_plus(out_file)[0]
        with file_transaction(paired.tumor_data, raw_file) as tx_out_file:
            ref_dir = _prep_genome(os.path.dirname(tx_out_file), paired.tumor_data)
            out_dir = os.path.dirname(tx_out_file)
            cores = dd.get_num_cores(paired.tumor_data)
            emit_min_af = min_af / 10.0
            cmd = ("pisces --bampaths {paired.tumor_bam} --genomepaths {ref_dir} --intervalpaths {target} "
                   "--maxthreads {cores} --minvf {emit_min_af} --vffilter {min_af} "
                   "--ploidy somatic --gvcf false -o {out_dir}")
            # Recommended filtering for low frequency indels
            # https://github.com/bcbio/bcbio-nextgen/commit/49d0cbb1f6dcbea629c63749e2f9813bd06dcee3#commitcomment-29765373
            cmd += " -RMxNFilter 5,9,0.35"
            # For low frequency UMI tagged variants, set higher variant thresholds
            # https://github.com/Illumina/Pisces/issues/14#issuecomment-399756862
            if min_af < (1.0 / 100.0):
                cmd += " --minbasecallquality 30"
            do.run(cmd.format(**locals()), "Pisces tumor-only somatic calling")
            shutil.move(os.path.join(out_dir, "%s.vcf" % base_out_name),
                        tx_out_file)
        vcfutils.bgzip_and_index(raw_file, paired.tumor_data["config"],
                                 prep_cmd="sed 's#%s.bam#%s#' | %s" %
                                 (base_out_name, dd.get_sample_name(paired.tumor_data),
                                  vcfutils.add_contig_to_header_cl(dd.get_ref_file(paired.tumor_data), out_file)))
    return vcfutils.bgzip_and_index(out_file, paired.tumor_data["config"])
267454
def _get_variant_silent(parser, variant):
    """"""
    prev_log = config.LOG_NOT_FOUND
    config.LOG_NOT_FOUND = False
    results = parser.get_variant_genotypes(variant)
    config.LOG_NOT_FOUND = prev_log
    return results
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
1072342
def make_V_and_J_mask_mapping(self, genV, genJ):
        
        #construct mapping between allele/gene names and index for custom V_usage_masks
        V_allele_names = [V[0] for V in genV]
        V_mask_mapping = {}
        for v in set([x.split('*')[0] for x in V_allele_names]):
            V_mask_mapping[v] = []
        for v in set(['V'.join((x.split('*')[0]).split('V')[1:]) for x in V_allele_names]):
            V_mask_mapping[v] = []
        for i, v in enumerate(V_allele_names):
            V_mask_mapping[v] = [i]
            V_mask_mapping['V'.join((v.split('*')[0]).split('V')[1:])].append(i)
            V_mask_mapping[v.split('*')[0]].append(i)      
        
        #construct mapping between allele/gene names and index for custom J_usage_masks
        J_allele_names = [J[0] for J in genJ]
        J_mask_mapping = {}
        for j in set([x.split('*')[0] for x in J_allele_names]):
            J_mask_mapping[j] = []
        for j in set(['J'.join((x.split('*')[0]).split('J')[1:]) for x in J_allele_names]):
            J_mask_mapping[j] = []
        for i, j in enumerate(J_allele_names):
            J_mask_mapping[j] = [i]
            J_mask_mapping['J'.join((j.split('*')[0]).split('J')[1:])].append(i)
            J_mask_mapping[j.split('*')[0]].append(i)
            
        self.V_allele_names = V_allele_names
        self.V_mask_mapping = V_mask_mapping
        self.J_allele_names = J_allele_names
        self.J_mask_mapping = J_mask_mapping
724720
def add_reference_resources(data, remote_retriever=None):
    
    aligner = data["config"]["algorithm"].get("aligner", None)
    if remote_retriever:
        data["reference"] = remote_retriever.get_refs(data["genome_build"],
                                                      alignment.get_aligner_with_aliases(aligner, data),
                                                      data["config"])
    else:
        data["reference"] = genome.get_refs(data["genome_build"], alignment.get_aligner_with_aliases(aligner, data),
                                            data["dirs"]["galaxy"], data)
        _check_ref_files(data["reference"], data)
    # back compatible `sam_ref` target
    data["sam_ref"] = utils.get_in(data, ("reference", "fasta", "base"))
    ref_loc = utils.get_in(data, ("config", "resources", "species", "dir"),
                           utils.get_in(data, ("reference", "fasta", "base")))
    if remote_retriever:
        data = remote_retriever.get_resources(data["genome_build"], ref_loc, data)
    else:
        data["genome_resources"] = genome.get_resources(data["genome_build"], ref_loc, data)
    data["genome_resources"] = genome.add_required_resources(data["genome_resources"])
    if effects.get_type(data) == "snpeff" and "snpeff" not in data["reference"]:
        data["reference"]["snpeff"] = effects.get_snpeff_files(data)
    if "genome_context" not in data["reference"]:
        data["reference"]["genome_context"] = annotation.get_context_files(data)
    if "viral" not in data["reference"]:
        data["reference"]["viral"] = viral.get_files(data)
    if not data["reference"]["viral"]:
        data["reference"]["viral"] = None
    if "versions" not in data["reference"]:
        data["reference"]["versions"] = _get_data_versions(data)

    data = _fill_validation_targets(data)
    data = _fill_prioritization_targets(data)
    data = _fill_capture_regions(data)
    # Re-enable when we have ability to re-define gemini configuration directory
    if False:
        data["reference"]["gemini"] = population.get_gemini_files(data)
    return data
371941
def load(ctx, variant_source, family_file, family_type, root):
    
    root = root or ctx.obj.get('root') or os.path.expanduser("~/.puzzle")

    if os.path.isfile(root):
        logger.error("'root' can't be a file")
        ctx.abort()

    logger.info("Root directory is: {}".format(root))

    db_path = os.path.join(root, 'puzzle_db.sqlite3')
    logger.info("db path is: {}".format(db_path))
    
    if not os.path.exists(db_path):
        logger.warn("database not initialized, run 'puzzle init'")
        ctx.abort()


    if not os.path.isfile(variant_source):
        logger.error("Variant source has to be a file")
        ctx.abort()
    
    mode = get_file_type(variant_source)
    if mode == 'unknown':
        logger.error("Unknown file type")
        ctx.abort()
    #Test if gemini is installed
    elif mode == 'gemini':
        logger.debug("Initialzing GEMINI plugin")
        if not GEMINI:
            logger.error("Need to have gemini installed to use gemini plugin")
            ctx.abort()

    logger.debug('Set puzzle backend to {0}'.format(mode))
    
    variant_type = get_variant_type(variant_source)
    logger.debug('Set variant type to {0}'.format(variant_type))
    
    cases = get_cases(
        variant_source=variant_source,
        case_lines=family_file, 
        case_type=family_type, 
        variant_type=variant_type, 
        variant_mode=mode
    )
    
    if len(cases) == 0:
        logger.warning("No cases found")
        ctx.abort()

    logger.info("Initializing sqlite plugin")
    store = SqlStore(db_path)

    for case_obj in cases:
        if store.case(case_obj.case_id) is not None:
            logger.warn("{} already exists in the database"
                        .format(case_obj.case_id))
            continue
        # extract case information
        logger.debug("adding case: {} to puzzle db".format(case_obj.case_id))
        store.add_case(case_obj, vtype=variant_type, mode=mode)
1344741
def run_haplotyper(align_bams, items, ref_file, assoc_files,
                     region=None, out_file=None):
    
    if out_file is None:
        out_file = "%s-variants.vcf.gz" % utils.splitext_plus(align_bams[0])[0]
    if not utils.file_exists(out_file):
        variant_regions = bedutils.population_variant_regions(items, merged=True)
        interval = _get_interval(variant_regions, region, out_file, items)
        with file_transaction(items[0], out_file) as tx_out_file:
            dbsnp = "--dbsnp %s" % (assoc_files.get("dbsnp")) if "dbsnp" in assoc_files else ""
            bams = " ".join(["-i %s" % x for x in align_bams])
            license = license_export(items[0])
            cores = dd.get_num_cores(items[0])
            out_mode = "--emit_mode gvcf" if joint.want_gvcf(items) else ""
            cmd = ("{license}sentieon driver -t {cores} -r {ref_file} "
                   "{bams} {interval} --algo Haplotyper {out_mode} {dbsnp} {tx_out_file}")
            do.run(cmd.format(**locals()), "Sentieon Haplotyper")
    return out_file
149602
def store_disorder(self, sc=None, force_rerun=False):
        """"""
        log.info('Loading sequences to reference GEM-PRO...')
        from random import shuffle
        g_ids = [g.id for g in self.reference_gempro.functional_genes]
        shuffle(g_ids)

        def _store_disorder_sc(g_id, outdir=self.sequences_by_gene_dir,
                               g_to_pickle=self.gene_protein_pickles, force_rerun=force_rerun):
            """Load orthologous strain sequences to reference Protein object, save as new pickle"""
            import ssbio.utils
            import ssbio.io
            import os.path as op

            protein_seqs_pickle_path = op.join(outdir, '{}_protein_withseqs_dis.pckl'.format(g_id))

            if ssbio.utils.force_rerun(flag=force_rerun, outfile=protein_seqs_pickle_path):
                protein_pickle_path = g_to_pickle[g_id]
                protein_pickle = ssbio.io.load_pickle(protein_pickle_path)
                protein_pickle.get_all_disorder_predictions(representative_only=False)
                protein_pickle.save_pickle(outfile=protein_seqs_pickle_path)

            return g_id, protein_seqs_pickle_path

        if sc:
            genes_rdd = sc.parallelize(g_ids)
            result = genes_rdd.map(_store_disorder_sc).collect()
        else:
            result = []
            for g in tqdm(g_ids):
                result.append(self._load_sequences_to_reference_gene(g, force_rerun))

        log.info('Storing paths to new Protein objects in self.gene_protein_pickles...')
        for g_id, protein_pickle in result:
            self.gene_protein_pickles[g_id] = protein_pickle