Maps RNA-Seq reads to a reference genome.

    job_vars: tuple     Tuple of dictionaries: input_args and ids
555362
def mapsplice(job, job_vars):
    
    # Unpack variables
    input_args, ids = job_vars
    work_dir = job.fileStore.getLocalTempDir()
    cores = input_args['cpu_count']
    sudo = input_args['sudo']
    single_end_reads = input_args['single_end_reads']
    files_to_delete = ['R1.fastq']
    # I/O
    return_input_paths(job, work_dir, ids, 'ebwt.zip', 'chromosomes.zip')
    if single_end_reads:
        return_input_paths(job, work_dir, ids, 'R1.fastq')
    else:
        return_input_paths(job, work_dir, ids, 'R1.fastq', 'R2.fastq')
        files_to_delete.extend(['R2.fastq'])
    for fname in ['chromosomes.zip', 'ebwt.zip']:
        subprocess.check_call(['unzip', '-o', os.path.join(work_dir, fname), '-d', work_dir])
    # Command and call
    parameters = ['-p', str(cores),
                  '-s', '25',
                  '--bam',
                  '--min-map-len', '50',
                  '-x', '/data/ebwt',
                  '-c', '/data/chromosomes',
                  '-1', '/data/R1.fastq',
                  '-o', '/data']
    if not single_end_reads:
        parameters.extend(['-2', '/data/R2.fastq'])
    docker_call(tool='quay.io/ucsc_cgl/mapsplice:2.1.8--dd5ac549b95eb3e5d166a5e310417ef13651994e',
                tool_parameters=parameters, work_dir=work_dir, sudo=sudo)
    # Write to FileStore
    for fname in ['alignments.bam', 'stats.txt']:
        ids[fname] = job.fileStore.writeGlobalFile(os.path.join(work_dir, fname))
    for fname in files_to_delete:
        job.fileStore.deleteGlobalFile(ids[fname])
    # Run child job
    # map_id = job.addChildJobFn(mapping_stats, job_vars).rv()
    if input_args['upload_bam_to_s3'] and input_args['s3_dir']:
        job.addChildJobFn(upload_bam_to_s3, job_vars)
    output_ids = job.addChildJobFn(add_read_groups, job_vars, disk='30 G').rv()
    return output_ids
8632
def _load_sequences_to_strain(self, strain_id, force_rerun=False):
        """"""
        gp_seqs_path = op.join(self.model_dir, '{}_gp_withseqs.pckl'.format(strain_id))

        if ssbio.utils.force_rerun(flag=force_rerun, outfile=gp_seqs_path):
            gp_noseqs = ssbio.io.load_pickle(self.strain_infodict[strain_id]['gp_noseqs_path'])
            strain_sequences = SeqIO.index(self.strain_infodict[strain_id]['genome_path'], 'fasta')
            for strain_gene in gp_noseqs.functional_genes:
                # Pull the gene ID of the strain from the orthology matrix
                strain_gene_key = self.df_orthology_matrix.at[strain_gene.id, strain_id]
                # Load into the strain GEM-PRO
                new_id = '{}_{}'.format(strain_gene.id, strain_id)
                if strain_gene.protein.sequences.has_id(new_id):
                    continue
                strain_gene.protein.load_manual_sequence(seq=strain_sequences[strain_gene_key], ident=new_id,
                                                         set_as_representative=True)
            gp_noseqs.save_pickle(outfile=gp_seqs_path)

        return strain_id, gp_seqs_path
1159976
def liftover_to_genome(pass_pos, gtf):
    """"""

    fixed_pos = []
    for pos in pass_pos:
        if pos["chrom"] not in gtf:
            continue
        db_pos = gtf[pos["chrom"]][0]
        mut = _parse_mut(pos["sv"])
        print([db_pos, pos])
        if db_pos[3] == "+":
            pos['pre_pos'] = db_pos[1] + pos["pre_pos"] + 1
        else:
            pos['pre_pos'] = db_pos[2] - (pos["pre_pos"] - 1)
        pos['chrom'] = db_pos[0]
        pos['nt'] = list(mut[0])
        fixed_pos.append(pos)

    _print_header(fixed_pos)
    for pos in fixed_pos:
        print_vcf(pos)
