DMDID and ADMID relationships found in the METS file.
183494
private void setRDFAsDatastream(StringBuffer buf) {

        DatastreamXMLMetadata ds = new DatastreamXMLMetadata();
        // set the attrs common to all datastream versions
        ds.DatastreamID = "RELS-INT";
        ds.DSVersionable = false;
        ds.DSFormatURI = m_dsFormatURI;
        ds.DatastreamAltIDs = m_dsAltIDs;
        ds.DSVersionID = "RELS-INT.0";
        ds.DSLabel =
                "DO NOT EDIT: System-generated datastream to preserve METS DMDID/ADMID relationships.";
        ds.DSCreateDT = new Date();
        ds.DSMIME = "application/rdf+xml";
        // set the attrs specific to datastream version
        ds.DSControlGrp = "X";
        ds.DSState = "A";
        ds.DSLocation =
                m_obj.getPid() + "+" + ds.DatastreamID + "+" + ds.DSVersionID;
        ds.DSLocationType = Datastream.DS_LOCATION_TYPE_INTERNAL;
        ds.DSInfoType = "DATA";
        ds.DSMDClass = DatastreamXMLMetadata.TECHNICAL;

        // now set the xml content stream itself...
        try {
            ds.xmlContent = buf.toString().getBytes(m_characterEncoding);
            ds.DSSize = ds.xmlContent.length;
        } catch (UnsupportedEncodingException uee) {
            logger.error("Encoding error when creating RELS-INT datastream", uee);
        }
        // FINALLY! add the RDF and an inline xml datastream in the digital object
        m_obj.addDatastreamVersion(ds, true);
    }
1020334
def _adj_for_assc(self):
        """"""
        if self.gene2gos:
            gos_assoc = set(get_b2aset(self.gene2gos).keys())
            if 'item_marks' not in self.kws:
                self.kws['item_marks'] = {go:'>' for go in gos_assoc}
            if 'include_only' not in self.kws:
                gosubdag = GoSubDag(gos_assoc, self.gosubdag.go2obj,
                                    self.gosubdag.relationships)
                self.kws['include_only'] = gosubdag.go2obj
125799
def match_RCSB_pdb_chains(pdb_id1, pdb_id2, cut_off = 60.0, allow_multiple_matches = False, multiple_match_error_margin = 3.0, use_seqres_sequences_if_possible = True, strict = True):
    ''''''
    try:
        stage = pdb_id1
        pdb_1 = PDB(retrieve_pdb(pdb_id1), strict = strict)
        stage = pdb_id2
        pdb_2 = PDB(retrieve_pdb(pdb_id2), strict = strict)
    except (PDBParsingException, NonCanonicalResidueException, PDBValidationException), e:
        raise PDBParsingException("An error occurred while loading %s: '%s'" % (stage, str(e)))

    return match_pdb_chains(pdb_1, pdb_id1, pdb_2, pdb_id2, cut_off = cut_off, allow_multiple_matches = allow_multiple_matches, multiple_match_error_margin = multiple_match_error_margin, use_seqres_sequences_if_possible = use_seqres_sequences_if_possible)
