Using the data from the BLAST analyses, set the targets folder, and create the 'mapping file'. This is the
        genera-specific FASTA file that will be used for all the reference mapping; it replaces the 'bait file' in the
        code
496133
def targets(self):
        
        logging.info('Performing analysis with {} targets folder'.format(self.analysistype))
        for sample in self.runmetadata:
            if sample.general.bestassemblyfile != 'NA':
                sample[self.analysistype].targetpath = \
                    os.path.join(self.targetpath, 'genera', sample[self.analysistype].genus, '')
                # There is a relatively strict databasing scheme necessary for the custom targets. Eventually,
                # there will be a helper script to combine individual files into a properly formatted combined file
                try:
                    sample[self.analysistype].mappingfile = glob('{}*.fa'
                                                                 .format(sample[self.analysistype].targetpath))[0]
                # If the fasta file is missing, raise a custom error
                except IndexError as e:
                    # noinspection PyPropertyAccess
                    e.args = ['Cannot find the combined fasta file in {}. Please note that the file must have a '
                              '.fasta extension'.format(sample[self.analysistype].targetpath)]
                    if os.path.isdir(sample[self.analysistype].targetpath):
                        raise
                    else:
                        sample.general.bestassemblyfile = 'NA'
1008830
def build_blast_db_data
      # Create necessary directories and return.
      if @empty
        make_directories
        return
      end

      if @dir.blank?
        raise "DIR must be set to continue. Execute `rake -D` for more information."
      end

      unless VALID_TYPES.include?(@type)
        raise "Unknown type: #{@type}. " <<
          "Please provide one: both, nucl or prot."
      end

      check_dependencies

      make_directories

      begin
        @dir.split(':').each do |d|
          unless File.directory?(d)
            raise "Directory not found: #{d}"
          end

          @data = Dir.glob("#{d}/*.{tgz,tar.gz,tbz,tar.bz2}")

          if @data.blank?
            raise "Data not found. Please check your directory and try " <<
            "again.\nDirectory Entered: #{d}"
          end

          @data.each do |s|
            if s =~ GZIP
              files = `tar -tzf #{s}`
              flag  = "z"
            elsif s =~ BZIP
              files = `tar -tjf #{s}`
              flag  = "j"
            end
            files = files.split(/\n/)
            files.each do |f|
              if f.include?(@prot_file)
                file_name = create_file_name(f, @blastdb_dir)
                extract_files(s, f, flag, File.join(@blastdb_dir, file_name, "peptides.fa"))
              elsif f.include?(@nucl_file)
                file_name = create_file_name(f, @blastdb_dir)
                extract_files(s, f, flag, File.join(@blastdb_dir, file_name, "contigs.fa"))
              elsif f.include?("gff")
                file_name = create_file_name(f, @gff_dir)
                extract_files(s, f, flag, File.join(@gff_dir, file_name, "annots.gff"))
              end
            end
          end
        end
        build_blast_db(@blastdb_dir)
      rescue Exception => e
        # Remove empty directories.
        `rm -rf #{@blastdb_dir}/*` if File.directory?(@blastdb_dir)
        `rm -rf #{@gff_dir}/*` if File.directory?(@gff_dir)
        raise e
      end
      readme
    end
1151992
def _prep_genome(out_dir, data):
    
    genome_name = utils.splitext_plus(os.path.basename(dd.get_ref_file(data)))[0]
    out_dir = utils.safe_makedir(os.path.join(out_dir, genome_name))
    ref_file = dd.get_ref_file(data)
    utils.symlink_plus(ref_file, os.path.join(out_dir, os.path.basename(ref_file)))
    with open(os.path.join(out_dir, "GenomeSize.xml"), "w") as out_handle:
        out_handle.write('<sequenceSizes genomeName="%s">' % genome_name)
        for c in pysam.AlignmentFile("%s.dict" % utils.splitext_plus(ref_file)[0]).header["SQ"]:
            cur_ploidy = ploidy.get_ploidy([data], region=[c["SN"]])
            out_handle.write('<chromosome fileName="%s" contigName="%s" totalBases="%s" knownBases="%s" '
                             'isCircular="false" ploidy="%s" md5="%s"/>' %
                             (os.path.basename(ref_file), c["SN"], c["LN"], c["LN"], cur_ploidy, c["M5"]))
        out_handle.write('</sequenceSizes>')
    return out_dir