255278
def run_vardict(align_bams, items, ref_file, assoc_files, region=None,
                out_file=None):
    
    items = shared.add_highdepth_genome_exclusion(items)
    if vcfutils.is_paired_analysis(align_bams, items):
        call_file = _run_vardict_paired(align_bams, items, ref_file,
                                        assoc_files, region, out_file)
    else:
        vcfutils.check_paired_problems(items)
        call_file = _run_vardict_caller(align_bams, items, ref_file,
                                        assoc_files, region, out_file)
    return call_file
1323057
def InputCodon(seq_length, ignore_stop_codons=True, name=None, **kwargs):
    
    if ignore_stop_codons:
        vocab = CODONS
    else:
        vocab = CODONS + STOP_CODONS

    assert seq_length % 3 == 0
    return Input((seq_length / 3, len(vocab)), name=name, **kwargs)
151157
def batch(samples):
    
    samples = [utils.to_single_data(x) for x in samples]
    sample_order = [dd.get_sample_name(x) for x in samples]
    batch_groups = collections.defaultdict(list)
    for data in samples:
        batch_samples = tuple(data.get("batch_samples", [dd.get_sample_name(data)]))
        batch_groups[(batch_samples, dd.get_phenotype(data))].append(data)

    out = []
    for (batch_samples, phenotype), gsamples in batch_groups.items():
        if len(gsamples) > 1:
            batches = set([])
            for d in gsamples:
                batches |= set(dd.get_batches(d))
            gsamples.sort(key=dd.get_variantcaller_order)
            cur = copy.deepcopy(gsamples[0])
            cur.update({"batch_id": sorted(list(batches))[0] if batches else "_".join(batch_samples),
                        "batch_samples": batch_samples,
                        "variants": {"variantcallers": [dd.get_variantcaller(d) for d in gsamples],
                                     "calls": [d.get("vrn_file") for d in gsamples]}})
            out.append(cur)

    def by_original_order(d):
        return min([sample_order.index(s) for s in d["batch_samples"] if s in sample_order])
    return sorted(out, key=by_original_order)
799136
def sample_chromosomes(job, genome_fai_file):
    
    work_dir = os.getcwd()
    genome_fai = untargz(job.fileStore.readGlobalFile(genome_fai_file), work_dir)
    return chromosomes_from_fai(genome_fai)
248485
def strainer(sequencepath):
    
    metadata_list = list()
    assert os.path.isdir(sequencepath), 'Cannot locate sequence path as specified: {}' \
        .format(sequencepath)
    # Get the sequences in the sequences folder into a list. Note that they must have a file extension that
    # begins with .fa
    strains = sorted(glob.glob(os.path.join(sequencepath, '*.fa*')))
    # Populate the metadata object. This object will be populated to mirror the objects created in the
    # genome assembly pipeline. This way this script will be able to be used as a stand-alone, or as part
    # of a pipeline
    assert strains, 'Could not find any files with an extension starting with "fa" in {}. Please check ' \
                    'to ensure that your sequence path is correct'.format(sequencepath)
    for sample in strains:
        # Create the object
        metadata = MetadataObject()
        # Set the base file name of the sequence. Just remove the file extension
        filename = os.path.splitext(os.path.split(sample)[1])[0]
        # Set the .name attribute to be the file name
        metadata.name = filename
        # Create the .general attribute
        metadata.general = GenObject()
        metadata.commands = GenObject()
        metadata.general.outputdirectory = os.path.join(sequencepath, filename)
        # Set the .general.bestassembly file to be the name and path of the sequence file
        metadata.general.bestassemblyfile = os.path.join(metadata.general.outputdirectory, '{sn}.fasta'
                                                         .format(sn=filename))
        make_path(metadata.general.outputdirectory)
        # Create a symlink to the directory
        relative_symlink(sample,
                         metadata.general.outputdirectory)
        metadata.general.logout = os.path.join(metadata.general.outputdirectory, 'out')
        metadata.general.logerr = os.path.join(metadata.general.outputdirectory, 'err')
        # Append the metadata for each sample to the list of samples
        metadata_list.append(metadata)
    return strains, metadata_list