250662
private List<List<Atom[]>> getMatchedAtomsOfLinkages(
			ModificationCondition condition, Map<Component, Set<Group>> mapCompGroups) {
		List<ModificationLinkage> linkages = condition.getLinkages();
		int nLink = linkages.size();

		List<List<Atom[]>> matchedAtomsOfLinkages =
				new ArrayList<List<Atom[]>>(nLink);

		for (int iLink=0; iLink<nLink; iLink++) {
			ModificationLinkage linkage = linkages.get(iLink);
			Component comp1 = linkage.getComponent1();
			Component comp2 = linkage.getComponent2();

//			boolean isAA1 = comp1.;
//			boolean isAA2 = comp2.getType()==true;

			Set<Group> groups1 = mapCompGroups.get(comp1);
			Set<Group> groups2 = mapCompGroups.get(comp2);

			List<Atom[]> list = new ArrayList<Atom[]>();

			List<String> potentialNamesOfAtomOnGroup1 = linkage.getPDBNameOfPotentialAtomsOnComponent1();
			for (String name : potentialNamesOfAtomOnGroup1) {
				if (name.equals("*")) {
					// wildcard
					potentialNamesOfAtomOnGroup1 = null; // search all atoms
					break;
				}
			}

			List<String> potentialNamesOfAtomOnGroup2 = linkage.getPDBNameOfPotentialAtomsOnComponent2();
			for (String name : potentialNamesOfAtomOnGroup2) {
				if (name.equals("*")) {
					// wildcard
					potentialNamesOfAtomOnGroup2 = null; // search all atoms
					break;
				}
			}

			for (Group g1 : groups1) {
				for (Group g2 : groups2) {
					if (g1.equals(g2)) {
						continue;
					}

					// only for wildcard match of two residues
					boolean ignoreNCLinkage =
						potentialNamesOfAtomOnGroup1 == null &&
						potentialNamesOfAtomOnGroup2 == null &&
						residues.contains(g1) &&
						residues.contains(g2);

					Atom[] atoms = StructureUtil.findNearestAtomLinkage(
							g1, g2,
							potentialNamesOfAtomOnGroup1,
							potentialNamesOfAtomOnGroup2,
							ignoreNCLinkage,
							bondLengthTolerance);
					if (atoms!=null) {
						list.add(atoms);
					}
				}
			}

			if (list.isEmpty()) {
				// broken linkage
				break;
			}

			matchedAtomsOfLinkages.add(list);
		}

		return matchedAtomsOfLinkages;
	}
749780
def dannots2dalignbed2dannotsagg(cfg):
    
    datatmpd=cfg['datatmpd']
    
    daannotp=f'{datatmpd}/08_dannot.tsv'
    cfg['daannotp']=daannotp
    dannotsaggp=cfg['dannotsaggp']
    logging.info(basename(daannotp))
    if ((not exists(daannotp)) and (not exists(dannotsaggp))) or cfg['force']:
        gff_renamed_cols=[c+' annotation' if c in set(bed_colns).intersection(gff_colns) else c for c in gff_colns]
        dannots=pd.read_csv(cfg['annotationsbedp'],sep='\t',
                   names=bed_colns+gff_renamed_cols,
                           low_memory=False)
        dannots=del_Unnamed(dannots)

        dannots=dannots.set_index('id')
        dannots['annotations count']=1
        # separate ids from attribute columns
        dannots=lambda2cols(dannots,lambdaf=gffatributes2ids,
                    in_coln='attributes',
                to_colns=['gene name','gene id','transcript id','protein id','exon id'])
        dannots=dannots.drop(['attributes']+[c for c in gff_renamed_cols if 'annotation' in c],axis=1)
        logging.debug('or this step takes more time?')
        to_table(dannots,daannotp)
#         to_table_pqt(dannots,daannotp)
    else:
#         dannots=read_table_pqt(daannotp)
        dannots=read_table(daannotp)
        dannots=del_Unnamed(dannots)
        
    logging.info(basename(dannotsaggp))
    if not exists(dannotsaggp) or cfg['force']:
        if not 'dannots' in locals():
#             dannots=read_table_pqt(daannotp)
            dannots=pd.read_table(daannotp,low_memory=False)
        dannots=del_Unnamed(dannots)
        dannots=dannots.reset_index()
        logging.debug('aggregating the annotations')
        from rohan.dandage.io_sets import unique 
        cols2aggf={'annotations count':np.sum,
                  'type': unique,
                  'gene name': unique,
                  'gene id': unique,
                  'transcript id': unique,
                  'protein id': unique,
                  'exon id': unique}
        dannotsagg=dannots.groupby('id').agg(cols2aggf)
        dannotsagg['annotations count']=dannotsagg['annotations count']-1
        dannotsagg.loc[dannotsagg['annotations count']==0,'region']='intergenic'
        dannotsagg.loc[dannotsagg['annotations count']!=0,'region']='genic'
        logging.debug('end of the slowest step')            
        del dannots    
        dannotsagg=dannotsagg.reset_index()