605244
def add_new_target(self, address, target_type, target_base=None, dependencies=None,
                     derived_from=None, **kwargs):
    
    rel_target_base = target_base or address.spec_path
    abs_target_base = os.path.join(get_buildroot(), rel_target_base)
    if not os.path.exists(abs_target_base):
      os.makedirs(abs_target_base)
      # TODO: Adding source roots on the fly like this is yucky, but hopefully this
      # method will go away entirely under the new engine. It's primarily used for injecting
      # synthetic codegen targets, and that isn't how codegen will work in the future.
    if not self.source_roots.find_by_path(rel_target_base):
      # TODO: Set the lang and root category (source/test/thirdparty) based on the target type?
      self.source_roots.add_source_root(rel_target_base)
    if dependencies:
      dependencies = [dep.address for dep in dependencies]

    self.build_graph.inject_synthetic_target(address=address,
                                             target_type=target_type,
                                             dependencies=dependencies,
                                             derived_from=derived_from,
                                             **kwargs)
    new_target = self.build_graph.get_target(address)

    return new_target
1174595
def _varscan_work(align_bams, ref_file, items, target_regions, out_file):
    
    config = items[0]["config"]

    orig_out_file = out_file
    out_file = orig_out_file.replace(".vcf.gz", ".vcf")

    max_read_depth = "1000"
    sample_list = _create_sample_list(align_bams, out_file)
    mpileup = samtools.prep_mpileup(align_bams, ref_file, config, max_read_depth,
                                    target_regions=target_regions, want_bcf=False)
    # VarScan fails to generate a header on files that start with
    # zerocoverage calls; strip these with grep, we're not going to
    # call on them
    remove_zerocoverage = r"{ ifne grep -v -P '\t0\t\t$' || true; }"
    # we use ifne from moreutils to ensure we process only on files with input, skipping otherwise
    # http://manpages.ubuntu.com/manpages/natty/man1/ifne.1.html
    with tx_tmpdir(items[0]) as tmp_dir:
        jvm_opts = _get_jvm_opts(config, tmp_dir)
        opts = " ".join(_varscan_options_from_config(config))
        min_af = float(utils.get_in(config, ("algorithm", "min_allele_fraction"), 10)) / 100.0
        fix_ambig_ref = vcfutils.fix_ambiguous_cl()
        fix_ambig_alt = vcfutils.fix_ambiguous_cl(5)
        py_cl = os.path.join(os.path.dirname(sys.executable), "py")
        export = utils.local_path_export()
        cmd = ("{export} {mpileup} | {remove_zerocoverage} | "
               "ifne varscan {jvm_opts} mpileup2cns {opts} "
               "--vcf-sample-list {sample_list} --min-var-freq {min_af} --output-vcf --variants | "
               """{py_cl} -x 'bcbio.variation.vcfutils.add_contig_to_header(x, "{ref_file}")' | """
               "{py_cl} -x 'bcbio.variation.varscan.fix_varscan_output(x)' | "
               "{fix_ambig_ref} | {fix_ambig_alt} | ifne vcfuniqalleles > {out_file}")
        do.run(cmd.format(**locals()), "Varscan", None,
                [do.file_exists(out_file)])
    os.remove(sample_list)
    # VarScan can create completely empty files in regions without
    # variants, so we create a correctly formatted empty file
    if os.path.getsize(out_file) == 0:
        write_empty_vcf(out_file)

    if orig_out_file.endswith(".gz"):
        vcfutils.bgzip_and_index(out_file, config)
1188571
def _copy_sources(self, target, results_dir):
    
    buildroot = get_buildroot()
    source_relative_to = target.address.spec_path
    for source in target.sources_relative_to_buildroot():
      dest = os.path.join(results_dir, os.path.relpath(source, source_relative_to))
      safe_mkdir(os.path.dirname(dest))
      shutil.copyfile(os.path.join(buildroot, source), dest)