1150902
def main(args):
   """"""
   txome = Transcriptome()
   #read the reference gpd if one is gven
   if args.reference: 
      rinf = None
      if re.search('\.gz$',args.reference):
         rinf = gzip.open(args.reference)
      else:
         rinf = open(args.reference)
      sys.stderr.write("Reading in reference\n")
      z = 0
      # populate txome with reference transcripts for each chromosome
      for line in rinf:
         z += 1
         gpd = GPD(line)
         gpd.set_payload(z)
         if z%100 == 0:  sys.stderr.write(str(z)+"          \r")
         txome.add_transcript(gpd)
      rinf.close()
      sys.stderr.write(str(z)+"          \r")
      sys.stderr.write("\n")
   txome.sort_transcripts()
   sys.stderr.write("Buffering mappings\n")
   inf = sys.stdin
   if args.input != '-':
      if re.search('\.gz$',args.input):
         inf = gzip.open(args.input)
      else:
         inf = open(args.input)
   tof = gzip.open(args.tempdir+'/reads.gpd.gz','w')
   for line in inf: tof.write(line.rstrip()+"\n")
   tof.close()
  
   sys.stderr.write("1. Process by overlapping locus.\n")
   annotated_reads(txome,gzip.open(args.tempdir+'/reads.gpd.gz'),'initial',args)

   #of.close()

   """Now sequences we should be able to annotate as partial have
      consensus transcripts.  Lets put those to annotation scrutiny again.
   """
   sys.stderr.write("2. Reannotated partially annotated by overlapping locus.\n")
   annotated_reads(txome,gzip.open(args.tempdir+'/initial/partial_annotated_multiexon.sorted.gpd.gz'),'repartial',args)


   txs = {}
   for transcript in txome.transcripts:
      txs[transcript.name] = transcript
   detected = {}
   tinf = gzip.open(args.tempdir+'/initial/annotated.txt.gz')
   for line in tinf:
      f = line.rstrip().split("\t")
      detected[f[1]] = True
   tinf.close()
   tinf = gzip.open(args.tempdir+'/repartial/annotated.txt.gz')
   remove_partial = {}
   for line in tinf:
      f = line.rstrip().split("\t")
      detected[f[1]] = True
      remove_partial[f[1]] = True
   tinf.close()
   tof = gzip.open(args.tempdir+'/candidate.gpd.gz','w')
   for name in detected:
      tof.write(txs[name].get_gpd_line()+"\n")
   tinf = gzip.open(args.tempdir+'/initial/partial_annotated_multiexon.sorted.gpd.gz')
   for line in tinf:
      f = line.rstrip().split("\t")
      if f[1] not in remove_partial:
         tof.write(line)
   tinf = gzip.open(args.tempdir+'/initial/unannotated_singleexon.sorted.gpd.gz')
   for line in tinf:
      f = line.rstrip().split("\t")
      f[0] = str(uuid.uuid4())
      tof.write("\t".join(f)+"\n")
   tinf.close()
   tof.close()
   sort_gpd(args.tempdir+'/candidate.gpd.gz',args.tempdir+'/candidate.sorted.gpd.gz',args)
   """We can ignore the partial annotations that have been detected"""
   ntxome = Transcriptome()
   tinf = gzip.open(args.tempdir+'/candidate.sorted.gpd.gz')
   for line in tinf: ntxome.add_transcript(GPD(line))     
   annotated_reads(ntxome,gzip.open(args.tempdir+'/initial/unannotated_multiexon.sorted.gpd.gz'),'unannotated',args)
   """now we know which unannotated reads actually have annotations"""
   tinf.close()

   tof = gzip.open(args.tempdir+'/final.gpd.gz','w')
   tinf = gzip.open(args.tempdir+'/candidate.sorted.gpd.gz')   
   for line in tinf: tof.write(line)
   tinf.close()
   tinf = gzip.open(args.tempdir+'/unannotated/unannotated_multiexon.sorted.gpd.gz')   
   for line in tinf: tof.write(line)
   tinf.close()
   tinf = gzip.open(args.tempdir+'/unannotated/partial_annotated_multiexon.sorted.gpd.gz')   
   for line in tinf: tof.write(line)
   tinf.close()
   tof.close()
   sort_gpd(args.tempdir+'/final.gpd.gz',args.tempdir+'/final.sorted.gpd.gz',args)
   """Prepare outputs"""
   of = sys.stdout
   if args.output:
      if re.search('\.gz$',args.output):
         of = gzip.open(args.output,'w')
      else:
         of = open(args.output,'w')
   tinf = gzip.open(args.tempdir+'/final.sorted.gpd.gz')
   for line in tinf: of.write(line)
   tinf.close()
   of.close()