#         to_table_pqt(dannotsagg,dannotsaggp)
        dannotsagg.to_csv(dannotsaggp,sep='\t')        
    return cfg
769342
function () {
      var doUniqueSlotIdsCheck = function (compound) {
        var duplicateSlotIds = [];
        compound.slots.forEach(function (slot) {
          var slotId = slot.slotId;

          function filterBySlotId (slot) {
            return slot.slotId === slotId;
          }

          var arrayByUniqueSlotIds = compound.slots.filter(filterBySlotId);
          if (arrayByUniqueSlotIds.length > 1 && duplicateSlotIds.indexOf(slotId) < 0) {
            duplicateSlotIds.push(slotId);
          }
        });
        duplicateSlotIds.forEach(function (slotId) {
          errors.push('RuleViolation identified: Artifact \'' + compound.artifactId +
            '\' contains multiple slots with \'slotId\'=\'' + slotId + '\'');
        });
      };
      //
      if (doc.artifacts && doc.artifacts.compoundComponents) {
        doc.artifacts.compoundComponents.forEach(function (compound) {
          if (compound.slots) {
            doUniqueSlotIdsCheck(compound);
          }
        });
      }
    }
42212
function dmd (templateData, options) {
  options = new DmdOptions(options)
  if (skipCache(options)) {
    return generate(templateData, options)
  } else {
    const cached = dmd.cache.readSync([ templateData, options, dmdVersion ])
    if (cached) {
      return cached
    } else {
      return generate(templateData, options)
    }
  }
}
67576
def geneseekrsummary(self):
        
        logging.info('Creating predicted plasmid-borne gene summary table')
        with open(os.path.join(self.reportpath, 'plasmid_borne_summary.csv'), 'w') as pbs:
            data = 'Strain,Gene,PercentIdentity,Contig,Location,PlasmidIncompatibilitySets\n'
            for sample in self.metadata:
                # Create a flag to determine whether the strain name needs to be added to the data string if there
                # were no results
                result_bool = False
                # Initialise a dictionary to store a set of all the incompatibility types listed for a contig.
                # As the inc type will only be located on one of possibly several contigs associated with a predicted
                # plasmid, it is nice to know details about the plasmid
                inc_dict = dict()
                # Iterate through all the MOB recon outputs to populate the incompatibility set
                for primarykey, results in sample[self.analysistype].report_dict.items():
                    try:
                        inc = results['cluster_id']
                        # Convert the rep_type field (predicted incompatibilities) into a more a consistent
                        # format - pandas will call empty fields 'nan', which is a float
                        rep = str(results['rep_type']).replace(',', ';') if str(results['rep_type']) != 'nan' else 'ND'
                        # Add the incompatibility to the set
                        try:
                            inc_dict[inc].add(rep)
                        except KeyError:
                            inc_dict[inc] = set()
                            inc_dict[inc].add(rep)
                    except KeyError:
                        pass
                for primarykey, results in sample[self.analysistype].report_dict.items():
                    try:
                        contig = results['contig_id'].split('|')[1]
                        # Unicycler gives contigs names such as: 3_length=187116_depth=1.60x_circular=true - test
                        # to see if the contig name looks unicycler-like, and set the name appropriately (in this
                        # case, it would be 3)
                        if contig.split('_')[1].startswith('length'):
                            contig = contig.split('_')[0]
                        for gene, result_dict in sample.geneseekr_results.sampledata.items():
                            if contig == result_dict['query_id']:
                                percent_identity = result_dict['PercentIdentity']
                                # Set up the output string if the percent identity of the match is greater than the
                                # cutoff
                                if float(result_dict['PercentIdentity']) >= self.cutoff:
                                    # As there was at least a single gene passing the threshold, set the boolean to True
                                    result_bool = True
                                    data += '{sn},'.format(sn=sample.name)
                                    data += '{gene},{pi},{contig},{cid},{inc}\n'\
                                        .format(gene=gene,
                                                pi=percent_identity,
                                                contig=contig,
                                                cid=results['cluster_id'],
                                                inc=';'.join(sorted(inc_dict[str(results['cluster_id'])])))
                    except KeyError:
                        pass
                # If there were no results associated with the strain, make the row the strain name only
                if not result_bool:
                    data += '{sn}\n'.format(sn=sample.name)
            # Write the string to the report
            pbs.write(data)