1365330
def x_build_targets_target( self, node ):
        '''
        
        '''
        target_node = node
        name = self.get_child_data(target_node,tag='name',strip=True)
        path = self.get_child_data(target_node,tag='path',strip=True)
        jam_target = self.get_child_data(target_node,tag='jam-target',strip=True)
        #~ Map for jam targets to virtual targets.
        self.target[jam_target] = {
            'name' : name,
            'path' : path
            }
        #~ Create the ancestry.
        dep_node = self.get_child(self.get_child(target_node,tag='dependencies'),tag='dependency')
        while dep_node:
            child = self.get_data(dep_node,strip=True)
            child_jam_target = '<p%s>%s' % (path,child.split('//',1)[1])
            self.parent[child_jam_target] = jam_target
            dep_node = self.get_sibling(dep_node.nextSibling,tag='dependency')
        return None
1252421
def loadTargetsFromFile(self, filename, cols = None, everyNrows = 1,
                            delim = ' ', checkEven = 1):
        
        self.targets = self.loadVectors(filename, cols, everyNrows,
                                        delim, checkEven)
345429
def write_prepareref_fasta_file(outfile, gene_coords, genes_need_upstream, genes_non_upstream, upstream_before=100, upstream_after=100):
    ''''''
    tmp_dict = {}
    fasta_in = os.path.join(data_dir, 'NC_000962.3.fa.gz')
    pyfastaq.tasks.file_to_dict(fasta_in, tmp_dict)
    ref_seq = tmp_dict['NC_000962.3']

    with open(outfile, 'w') as f:
        for gene in genes_non_upstream:
            start = gene_coords[gene]['start']
            end = gene_coords[gene]['end']
            if start < end:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[start:end+1])
            else:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[end:start+1])
                gene_fa.revcomp()

            print(gene_fa, file=f)

        for gene in genes_need_upstream:
            start = gene_coords[gene]['start']
            end = gene_coords[gene]['end']
            if start < end:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[start - upstream_before:start + upstream_after])
            else:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[start - upstream_after + 1:start + upstream_before + 1])
                gene_fa.revcomp()

            gene_fa.id += '_upstream'
            print(gene_fa, file=f)
286621
def populateFromFile(self, dataUrl):
        
        self._dataUrl = dataUrl
        fastaFile = self.getFastaFile()
        for referenceName in fastaFile.references:
            reference = HtslibReference(self, referenceName)
            # TODO break this up into chunks and calculate the MD5
            # in bits (say, 64K chunks?)
            bases = fastaFile.fetch(referenceName)
            md5checksum = hashlib.md5(bases).hexdigest()
            reference.setMd5checksum(md5checksum)
            reference.setLength(len(bases))
            self.addReference(reference)
1215622
def blastnthreads(self):
        """"""
        # Create the threads for the BLAST analysis
        for i in range(self.cpus):
            threads = Thread(target=self.runblast, args=())
            threads.setDaemon(True)
            threads.start()
        # Populate threads for each gene, genome combination
        for sample in self.metadata:
            make_path(sample[self.analysistype].reportdir)
            sample[self.analysistype].report = os.path.join(
                sample[self.analysistype].reportdir, '{}.csv'.format(sample.name))
            if sample[self.analysistype].combinedtargets != 'NA':
                # Add each fasta file combination to the threads
                self.blastqueue.put((sample.general.bestassemblyfile, sample[self.analysistype].combinedtargets,
                                     sample))
        # Join the threads
        self.blastqueue.join()
1023486
def _prepare_sample(data, run_folder):
    
    want = set(["description", "files", "genome_build", "name", "analysis", "upload", "algorithm"])
    out = {}
    for k, v in data.items():
        if k in want:
            out[k] = _relative_paths(v, run_folder)
    if "algorithm" not in out:
        analysis, algorithm = _select_default_algorithm(out.get("analysis"))
        out["algorithm"] = algorithm
        out["analysis"] = analysis
    description = "%s-%s" % (out["name"], clean_name(out["description"]))
    out["name"] = [out["name"], description]
    out["description"] = description
    return out