1007355
def sam2rnf(args):
    

    rnftools.mishmash.Source.recode_sam_reads(
        sam_fn=args.sam_fn,
        fastq_rnf_fo=args.fq_fo,
        fai_fo=args.fai_fo,
        genome_id=args.genome_id,
        number_of_read_tuples=10**9,
        simulator_name=args.simulator_name,
        allow_unmapped=args.allow_unmapped,
    )
470748
def referencegenomefinder(self):
        
        # Initialise dictionaries
        referencematch = defaultdict(make_dict)
        referencehits = defaultdict(make_dict)
        # Set the name of the reference profile file
        referencegenomeprofile = '{}rMLST_referenceprofile.json'.format(self.referenceprofilepath)
        # Open the reference profile and load the profile into memory
        with open(referencegenomeprofile) as referencefile:
            referencetypes = json.load(referencefile)
        # Iterate through the samples
        for sample in self.metadata:
            if sample[self.analysistype].reportdir != 'NA':
                # Iterate through the reference genomes in the profile
                for genome in referencetypes:
                    # Initialise the number of identical alleles between the assembly of interest and the
                    # reference genome to 0
                    referencehits[sample.name][genome] = 0
                    # Iterate through all the genes in the analysis
                    for gene in self.bestdict[sample.name]:
                        # If the alleles match between the assembly of interest and the reference genome, increment
                        # the number of matches by 1
                        if list(self.bestdict[sample.name][gene].keys())[0] == referencetypes[genome][gene]:
                            referencematch[sample.name][genome][gene] = 1
                            referencehits[sample.name][genome] += 1
                        else:
                            referencematch[sample.name][genome][gene] = 0
        #
        for sample in self.metadata:
            if sample[self.analysistype].reportdir != 'NA':
                # Get the best number of matches
                # From: https://stackoverflow.com/questions/613183/sort-a-python-dictionary-by-value
                try:
                    matches = sorted(referencehits[sample.name].items(),
                                     key=operator.itemgetter(1), reverse=True)
                    most_matches = matches[0][1]
                    i = 0
                    match_list = list()
                    while matches[i][1] == most_matches:
                        match_list.append(matches[i])
                        i += 1
                    sorted_list = sorted(match_list)
                    sortedmatches = sorted_list[0]
                except IndexError:
                    sortedmatches = (0, 0)
                # If there are fewer matches than the total number of genes in the typing scheme
                if 0 < int(sortedmatches[1]) < len(sample[self.analysistype].allelenames):
                    mismatches = []
                    # Iterate through the gene in the analysis
                    for gene, allele in referencetypes[sortedmatches[0]].items():
                        # Populate :self.referencegenome with the genome name, best reference match, number of matches,
                        # gene, query allele(s), and percent identity
                        percentidentity = '{:.2f}'.format(list(self.bestdict[sample.name][gene].values())[0])
                        self.referencegenome[sample.name][sortedmatches[0]][sortedmatches[1]][gene][list(self.bestdict[
                            sample.name][gene].keys())[0]] = percentidentity
                        if list(self.bestdict[sample.name][gene].keys())[0] != allele:
                            sample[self.analysistype].referencegenome = sortedmatches[0]
                            sample.general.referencegenus = sortedmatches[0].split('_')[0]
                            sample[self.analysistype].referencegenomepath = '{}{}.fa' \
                                .format(self.referenceprofilepath, sortedmatches[0])
                            sample[self.analysistype].matchestoreferencegenome = sortedmatches[1]
                            mismatches.append(({gene: ('{} ({})'.format(list(self.bestdict[sample.name][gene]
                                                                        .keys())[0], allele))}))
                        # sample[self.analysistype].mismatchestoreferencegenome = sorted(mismatches)
                        sample[self.analysistype].mismatchestoreferencegenome = mismatches
                elif sortedmatches == 0:
                    for gene in sample[self.analysistype].allelenames:
                        # Populate the profile of results with 'negative' values for sequence type and sorted matches
                        self.referencegenome[sample.name][sortedmatches[0]][0][gene]['NA'] = 0
                        sample[self.analysistype].referencegenome = 'NA'
                        sample.general.referencegenus = 'NA'
                        sample[self.analysistype].referencegenomepath = 'NA'
                        sample[self.analysistype].matchestoreferencegenome = 0
                        sample[self.analysistype].mismatchestoreferencegenome = [0]
                # Otherwise, the query profile matches the reference profile
                else:
                    for gene in referencetypes[sortedmatches[0]]:
                        # Populate self.referencegenome as above
                        self.referencegenome[sample.name][sortedmatches[0]][sortedmatches[1]][gene][list(self.bestdict[
                            sample.name][gene].keys())[0]] = '{:.2f}'.format(list(self.bestdict[
                                                                                sample.name][gene].values())[0])
                        sample[self.analysistype].referencegenome = sortedmatches[0]
                        sample[self.analysistype].referencegenomepath = '{}{}.fa' \
                            .format(self.referenceprofilepath, sortedmatches[0])
                        sample.general.referencegenus = sortedmatches[0].split('_')[0]
                        sample[self.analysistype].matchestoreferencegenome = sortedmatches[1]
                        sample[self.analysistype].mismatchestoreferencegenome = [0]
        # Print the results to file
        make_path(self.reportpath)
        with open('{}referencegenomes.csv'.format(self.reportpath), 'w') as referencegenomereport:
            row = 'Strain,referencegenome,numberofmatches\n'
            for sample in self.metadata:
                if sample[self.analysistype].reportdir != 'NA':
                    row += '{},{},{}\n'.format(sample.name, sample[self.analysistype].referencegenome,
                                               sample[self.analysistype].matchestoreferencegenome)
            referencegenomereport.write(row)
            dotter()