979688
function collectNodeIds(stmt) {
  var ids = {},
      stack = [],
      curr;

  var push = stack.push.bind(stack);

  push(stmt);
  while(stack.length) {
    curr = stack.pop();
    switch(curr.type) {
      case "node": ids[curr.id] = true; break;
      case "edge": _.each(curr.elems, push); break;
      case "subgraph": _.each(curr.stmts, push); break;
    }
  }

  return _.keys(ids);
}
1281098
def _write_pair_information(gsd_file, structure):
    
    pair_types = []
    pair_typeid = []
    pairs = []
    for ai in structure.atoms:
        for aj in ai.dihedral_partners:
            #make sure we don't double add
            if ai.idx > aj.idx:
                ps = '-'.join(sorted([ai.type, aj.type], key=natural_sort))
                if ps not in pair_types:
                    pair_types.append(ps)
                pair_typeid.append(pair_types.index(ps))
                pairs.append((ai.idx, aj.idx))
    gsd_file.pairs.types = pair_types
    gsd_file.pairs.typeid = pair_typeid
    gsd_file.pairs.group = pairs
    gsd_file.pairs.N = len(pairs)
984646
def df_homology_models(self):
        """"""

        # TODO: add definitions of column names
        # model_date	Date the model was created
        # difficulty	Difficulty level of the modeling run (easy, medium, hard)
        # top_template_pdb	Top template used to model the protein (PDB ID)
        # top_template_chain	Chain of the top template used to model the protein
        # c_score	Confidence score of the homology model from [-5,2]
        # tm_score	Structural similarity score to best template
        # tm_score_err	Standard error of TM score
        # rmsd	RMSD to best template
        # rmsd_err	Standard error of RMSD
        # top_bsite_site_num	Cluster which contained the consensus binding site
        # top_bsite_c_score	Confidence score of the consensus binding site prediction from [0,1]
        # top_bsite_cluster_size	Number of predictions within this cluster
        # top_bsite_binding_residues	Residue numbers of top ranked binding site prediction
        # top_bsite_binding_location_coords	Cartesian coordinates of top ranked binding site prediction
        # top_bsite_pdb_ligand	Top predicted ligand to bind
        # top_bsite_ligand_cluster_counts	Number of predictions for all ligands
        # top_bsite_algorithm	Algorithm used for top ranked binding site prediction
        # top_bsite_c_score_method	Confidence score of just the top algorithm's binding site prediction from [0,1]. >0.35 is somewhat reliable.
        # top_bsite_pdb_template_id	PDB ID of the template used to make the prediction
        # top_bsite_pdb_template_chain	Chain of the PDB which has the ligand
        # top_ec_ec_number	Top predicted EC number
        # top_ec_c_score	Confidence score of the consensus EC number prediction from [0,1]
        # top_ec_pdb_template_id	Top PDB ID used as template for EC number prediction
        # top_ec_pdb_template_chain	Chain of PDB ID used as template for EC number prediction
        # top_ec_rmsd	RMSD to best template
        # top_ec_seq_coverage	Sequence coverage percentage to best template
        # top_ec_seq_ident	Sequence identity percentage to best template
        # top_ec_tm_score	Structural similarity score to best template
        # top_go_bp_go_id	Predicted GO ID for biological process
        # top_go_bp_go_term	Predicted GO term for biological process
        # top_go_bp_c_score	Confidence score of the GO BP prediction from [0,1]
        # top_go_cc_go_id	Predicted GO ID for cellular compartment
        # top_go_cc_go_term	Predicted GO term for cellular compartment
        # top_go_cc_c_score	Confidence score of the GO CC prediction from [0,1]
        # top_go_mf_go_id	Predicted GO ID for molecular function
        # top_go_mf_go_term	Predicted GO term for molecular function
        # top_go_mf_c_score	Confidence score of the GO MF prediction from [0,1]

        itasser_pre_df = []

        df_cols = ['id', 'structure_file', 'model_date', 'difficulty',
                   'top_template_pdb', 'top_template_chain', 'c_score',
                   'tm_score', 'tm_score_err', 'rmsd', 'rmsd_err',
                   'top_bsite_site_num', 'top_bsite_c_score', 'top_bsite_cluster_size', 'top_bsite_algorithm',
                   'top_bsite_pdb_template_id', 'top_bsite_pdb_template_chain', 'top_bsite_pdb_ligand',
                   'top_bsite_binding_location_coords', 'top_bsite_c_score_method', 'top_bsite_binding_residues',
                   'top_bsite_ligand_cluster_counts',
                   'top_ec_pdb_template_id', 'top_ec_pdb_template_chain', 'top_ec_tm_score', 'top_ec_rmsd',
                   'top_ec_seq_ident', 'top_ec_seq_coverage', 'top_ec_c_score', 'top_ec_ec_number',
                   'top_ec_binding_residues',
                   'top_go_mf_go_id', 'top_go_mf_go_term', 'top_go_mf_c_score', 'top_go_bp_go_id', 'top_go_bp_go_term',
                   'top_go_bp_c_score', 'top_go_cc_go_id', 'top_go_cc_go_term', 'top_go_cc_c_score']

        for p in self.get_homology_models():
            # Summary dataframe
            new_homology_dict = p.get_dict(df_format=True, only_attributes=df_cols)
            itasser_pre_df.append(new_homology_dict)

        df = pd.DataFrame.from_records(itasser_pre_df, columns=df_cols).set_index('id')
        return ssbio.utils.clean_df(df)