674996
def build_blast_db(blastdb)
      Dir.glob(File.expand_path(blastdb) + "/*").each do |d|
        if File.directory?(d)
          contigs  = File.join(d, "contigs.fa")
          peptides = File.join(d, "peptides.fa")

          found = false

          if File.exists?(contigs) && File.readable?(contigs)
            execute_makeblastdb("nucl", d, contigs)
            found = true
          end
          if File.exists?(peptides) && File.readable?(peptides)
            execute_makeblastdb("prot", d, peptides)
            found = true
          end

          unless found
            raise "Extracted data not found for #{contigs} or #{peptides}. " <<
            "Make sure you supplied the correct data directory and file names."
          end
        end
      end
    end
1031388
def prepare_intervals(data, region_file, work_dir):
    
    target_file = os.path.join(work_dir, "%s-target.interval_list" % dd.get_sample_name(data))
    if not utils.file_uptodate(target_file, region_file):
        with file_transaction(data, target_file) as tx_out_file:
            params = ["-T", "PreprocessIntervals", "-R", dd.get_ref_file(data),
                      "--interval-merging-rule", "OVERLAPPING_ONLY",
                      "-O", tx_out_file]
            if dd.get_coverage_interval(data) == "genome":
                params += ["--bin-length", "1000", "--padding", "0"]
            else:
                params += ["-L", region_file, "--bin-length", "0", "--padding", "250"]
            _run_with_memory_scaling(params, tx_out_file, data)
    return target_file
1179713
def mapper(self, sample):
        
        if sample.general.bestassemblyfile != "NA":
            # Define the Qualimap log and report files
            reportfile = os.path.join(sample.general.QualimapResults, 'genome_results.txt')
            # Define the Qualimap call
            qualimapcall = 'qualimap bamqc -bam {} -outdir {}'.format(sample.general.sortedbam,
                                                                      sample.general.QualimapResults)
            sample.commands.qualimap = qualimapcall
            # Initialise a dictionary to hold the Qualimap results
            qdict = dict()
            # If the report file doesn't exist, run Qualimap, and print logs to the log file
            if not os.path.isfile(reportfile):
                tlock = threading.Lock()
                out, err = run_subprocess(sample.commands.qualimap)
                tlock.acquire()
                write_to_logfile(sample.commands.qualimap, sample.commands.qualimap, self.logfile,
                                 sample.general.logout, sample.general.logerr, None, None)
                write_to_logfile(out, err, self.logfile, sample.general.logout, sample.general.logerr, None, None)
                tlock.release()
            # Initialise a genobject to store the coverage dictionaries
            sample.depth = GenObject()
            sample.depth.length = dict()
            sample.depth.bases = dict()
            sample.depth.coverage = dict()
            sample.depth.stddev = dict()
            try:
                with open(reportfile) as report:
                    # Read the report
                    for line in report:
                        # Sanitise the keys and values using self.analyze
                        key, value = self.analyze(line)
                        # If the keys and values exist, enter them into the dictionary
                        if (key, value) != (None, None):
                            qdict[key] = value
                        if 'Coverage per contig' in line:
                            for contigline in report:
                                try:
                                    _, name, length, bases, coverage, stddev = contigline.rstrip().split('\t')
                                    sample.depth.length.update({name: length})
                                    sample.depth.bases.update({name: bases})
                                    sample.depth.coverage.update({name: coverage})
                                    sample.depth.stddev.update({name: stddev})
                                except ValueError:
                                    pass

            except (IOError, FileNotFoundError):
                pass
            # If there are values in the dictionary
            if qdict:
                # Make new category for Qualimap results and populate this category with the report data
                for attribute in qdict:
                    # Remove the 'X' from the depth values e.g. 40.238X
                    setattr(sample.mapping, attribute, qdict[attribute].rstrip('X'))
532517
def filer(args):
        
        # List to store all the metadata objects
        samples = list()
        # Find all the sequence files in the path
        fastas = sorted(glob(os.path.join(args.sequencepath, '*.fa*')))
        for fasta in fastas:
            # Create a metadata object for each sample
            metadata = MetadataObject()
            # Populate the metadata object with the required attributes
            metadata.name = os.path.splitext(os.path.basename(fasta))[0]
            metadata.general = GenObject()
            metadata.commands = GenObject()
            metadata.general.bestassemblyfile = fasta
            metadata.general.outputdirectory = os.path.join(args.sequencepath, metadata.name)
            samples.append(metadata)
        return samples