394069
def profiler(self):
        
        # Initialise variables
        profiledata = defaultdict(make_dict)
        profileset = set()
        genedict = {}
        # Find all the unique profiles to use with a set
        for sample in self.metadata.samples:
            if sample[self.analysistype].profile != 'NA':
                profileset.add(sample[self.analysistype].profile[0])
        # Extract the profiles for each set
        for sequenceprofile in profileset:
            # Clear the list of genes
            genelist = []
            for sample in self.metadata.samples:
                if sequenceprofile == sample[self.analysistype].profile[0]:
                    genelist = [os.path.split(x)[1].split('.')[0] for x in sample[self.analysistype].alleles]
            try:
                # Open the sequence profile file as a dictionary
                profile = DictReader(open(sequenceprofile))
            # Revert to standard comma separated values
            except KeyError:
                # Open the sequence profile file as a dictionary
                profile = DictReader(open(sequenceprofile))
            # Iterate through the rows
            for row in profile:
                # Iterate through the genes
                for gene in genelist:
                    # Add the sequence profile, and type, the gene name and the allele number to the dictionary
                    try:
                        profiledata[sequenceprofile][row['ST']][gene] = row[gene]
                    except KeyError:
                        pass
            # Add the gene list to a dictionary
            genedict[sequenceprofile] = sorted(genelist)
            # Add the profile data, and gene list to each sample
            for sample in self.metadata.samples:
                if sample.general.bestassemblyfile != 'NA':
                    if sequenceprofile == sample[self.analysistype].profile[0]:
                        # Populate the metadata with the profile data
                        sample[self.analysistype].profiledata = profiledata[sample[self.analysistype].profile[0]]
                        # Add the allele directory to a list of directories used in this analysis
                        self.allelefolders.add(sample[self.analysistype].alleledir)
                        dotter()