1074404
def align_transcriptome(fastq_file, pair_file, ref_file, data):
    
    work_bam = dd.get_work_bam(data)
    base, ext = os.path.splitext(work_bam)
    out_file = base + ".transcriptome" + ext
    if utils.file_exists(out_file):
        data = dd.set_transcriptome_bam(data, out_file)
        return data
    # bwa mem needs phred+33 quality, so convert if it is Illumina
    if dd.get_quality_format(data).lower() == "illumina":
        logger.info("bwa mem does not support the phred+64 quality format, "
                    "converting %s and %s to phred+33.")
        fastq_file = fastq.groom(fastq_file, data, in_qual="fastq-illumina")
        if pair_file:
            pair_file = fastq.groom(pair_file, data, in_qual="fastq-illumina")
    bwa = config_utils.get_program("bwa", data["config"])
    gtf_file = dd.get_gtf_file(data)
    gtf_fasta = index_transcriptome(gtf_file, ref_file, data)
    args = " ".join(_bwa_args_from_config(data["config"]))
    num_cores = data["config"]["algorithm"].get("num_cores", 1)
    samtools = config_utils.get_program("samtools", data["config"])
    cmd = ("{bwa} mem {args} -a -t {num_cores} {gtf_fasta} {fastq_file} "
           "{pair_file} ")
    with file_transaction(data, out_file) as tx_out_file:
        message = "Aligning %s and %s to the transcriptome." % (fastq_file, pair_file)
        cmd += "| " + postalign.sam_to_sortbam_cl(data, tx_out_file, name_sort=True)
        do.run(cmd.format(**locals()), message)
    data = dd.set_transcriptome_bam(data, out_file)
    return data