1359448
def alleleupdater1(self, sample, gene, targetallele):
        
        from io import StringIO
        from Bio.Blast import NCBIXML
        from Bio.SeqRecord import SeqRecord
        from Bio.Seq import Seq
        from Bio.Alphabet import generic_dna
        # Find the name of the allele file based on whether the gene name is found in the file name
        genefile = [x for x in sample[self.analysistype].alleles if gene in x][0]
        # Remove the extension of the gene file for use in makeblastdb and blastn
        genefilenoext = genefile.split(".")[0]
        # Create the blast databases every time the method is called
        subprocess.call(shlex.split('makeblastdb -in {} -parse_seqids -max_file_sz 2GB -dbtype nucl -out {}'
                                    .format(genefile, genefilenoext)), stdout=self.fnull, stderr=self.fnull)
        # Re-perform BLAST analyses, but using an XML output that stores the alignment
        blastn = NcbiblastnCommandline(query=sample.general.bestassemblyfile, db=genefilenoext, evalue=0.1, outfmt=5)
        # Note that there is no output file specified -  the search results are currently stored in stdout
        stdout, stderr = blastn()
        if stdout.find('Hsp') != -1:
            # Map the blast record to memory
            blast_handle = StringIO(stdout)
            # Open record from memory-mapped file
            records = NCBIXML.parse(blast_handle)
            # Iterate through the records in the blast results
            for record in records:  # This process is just to retrieve HSPs from xml files
                for alignment in record.alignments:
                    for hsp in alignment.hsps:
                        # Extract the allele name from the blast result
                        allele = str(alignment.accession.split("_")[-1]) if "_" in alignment.accession else \
                            str(alignment.accession.split("-")[-1])
                        # The point of this blast is to extract the sequence of the allele
                        # If the allele name matches the target allele (the closest match e.g. BACT000063_20 -> 20)
                        if allele == str(targetallele):
                            # As there is some discrepancy with the capitalisation of terms, make sure it is consistent
                            analysistype = 'rMLST' if self.analysistype.lower() == 'rmlst' else 'MLST'
                            # Set the directory containing the profile and alleles
                            alleledir = self.referencefilepath + analysistype if self.pipeline else \
                                self.referencefilepath
                            # The name of the supplemental allele file (without and with the .fa extension)
                            allelefilenoext = '{}/OLC_{}_alleles'.format(alleledir, analysistype) if self.pipeline \
                                else '{}/rMLST_combined'.format(alleledir)
                            allelefile = allelefilenoext + '.fa' if self.pipeline else allelefilenoext + '.fasta'
                            # Create the file if it doesn't exist
                            open(allelefile, 'a').close()
                            # Create a list of all the blast database files in the folder
                            dbfiles = glob('{}.n*'.format(allelefilenoext))
                            # Remove the database files
                            map(lambda y: os.remove(y), dbfiles)
                            # Create the necessary blast database files
                            subprocess.call(
                                shlex.split('makeblastdb -in {} -parse_seqids -max_file_sz 2GB -dbtype nucl -out {}'
                                            .format(allelefile, allelefilenoext)), stdout=self.fnull, stderr=self.fnull)

                            # Perform BLAST analysis using the supplemental allele file as the database
                            nestedblastn = NcbiblastnCommandline(db=allelefilenoext, evalue=0.1, outfmt=5)
                            # There is no output file specified; the search results are currently stored in stdout
                            # Additionally, the sequence from the previous BLAST query is used as stdin in this BLAST
                            import Bio.Application
                            try:
                                nestedstdout, nestedstderr = nestedblastn(stdin=hsp.query)
                            except Bio.Application.ApplicationError:
                                nestedstdout = ''
                                pass
                            # Search stdout for matches - if the term Hsp appears (the .find function will NOT
                            # return -1), a match has been found, and stdout is written to file
                            if nestedstdout.find('Hsp') != -1:
                                nested_blast_handle = StringIO(nestedstdout)
                                # Open record from memory-mapped file
                                nestedrecords = NCBIXML.parse(nested_blast_handle)
                                # Initialise variables
                                previousallele = ''
                                bitscore = ''
                                # Iterate through the records
                                for nestedrecord in nestedrecords:
                                    for nestedalignment in nestedrecord.alignments:
                                        for nestedhsp in nestedalignment.hsps:
                                            # Calculate the percent identity
                                            percentidentity = float("%.2f" % float(float(nestedhsp.identities) /
                                                                                   float(nestedalignment.length) * 100))
                                            # Only set the previous allele and the bitscore if the % identity is 100%
                                            if percentidentity == 100:
                                                # Get the allele number
                                                previousallele = nestedalignment.accession.split("_")[-1]
                                                bitscore = nestedhsp.score

                                # If this allele has already been found in the supplemental file, return these data
                                if previousallele:
                                    return gene, previousallele, 100.0, bitscore
                                # Otherwise, get the allele sequence into the supplemental file
                                else:
                                    # Initialise variables
                                    allelelist = []
                                    allelenumber = 1000000
                                    # Open the allele file to append
                                    with open(allelefile, 'ab+') as supplemental:
                                        for line in supplemental:
                                            # As there are multiple genes (e.g. BACT000001, BACT000063, etc.) check to
                                            # see if the gene of interest is in the line
                                            if gene in line:
                                                # Append the header to the list
                                                allelelist.append(line)
                                        # Find the last allele number from the header
                                        try:
                                            allelenumber = int(allelelist[-1].split("_")[1]) + 1
                                        # If there are no alleles for the gene of interest, then pass
                                        except IndexError:
                                            pass
                                        # Puts the HSP in the correct order -  hits to the negative strand will be
                                        # reversed compared to what we're looking for
                                        allelesequence = Seq(hsp.query, generic_dna)
                                        if hsp.sbjct_start < hsp.sbjct_end:
                                            end = hsp.sbjct_end
                                        else:
                                            end = hsp.sbjct_start
                                            allelesequence = allelesequence.reverse_complement()
                                        # Screen out hits that are shorter than the targets
                                        # Keeping this format though this statement could be re-written more efficiently
                                        if end < alignment.length:
                                            pass
                                        # The header will be >BACT00001_1000000
                                        definitionline = '{}_{}'.format(gene, allelenumber)
                                        # Create a sequence record using BioPython
                                        fasta = SeqRecord(allelesequence,
                                                          # Without this, the header will be improperly formatted
                                                          description='',
                                                          # Use >:definitionline as the header
                                                          id=definitionline)
                                        # Use the SeqIO module to properly format the new sequence record
                                        SeqIO.write(fasta, supplemental, "fasta")
                                        # Return the necessary values
                                        return gene, allelenumber, 100.0, hsp.score
                            # If there are hits in the supplemental allele file
                            else:
                                # Initialise variables
                                allelelist = []
                                allelenumber = 1000000
                                # Open the allele file to find the last allele associated with the gene of interest
                                with open(allelefile, 'a+') as supplemental:
                                    for line in supplemental:
                                        if gene in line:
                                            allelelist.append(line)
                                    try:
                                        allelenumber = int(allelelist[-1].split("_")[1]) + 1
                                    except IndexError:
                                        pass
                                    # Puts the HSP in the correct order -  hits to the negative strand will be
                                    # reversed compared to what we're looking for
                                    allelesequence = Seq(hsp.query, generic_dna)
                                    if hsp.sbjct_start < hsp.sbjct_end:
                                        end = hsp.sbjct_end
                                    else:
                                        end = hsp.sbjct_start
                                        allelesequence = allelesequence.reverse_complement()
                                    # Screen out hits that are shorter than the targets
                                    # Keeping this format, though this if statement could be re-written more efficiently
                                    if end < alignment.length:
                                        pass
                                    definitionline = '{}_{}'.format(gene, allelenumber)
                                    # Create a sequence record using BioPython
                                    fasta = SeqRecord(allelesequence,
                                                      # If this is not added, the header will not be formatted properly
                                                      description='',
                                                      # Use >:definitionline as the header
                                                      id=definitionline)
                                    # Use the SeqIO module to properly format the new sequence record
                                    SeqIO.write(fasta, supplemental, "fasta")
                                    # Return the appropriate information
                                    return gene, allelenumber, 100.0, hsp.score