832532
def build_job_configs(self, args):
        
        job_configs = {}

        ttype = args['ttype']
        (sim_targets_yaml, sim) = NAME_FACTORY.resolve_targetfile(args)
        targets = load_yaml(sim_targets_yaml)

        base_config = dict(ttype=ttype,
                           roi_baseline=args['roi_baseline'],
                           extracopy = args['extracopy'],
                           sim=sim)

        for target_name in targets.keys():
            targetdir = NAME_FACTORY.sim_targetdir(target_type=ttype,
                                                   target_name=target_name,
                                                   sim_name=sim)
            logfile = os.path.join(targetdir, 'copy_base_dir.log')
            job_config = base_config.copy()
            job_config.update(dict(target=target_name,
                                   logfile=logfile))
            job_configs[target_name] = job_config

        return job_configs
321553
def _variant_sv(checkpoints):
    
    if not checkpoints.get("sv"):
        return [], []
    sv = [s("detect_sv", "batch-single",
            [["sv_batch_rec"]],
            [cwlout("sv_rec", "record",
                    fields=[cwlout(["sv", "variantcaller"], ["string", "null"]),
                            cwlout(["sv", "vrn_file"], ["File", "null"], [".tbi"]),
                            cwlout(["sv", "supplemental"], {"type": "array", "items": ["File"]}),
                            cwlout(["svvalidate", "summary"], ["File", "null"]),
                            cwlout("inherit", exclude=[["align_bam"], ["work_bam_plus"],
                                                       ["reference", "snpeff"]])])],
            "bcbio-vc", ["bedtools", "cnvkit", "delly", "duphold", "extract-sv-reads", "gsort",
                         "lumpy-sv;env=python2", "manta;env=python2", "break-point-inspector", "mosdepth", "samtools",
                         "smoove;env=python2", "pysam>=0.13.0",
                         "seq2c", "simple_sv_annotation;env=python2", "survivor", "svtools;env=python2",
                         "svtyper;env=python2",
                         "r=3.5.1", "r-base", "xorg-libxt", "vawk;env=python2"],
            disk={"files": 2.0})]
    sv_batch_inputs = [["analysis"], ["genome_build"],
                       ["work_bam_plus", "disc"], ["work_bam_plus", "sr"],
                       ["config", "algorithm", "background", "cnv_reference"],
                       ["config", "algorithm", "tools_on"],
                       ["config", "algorithm", "tools_off"],
                       ["config", "algorithm", "svprioritize"],
                       ["config", "algorithm", "svvalidate"], ["regions", "sample_callable"],
                       ["genome_resources", "variation", "gc_profile"],
                       ["genome_resources", "variation", "germline_het_pon"],
                       ["genome_resources", "aliases", "snpeff"], ["reference", "snpeff", "genome_build"],
                       ["sv_coverage_rec"]]
    if checkpoints.get("vc"):
        sv_batch_inputs.append(["variants", "samples"])
    steps = [s("calculate_sv_bins", "multi-combined",
               [["align_bam"], ["reference", "fasta", "base"],
                ["metadata", "batch"], ["metadata", "phenotype"],
                ["config", "algorithm", "background", "cnv_reference"],
                ["config", "algorithm", "callable_regions"],
                ["config", "algorithm", "coverage_interval"],
                ["config", "algorithm", "exclude_regions"],
                ["config", "algorithm", "sv_regions"],
                ["config", "algorithm", "variant_regions"],
                ["config", "algorithm", "variant_regions_merged"],
                ["config", "algorithm", "seq2c_bed_ready"],
                ["config", "algorithm", "svcaller"],
                ["depth", "variant_regions", "regions"],
                ["genome_resources", "variation", "lcr"], ["genome_resources", "variation", "polyx"],
                ["genome_resources", "variation", "encode_blacklist"],
                ["genome_resources", "rnaseq", "gene_bed"]],
               [cwlout("sv_bin_rec", "record",
                       fields=[cwlout(["regions", "bins", "target"], ["File", "null"]),
                               cwlout(["regions", "bins", "antitarget"], ["File", "null"]),
                               cwlout(["regions", "bins", "gcannotated"], ["File", "null"]),
                               cwlout(["regions", "bins", "group"], ["string", "null"]),
                               cwlout("inherit")])],
               "bcbio-vc", ["bedtools", "cnvkit"],
               disk={"files": 1.5}, cores=1),
             s("calculate_sv_coverage", "multi-parallel",
               [["sv_bin_rec"]],
               [cwlout("sv_rawcoverage_rec", "record",
                       fields=[cwlout(["depth", "bins", "target"], ["File", "null"]),
                               cwlout(["depth", "bins", "antitarget"], ["File", "null"]),
                               cwlout(["depth", "bins", "seq2c"], ["File", "null"]),
                               cwlout("inherit")])],
               "bcbio-vc", ["mosdepth", "cnvkit", "seq2c"],
               disk={"files": 1.5}),
             s("normalize_sv_coverage", "multi-combined",
               [["sv_rawcoverage_rec"]],
               [cwlout("sv_coverage_rec", "record",
                       fields=[cwlout(["depth", "bins", "normalized"], ["File", "null"]),
                               cwlout(["depth", "bins", "background"], ["File", "null"]),
                               cwlout("inherit")])],
               "bcbio-vc", ["cnvkit"],
               disk={"files": 1.5}),
             s("batch_for_sv", "multi-batch", sv_batch_inputs,
               [cwlout("sv_batch_rec", "record")],
               "bcbio-vc",
               unlist=[["config", "algorithm", "svcaller"]]),
             w("svcall", "multi-parallel", sv, []),
             s("summarize_sv", "multi-combined",
               [["sv_rec"]],
               [cwlout(["sv", "calls"], {"type": "array", "items": ["File", "null"]}),
                cwlout(["sv", "supplemental"], {"type": "array", "items": ["File"]}),
                cwlout(["sv", "prioritize", "tsv"], {"type": "array", "items": ["File", "null"]}),
                cwlout(["sv", "prioritize", "raw"], {"type": "array", "items": ["File", "null"]}),
                cwlout(["svvalidate", "grading_summary"], ["File", "null"]),
                cwlout(["svvalidate", "grading_plots"], {"type": "array", "items": ["File", "null"]})],
               "bcbio-vc", ["bcbio-prioritize"], disk={"files": 1.0}, cores=1)]
    final_outputs = [["sv", "calls"], ["svvalidate", "grading_summary"], ["sv", "prioritize", "tsv"],
                     ["sv", "prioritize", "raw"], ["sv", "supplemental"]]
    return steps, final_outputs