764843
function () {
      var doUniqueArtifactIdsCheck = function (artifactsIds) {
        var duplicateArtifactIds = [];
        artifactsIds.forEach(function (artifactId) {
          var artifactIdUnderTest = artifactId;

          function filterById (id) {
            return id === artifactIdUnderTest;
          }

          var arrayByUniqueArtifactIds = artifactsIds.filter(filterById);
          if (arrayByUniqueArtifactIds.length > 1 && duplicateArtifactIds.indexOf(artifactIdUnderTest) < 0) {
            duplicateArtifactIds.push(artifactIdUnderTest);
          }
        });
        duplicateArtifactIds.forEach(function (artifactId) {
          errors.push(
            'RuleViolation identified: WebPackage contains multiple artifacts with \'artifactId\'=\'' +
            artifactId + '\'');
        });
      };
      //
      var collectArtifactIds = function (next) {
        var artifactsIds = [];
        if (doc.artifacts) {
          Object.keys(doc.artifacts).forEach(function (artifactType) {
            doc.artifacts[ artifactType ].forEach(function (artifact) {
              artifactsIds.push(artifact.artifactId);
            });
          });
        }
        next(artifactsIds);
      };
      collectArtifactIds(doUniqueArtifactIdsCheck);
    }
400458
def combine_express(samples, combined):
    """"""
    if not combined:
        return None
    to_combine = [dd.get_express_counts(x) for x in
                  dd.sample_data_iterator(samples) if dd.get_express_counts(x)]
    gtf_file = dd.get_gtf_file(samples[0][0])
    isoform_to_gene_file = os.path.join(os.path.dirname(combined), "isoform_to_gene.txt")
    isoform_to_gene_file = express.isoform_to_gene_name(
        gtf_file, isoform_to_gene_file, next(dd.sample_data_iterator(samples)))
    if len(to_combine) > 0:
        eff_counts_combined_file = os.path.splitext(combined)[0] + ".isoform.express_counts"
        eff_counts_combined = count.combine_count_files(to_combine, eff_counts_combined_file, ext=".counts")
        to_combine = [dd.get_express_tpm(x) for x in
                      dd.sample_data_iterator(samples) if dd.get_express_tpm(x)]
        tpm_counts_combined_file = os.path.splitext(combined)[0] + ".isoform.express_tpm"
        tpm_counts_combined = count.combine_count_files(to_combine, tpm_counts_combined_file)
        to_combine = [dd.get_express_fpkm(x) for x in dd.sample_data_iterator(samples)
                      if dd.get_express_fpkm(x)]
        fpkm_counts_combined_file = os.path.splitext(combined)[0] + ".isoform.express_fpkm"
        fpkm_counts_combined = count.combine_count_files(to_combine, fpkm_counts_combined_file, ext=".fpkm")
        return {'counts': eff_counts_combined, 'tpm': tpm_counts_combined,
                'fpkm': fpkm_counts_combined, 'isoform_to_gene': isoform_to_gene_file}
    return {}
598785
def prt_gene_aart(self, geneids, prt=sys.stdout):
        """"""
        patgene = self.datobj.kws["fmtgene"]
        itemid2name = self.datobj.kws.get("itemid2name")
        prt.write("\n{HDR}\n".format(HDR=self.str_hdr()))
        for geneid in geneids:
            symbol = "" if itemid2name is None else itemid2name.get(geneid, "")
            prt.write(patgene.format(AART=self.gene2aart[geneid], ID=geneid, NAME=symbol))