1125408
def rsem(job, job_vars):
    
    input_args, ids = job_vars
    work_dir = job.fileStore.getLocalTempDir()
    cpus = input_args['cpu_count']
    sudo = input_args['sudo']
    single_end_reads = input_args['single_end_reads']
    # I/O
    filtered_bam, rsem_ref = return_input_paths(job, work_dir, ids, 'filtered.bam', 'rsem_ref.zip')
    subprocess.check_call(['unzip', '-o', os.path.join(work_dir, 'rsem_ref.zip'), '-d', work_dir])
    output_prefix = 'rsem'
    # Make tool call to Docker
    parameters = ['--quiet',
                  '--no-qualities',
                  '-p', str(cpus),
                  '--forward-prob', '0.5',
                  '--seed-length', '25',
                  '--fragment-length-mean', '-1.0',
                  '--bam', docker_path(filtered_bam)]
    if not single_end_reads:
        parameters.extend(['--paired-end'])
    parameters.extend(['/data/rsem_ref/hg19_M_rCRS_ref', output_prefix])

    docker_call(tool='quay.io/ucsc_cgl/rsem:1.2.25--4e8d1b31d4028f464b3409c6558fb9dfcad73f88',
                tool_parameters=parameters, work_dir=work_dir, sudo=sudo)
    os.rename(os.path.join(work_dir, output_prefix + '.genes.results'), os.path.join(work_dir, 'rsem_gene.tab'))
    os.rename(os.path.join(work_dir, output_prefix + '.isoforms.results'), os.path.join(work_dir, 'rsem_isoform.tab'))
    # Write to FileStore
    ids['rsem_gene.tab'] = job.fileStore.writeGlobalFile(os.path.join(work_dir, 'rsem_gene.tab'))
    ids['rsem_isoform.tab'] = job.fileStore.writeGlobalFile(os.path.join(work_dir, 'rsem_isoform.tab'))
    # Run child jobs
    return job.addChildJobFn(rsem_postprocess, job_vars).rv()
416554
def sort_bam_by_reference(job, job_vars):
    
    # Unpack variables
    input_args, ids = job_vars
    work_dir = job.fileStore.getLocalTempDir()
    # I/O
    sorted_bam, sorted_bai = return_input_paths(job, work_dir, ids, 'sorted.bam', 'sorted.bam.bai')
    output = os.path.join(work_dir, 'sort_by_ref.bam')
    # Call: Samtools
    ref_seqs = []
    handle = subprocess.Popen(["samtools", "view", "-H", sorted_bam], stdout=subprocess.PIPE).stdout
    for line in handle:
        if line.startswith("@SQ"):
            tmp = line.split("\t")
            chrom = tmp[1].split(":")[1]
            ref_seqs.append(chrom)
    handle.close()
    # Iterate through chromosomes to create mini-bams
    for chrom in ref_seqs:
        # job.addChildJobFn(sbbr_child, chrom, os.path.join(work_dir, chrom), sorted_bam)
        cmd_view = ["samtools", "view", "-b", sorted_bam, chrom]
        cmd_sort = ["samtools", "sort", "-m", "3000000000", "-n", "-", os.path.join(work_dir, chrom)]
        p1 = subprocess.Popen(cmd_view, stdout=subprocess.PIPE)
        subprocess.check_call(cmd_sort, stdin=p1.stdout)
    sorted_files = [os.path.join(work_dir, chrom) + '.bam' for chrom in ref_seqs]
    cmd = ["samtools", "cat", "-o", output] + sorted_files
    subprocess.check_call(cmd)
    # Write to FileStore
    ids['sort_by_ref.bam'] = job.fileStore.writeGlobalFile(output)
    rsem_id = job.addChildJobFn(transcriptome, job_vars, disk='30 G', memory='30 G').rv()
    exon_id = job.addChildJobFn(exon_count, job_vars, disk='30 G').rv()
    return exon_id, rsem_id
707707
def map_to_precursors_on_fly(seqs, names, loci, args):
    """"""
    precursor = precursor_sequence(loci, args.ref).upper()
    dat = dict()
    for s, n in itertools.izip(seqs, names):
        res = pyMatch.Match(precursor, str(s), 1, 3)
        if res > -1:
            dat[n] = [res, res + len(s)]
    logger.debug("mapped in %s: %s out of %s" % (loci, len(dat), len(seqs)))
    return dat