304549
def kma(inputfile_1, out_path, databases, db_path_kma, min_cov=0.6,
            threshold=0.9, kma_path="cge/kma/kma", sample_name="",
            inputfile_2=None, kma_mrs=None, kma_gapopen=None,
            kma_gapextend=None, kma_penalty=None, kma_reward=None, kma_pm=None,
            kma_fpm=None, kma_memmode=False, kma_nanopore=False, debug=False,
            kma_add_args=None):
        
        threshold = threshold * 100
        min_cov = min_cov * 100

        kma_results = dict()
        kma_results["excluded"] = dict()

        if(sample_name):
           sample_name = "_" + sample_name

        # Initiate output dicts.
        gene_align_sbjct = {}
        gene_align_query = {}
        gene_align_homo = {}

        for db in databases:
            kma_db = db_path_kma + "/" + db
            kma_outfile = out_path + "/kma_" + db + sample_name
            kma_cmd = ("%s -t_db %s -o %s -e 1.0" % (kma_path,
                       kma_db, kma_outfile))
            if(inputfile_2 is not None):
                kma_cmd += " -ipe " + inputfile_1 + " " + inputfile_2
            else:
                kma_cmd += " -i " + inputfile_1
            if(kma_mrs is not None):
                kma_cmd += " -mrs " + str(kma_mrs)
            if(kma_gapopen is not None):
                kma_cmd += " -gapopen " + str(kma_gapopen)
            if(kma_gapextend is not None):
                kma_cmd += " -gapextend " + str(kma_gapextend)
            if(kma_penalty is not None):
                kma_cmd += " -penalty " + str(kma_penalty)
            if(kma_reward is not None):
                kma_cmd += " -reward " + str(kma_reward)
            if(kma_pm is not None):
                kma_cmd += " -pm " + kma_pm
            if(kma_fpm is not None):
                kma_cmd += " -fpm " + kma_fpm
            if (kma_memmode):
                kma_cmd += " -mem_mode "
            if (kma_nanopore):
                kma_cmd += " -bcNano "
                kma_cmd += " -mp 20 "
            if (kma_add_args is not None):
                kma_cmd += " " + kma_add_args + " "

            # kma output files
            align_filename = kma_outfile + ".aln"
            res_filename = kma_outfile + ".res"

            # If .res file exists then skip mapping
            if os.path.isfile(res_filename) and os.access(res_filename, os.R_OK):
                print("Found " + res_filename + " skipping DB.")
            else:
                # Call KMA
                if(debug):
                    print("KMA cmd: " + kma_cmd)

                process = subprocess.Popen(kma_cmd, shell=True,
                                           stdout=subprocess.PIPE,
                                           stderr=subprocess.PIPE)
                out, err = process.communicate()

            kma_results[db] = 'No hit found'

            # Open res file
            try:
                res_file = open(res_filename, "r")
                header = res_file.readline()
            except IOError as error:
                sys.exit("Error: KMA did not run as expected.\n" +
                         "KMA finished with the following response:" +
                         "\n{}\n{}".format(out.decode("utf-8"),
                                           err.decode("utf-8")))

            for line in res_file:

                if kma_results[db] == 'No hit found':
                    kma_results[db] = dict()
                    # kma_results[db]["excluded"] = dict()
                    # continue

                data = [data.strip() for data in line.split("\t")]
                gene = data[0]

                sbjct_len = int(data[3])
                sbjct_ident = float(data[4])
                coverage = float(data[5])
                depth = float(data[-3])
                q_value = float(data[-2])
                p_value = float(data[-1])


                if gene not in kma_results[db]:
                    hit = gene
                else:
                    hit = gene + "_" + str(len(kma_results[db][gene]) + 1)

                exclude_reasons = []

                if(coverage < min_cov or sbjct_ident < threshold):
                    exclude_reasons.append(coverage)
                    exclude_reasons.append(sbjct_ident)

                if(exclude_reasons):
                    # kma_results[db]["excluded"][hit] = exclude_reasons
                    kma_results["excluded"][hit] = exclude_reasons

                kma_results[db][hit] = dict()
                kma_results[db][hit]['sbjct_length'] = sbjct_len
                kma_results[db][hit]["perc_coverage"] = coverage
                kma_results[db][hit]["sbjct_string"] = []
                kma_results[db][hit]["query_string"] = []
                kma_results[db][hit]["homo_string"] = []
                kma_results[db][hit]["sbjct_header"] = gene
                kma_results[db][hit]["perc_ident"] = sbjct_ident
                kma_results[db][hit]["query_start"] = "NA"
                kma_results[db][hit]["query_end"] = "NA"
                kma_results[db][hit]["contig_name"] = "NA"
                kma_results[db][hit]["HSP_length"] = ""
                kma_results[db][hit]["cal_score"] = q_value
                kma_results[db][hit]["depth"] = depth
                kma_results[db][hit]["p_value"] = p_value
            res_file.close()


            if kma_results[db] == 'No hit found':
                continue

            # Open align file
            with open(align_filename, "r") as align_file:
                hit_no = dict()
                gene = ""
                # Parse through alignments
                for line in align_file:
                    # Skip empty lines
                    if(not line.strip()):
                        continue

                    # Check when a new gene alignment start
                    if line.startswith("#"):
                        gene = line[1:].strip()

                        if gene not in hit_no:
                            hit_no[gene] = str(1)
                        else:
                            hit_no[gene] += str(int(hit_no[gene]) + 1)

                    else:
                        # Check if gene one of the user specified genes
                        if hit_no[gene] == '1':
                            hit = gene
                        else:
                            hit = gene + "_" + hit_no[gene]

                        if hit in kma_results[db]:
                            line_data = line.split("\t")[-1].strip()
                            if line.startswith("template"):
                                kma_results[db][hit]["sbjct_string"] += (
                                    [line_data])
                            elif line.startswith("query"):
                                kma_results[db][hit]["query_string"] += (
                                    [line_data])
                            else:
                                kma_results[db][hit]["homo_string"] += (
                                    [line_data])
                        else:
                            print(hit + " not in results: ", kma_results)

            # concatinate all sequences lists and find subject start
            # and subject end

            gene_align_sbjct[db] = {}
            gene_align_query[db] = {}
            gene_align_homo[db] = {}

            for hit in kma_results[db]:
                # if(hit == "excluded"):
                # continue
                align_sbjct = "".join(kma_results[db][hit]['sbjct_string'])
                align_query = "".join(kma_results[db][hit]['query_string'])
                align_homo = "".join(kma_results[db][hit]['homo_string'])

                # Extract only aligned sequences
                start = re.search("^-*(\w+)", align_query).start(1)
                end = re.search("\w+(-*)$", align_query).start(1)

                kma_results[db][hit]['sbjct_string'] = align_sbjct[start:end]
                kma_results[db][hit]['query_string'] = align_query[start:end]
                kma_results[db][hit]['homo_string'] = align_homo[start:end]

                # Save align start and stop positions relative to
                # subject sequence
                kma_results[db][hit]['sbjct_start'] = start + 1
                kma_results[db][hit]["sbjct_end"] = end + 1
                kma_results[db][hit]["HSP_length"] = end - start

                # Count gaps in the alignment
                kma_results[db][hit]["gaps"] = (
                    kma_results[db][hit]['sbjct_string'].count("-") +
                    kma_results[db][hit]['query_string'].count("-"))

                # Save sequences covering the entire subject sequence
                # in seperate variables
                gene_align_sbjct[db][hit] = align_sbjct
                gene_align_query[db][hit] = align_query
                gene_align_homo[db][hit] = align_homo

        return FinderResult(kma_results, gene_align_sbjct, gene_align_query,
                            gene_align_homo)
683655
def dnld_annofile(fin_anno, anno_type):
    """"""
    if os.path.exists(fin_anno):
        return
    anno_type = get_anno_desc(fin_anno, anno_type)
    if anno_type == 'gene2go':
        dnld_ncbi_gene_file(fin_anno)
    if anno_type in {'gaf', 'gpad'}:
        dnld_annotation(fin_anno)