Create a report of the ePCR-calculated toxin profiles
168361
def epcr_report(self):
        
        logging.info('Creating {at} report'.format(at=self.analysistype))
        with open(os.path.join(self.reportpath, '{at}.csv'.format(at=self.analysistype)), 'w') as report:
            data = 'Strain,ToxinProfile\n'
            for sample in self.metadata:
                data += '{sn},{tp}\n'.format(sn=sample.name,
                                             tp=sample[self.analysistype].toxinprofile)
            # Write the data to the report
            report.write(data)
754114
def methodreporter(self):
        
        # Ensure that the analyses are set to complete
        self.analysescomplete = True
        # Reset the report path to original value
        self.reportpath = os.path.join(self.path, 'reports')
        # Clear the runmetadata - it will be populated with all the metadata from completemetadata
        self.runmetadata = MetadataObject()
        self.runmetadata.samples = list()
        # As the samples were entered into self.completemetadata depending on when they passed the quality threshold,
        # this list is not ordered numerically/alphabetically like the original runmetadata. Reset the order.
        for strain in self.samples:
            for sample in self.completemetadata:
                if sample.name == strain:
                    # Append the sample to the ordered list of objects
                    self.runmetadata.samples.append(sample)
        # Create the reports
        self.reporter()
        self.genusspecific()
        self.sixteensreporter()
        self.gdcsreporter()
        self.confindr_reporter()
515596
def report(args):
    
    logger.info("reading sequeces")
    data = load_data(args.json)

    logger.info("create profile")
    data = make_profile(data, os.path.join(args.out, "profiles"), args)
    logger.info("create database")
    make_database(data, "seqcluster.db", args.out)

    logger.info("Done. Download https://github.com/lpantano/seqclusterViz/archive/master.zip to browse the output.")
183115
def create_html(self):
        """"""

        roc_dicts = []
        with open(self._roc_fn, "r") as roc:
            for line in roc:
                line = line.strip()
                if line != "" and line[0] != "#":
                    (q, M, w, m, P, U, u, T, t, x, a) = map(int, line.split("\t"))
                    roc_dict = {
                        "q": q,
                        "M": M,
                        "w": w,
                        "m": m,
                        "P": P,
                        "U": U,
                        "u": u,
                        "T": T,
                        "t": t,
                        "x": x,
                        "a": a,
                        "mapped": M + w + m + P,
                        "unmapped": U + u + T + t + x,
                    }
                    roc_dicts.append(roc_dict)
        tbody = os.linesep.join(
            [
                """
				<tr>
					<td>        {quality}                   </td>
					<td>        {mapped}                    </td>
					<td><small> {mapped_proc:.2f}   </small></td>
					<td>        {M}                         </td>
					<td><small> {M_proc:.2f}        </small></td>
					<td>        {w}                         </td>
					<td><small> {w_proc:.2f}        </small></td>
					<td>        {m}                         </td>
					<td><small> {m_proc:.2f}        </small></td>
					<td>        {P}                         </td>
					<td><small> {P_proc:.2f}        </small></td>
					<td>        {unmapped}                  </td>
					<td><small> {unmapped_proc:.2f} </small></td>
					<td>        {U}                         </td>
					<td><small> {U_proc:.2f}        </small></td>
					<td>        {u}                         </td>
					<td><small> {u_proc:.2f}        </small></td>
					<td>        {T}                         </td>
					<td><small> {T_proc:.2f}        </small></td>
					<td>        {t}                         </td>
					<td><small> {t_proc:.2f}        </small></td>
					<td>        {x}                         </td>
					<td><small> {x_proc:.2f}        </small></td>
					<td>        {sum}                       </td>
					<td>        {prec_proc:.3f}              </td>
				</tr>
			""".format(
                    quality=roc_dict["q"],
                    mapped=roc_dict["mapped"],
                    mapped_proc=100.0 * (roc_dict["mapped"]) / roc_dict["a"],
                    M=roc_dict["M"],
                    M_proc=100.0 * (roc_dict["M"]) / (roc_dict["mapped"]) if (roc_dict["mapped"]) != 0 else 0,
                    w=roc_dict["w"],
                    w_proc=100.0 * (roc_dict["w"]) / (roc_dict["mapped"]) if (roc_dict["mapped"]) != 0 else 0,
                    m=roc_dict["m"],
                    m_proc=100.0 * (roc_dict["m"]) / (roc_dict["mapped"]) if (roc_dict["mapped"]) != 0 else 0,
                    P=roc_dict["P"],
                    P_proc=100.0 * (roc_dict["P"]) / (roc_dict["mapped"]) if (roc_dict["mapped"]) != 0 else 0,
                    unmapped=roc_dict["unmapped"],
                    unmapped_proc=100.0 * (roc_dict["unmapped"]) / roc_dict["a"],
                    U=roc_dict["U"],
                    U_proc=100.0 * (roc_dict["U"]) / (roc_dict["unmapped"]) if (roc_dict["unmapped"]) != 0 else 0,
                    u=roc_dict["u"],
                    u_proc=100.0 * (roc_dict["u"]) / (roc_dict["unmapped"]) if (roc_dict["unmapped"]) != 0 else 0,
                    T=roc_dict["T"],
                    T_proc=100.0 * (roc_dict["T"]) / (roc_dict["unmapped"]) if (roc_dict["unmapped"]) != 0 else 0,
                    t=roc_dict["t"],
                    t_proc=100.0 * (roc_dict["t"]) / (roc_dict["unmapped"]) if (roc_dict["unmapped"]) != 0 else 0,
                    x=roc_dict["x"],
                    x_proc=100.0 * (roc_dict["x"]) / (roc_dict["unmapped"]) if (roc_dict["unmapped"]) != 0 else 0,
                    sum=roc_dict["a"],
                    prec_proc=100.0 * (roc_dict["M"]) / (roc_dict["mapped"]) if (roc_dict["mapped"]) != 0 else 0,
                ) for roc_dict in roc_dicts
            ]
        )

        with open(self._html_fn, "w+") as html:
            program_info = ["No information available (PG header is essing)."]
            for x in rnftools.utils.shell(
                '"{samtools}" view -H "{bam}"'.format(
                    samtools="samtools",
                    bam=self._bam_fn,
                ), iterable=True
            ):
                x = x.strip()

                if x[:3] == "@PG":
                    pg_header = x[4:].strip()
                    parts = pg_header.split("\t")

                    program_info = ["<table>"]
                    for part in parts:
                        (lvalue, _, rvalue) = part.partition(":")
                        program_info.append(
                            '<tr><td style="font-weight:bold">{}:</td><td style="text-align:left">{}</td></tr>'.format(
                                lvalue, rvalue
                            )
                        )

                    program_info.append("</table>")

            eval_info = ["<table>"]
            for (lvalue, rvalue) in [
                ("BAM file", self._bam_fn),
                ("Allowed delta", self.report.allowed_delta),
            ]:
                eval_info.append(
                    '<tr><td style="font-weight:bold">{}:</td><td style="text-align:left">{}</td></tr>'.format(
                        lvalue, rvalue
                    )
                )

            eval_info.append("</table>")

            css_src = textwrap.dedent(
                """
					table            {border-collapse:collapse;}
					td               {border: solid #aaaaff 1px;padding:4px;text-align:right;}
					colgroup, thead  {border: solid black 2px;padding 2px;}
					.link_to_top     {font-size:10pt;}
					.desc            {color:#aaa;}
					.formats         {text-align:left;margin-bottom:20px;}
				"""
            )

            html_src = """<!DOCTYPE html>
				<html>
				<head>
					<meta charset="UTF-8" />
					<title>{name}</title>
					<style type="text/css">
					{css}
					</style>
				</head>
				<body>
					<h1 id="top">
						{name}
						<span class="link_to_top">
							[<a href="{homepage}">Back to main report</a>]
						</span>
					</h1>
	
					<p>
						<a href="#info_eval">Information about evaluation</a> -
						<a href="#info_prog">Information about program</a> -
						<a href="#roctable">ROC table</a> -
						<a href="#graphs">Graphs</a>
					</p>
	
					<h2 id="info_eval">
						Information about evaluation
						{headline_links}
					</h2>
	
					{eval_info}
	
					<h2 id="info_prog">
						Information about program
						{headline_links}
					</h2>
	
	
					{program_info}
	
					<h2 id="roctable">
						ROC table
						{headline_links}
					</h2>
	
					<p style="font-size:80%">
						<strong>M</strong>:
							mapped correctly
							<span class="desc">
								(all segments of tuple are mapped once and correctly),
							</span>
						<strong>w</strong>:
							mapped to wrong position
							<span class="desc">
								(at least one segment was mapped to a wrong position),
							</span>
						<strong>m</strong>:
							mapped but should be unmapped
							<span class="desc">
								(at least one segment was mapped but the read should not be mapped),
							</span>
						<strong>P</strong>:
							multimapped
							<span class="desc">
								(read should be mapped but it at least one segment was mapped several
								times and all segments were mapped correctly at least once),
							</span>
						<strong>U</strong>:
							unmapped correctly
							<span class="desc">
								(all segments of the read were correctly marked as unmapped),
							</span>
						<strong>u</strong>:
							unmapped but should be mapped
							<span class="desc">
								(at least one segment was mapped but entire read should be unmapped),
							</span>
						<strong>T</strong>:
							thresholded correctly
							<span class="desc">
								(read shoud not be mapped),
							</span>
						<strong>t</strong>:
							thresholded incorrectly
							<span class="desc">
								(read should be mapped),
							</span>
						<strong>x</strong>:
							unknown
							<span class="desc">
								(read is probably not reported by mapper)
							</span>
					</p>
					<table>
						<colgroup span="1" style="">
						<colgroup span="2" style="background-color:#ddd">
						<colgroup span="8" style="">
						<colgroup span="2" style="background-color:#ddd">
						<colgroup span="10" style="">
						<colgroup span="1" style="background-color:#ddd">
						<colgroup span="1" style="">
						<thead style="font-weight:bold;background-color:#ddddff">
							<tr style="font-weight:bold;background-color:#ddddff">
								<td>q</td>
								<td>mapped</td>
								<td>%</td>
								<td>M</td>
								<td>%</td>
								<td>w</td>
								<td>%</td>
								<td>m</td>
								<td>%</td>
								<td>P</td>
								<td>%</td>
								<td>unmapped</td>
								<td>%</td>
								<td>U</td>
								<td>%</td>
								<td>u</td>
								<td>%</td>
								<td>T</td>
								<td>%</td>
								<td>t</td>
								<td>%</td>
								<td>x</td>
								<td>%</td>
								<td>sum</td>
								<td>prec. (%)</td>
							</tr>
						</thead>
						<tbody>
							{tbody}
						</tbody>
					</table>
	
					<h2 id="graphs">
						Graphs
						{headline_links}
					</h2>
	
					<div class="formats">
						<img src="{svg}" />
						<br />
						<a href="{svg}">SVG version</a>
						|
						<a href="{roc}" type="text/csv">ROC file</a>
						|
						<a href="{gp}" type="text/plain">GP file</a>
					</div>
	
	
				</body>
				</html>
				""".format(
                name=self.name, tbody=tbody, css=css_src,
                svg=os.path.relpath(self._svg_fn, os.path.dirname(self._html_fn)),
                roc=os.path.relpath(self._roc_fn, os.path.dirname(self._html_fn)),
                gp=os.path.relpath(self._gp_fn, os.path.dirname(self._html_fn)), eval_info=os.linesep.join(eval_info),
                program_info=os.linesep.join(program_info),
                homepage=os.path.relpath(self.report.html_fn(), os.path.dirname(self._html_fn)), headline_links='''
							<span class="link_to_top">
								[<a href="#top">Top of this page</a>]
								[<a href="{homepage}">Back to main report</a>]
							</span>
						'''.format(homepage=os.path.relpath(self.report.html_fn(), os.path.dirname(self._html_fn)), )
            )

            tidy_html_src = bs4.BeautifulSoup(html_src).prettify()
            html.write(tidy_html_src)
492294
def process_pdb(pdbfile, outpath, as_string=False, outputprefix='report'):
    """"""
    if not as_string:
        startmessage = '\nStarting analysis of %s\n' % pdbfile.split('/')[-1]
    else:
        startmessage = "Starting analysis from stdin.\n"
    write_message(startmessage)
    write_message('='*len(startmessage)+'\n')
    mol = PDBComplex()
    mol.output_path = outpath
    mol.load_pdb(pdbfile, as_string=as_string)
    # #@todo Offers possibility for filter function from command line (by ligand chain, position, hetid)
    for ligand in mol.ligands:
        mol.characterize_complex(ligand)

    create_folder_if_not_exists(outpath)

    # Generate the report files
    streport = StructureReport(mol, outputprefix=outputprefix)

    config.MAXTHREADS = min(config.MAXTHREADS, len(mol.interaction_sets))

    ######################################
    # PyMOL Visualization (parallelized) #
    ######################################

    if config.PYMOL or config.PICS:
        try:
            from plip.modules.visualize import visualize_in_pymol
        except ImportError:
            from modules.visualize import visualize_in_pymol
        complexes = [VisualizerData(mol, site) for site in sorted(mol.interaction_sets)
                     if not len(mol.interaction_sets[site].interacting_res) == 0]
        if config.MAXTHREADS > 1:
            write_message('\nGenerating visualizations in parallel on %i cores ...' % config.MAXTHREADS)
            parfn = parallel_fn(visualize_in_pymol)
            parfn(complexes, processes=config.MAXTHREADS)
        else:
            [visualize_in_pymol(plcomplex) for plcomplex in complexes]

    if config.XML:  # Generate report in xml format
        streport.write_xml(as_string=config.STDOUT)

    if config.TXT:  # Generate report in txt (rst) format
        streport.write_txt(as_string=config.STDOUT)
50296
def reporter(self):
        
        # Initialise variables
        combinedrow = ''
        reportdirset = set()
        # Populate a set of all the report directories to use. A standard analysis will only have a single report
        # directory, while pipeline analyses will have as many report directories as there are assembled samples
        for sample in self.metadata:
            # Ignore samples that lack a populated reportdir attribute
            if sample[self.analysistype].reportdir != 'NA':
                make_path(sample[self.analysistype].reportdir)
                # Add to the set - I probably could have used a counter here, but I decided against it
                reportdirset.add(sample[self.analysistype].reportdir)
        # Create a report for each sample from :self.resultprofile
        for sample in self.metadata:
            if sample[self.analysistype].reportdir != 'NA':
                if type(sample[self.analysistype].allelenames) == list:
                    # Populate the header with the appropriate data, including all the genes in the list of targets
                    row = 'Strain,Genus,SequenceType,Matches,{},\n' \
                        .format(','.join(sorted(sample[self.analysistype].allelenames)))
                    # Set the sequence counter to 0. This will be used when a sample has multiple best sequence types.
                    # The name of the sample will not be written on subsequent rows in order to make the report clearer
                    seqcount = 0
                    # Iterate through the best sequence types for the sample (only occurs if update profile is disabled)
                    for seqtype in self.resultprofile[sample.name]:
                        """
                        {
                            "OLF15230-1_2015-SEQ-0783": {
                                "1000004_CFIA": {
                                    "7": {
                                        "dnaE": {
                                            "47": "100.00"
                                        },
                                        "dtdS": {
                                            "19": "100.00"
                                        },
                                        "gyrB": {
                                            "359": "100.00"
                                        },
                                        "pntA": {
                                            "50": "100.00"
                                        },
                                        "pyrC": {
                                            "143": "100.00"
                                        },
                                        "recA": {
                                            "31": "100.00"
                                        },
                                        "tnaA": {
                                            "26": "100.00"
                                        }
                                    }
                                }
                            }
                        }
                        """
                        # Becomes
                        """
                        Strain,SequenceType,Matches,dnaE,gyrB,recA,dtdS,pntA,pyrC,tnaA
                        OLF15230-1_2015-SEQ-0783,1000004_CFIA,7,26 (100.00%),359 (100.00%),31 (100.00%),50 (100.00%),
                            19 (100.00%),47 (100.00%),143 (100.00%)
                        """
                        sample[self.analysistype].sequencetype = seqtype
                        # The number of matches to the profile
                        matches = list(self.resultprofile[sample.name][seqtype].keys())[0]
                        # If this is the first of one or more sequence types, include the sample name
                        if seqcount == 0:
                            row += '{},{},{},{},'.format(sample.name, sample.general.referencegenus, seqtype, matches)
                        # Otherwise, skip the sample name
                        else:
                            row += ',,{},{},'.format(seqtype, matches)
                        # Iterate through all the genes present in the analyses for the sample
                        for gene in sorted(sample[self.analysistype].allelenames):
                            # refallele = self.profiledata[self.analysistype][seqtype][gene]
                            refallele = sample[self.analysistype].profiledata[seqtype][gene]
                            # Set the allele and percent id from the dictionary's keys and values, respectively
                            allele = list(self.resultprofile[sample.name][seqtype][matches][gene].keys())[0]
                            percentid = list(self.resultprofile[sample.name][seqtype][matches][gene].values())[0]
                            if refallele and refallele != allele:
                                if 0 < float(percentid) < 100:
                                    row += '{} ({:.2f}%),'.format(allele, float(percentid))
                                else:
                                    row += '{} ({}),'.format(allele, refallele)
                            else:
                                # Add the allele and % id to the row (only add the percent identity if it is not 100%)
                                if 0 < float(percentid) < 100:
                                    row += '{} ({:.2f}%),'.format(allele, float(percentid))
                                else:
                                    row += '{},'.format(allele)
                            self.referenceprofile[sample.name][gene] = allele
                        # Add a newline
                        row += '\n'
                        # Increment the number of sequence types observed for the sample
                        seqcount += 1
                    combinedrow += row
                    # If the length of the # of report directories is greater than 1 (script is being run as part of
                    # the assembly pipeline) make a report for each sample
                    if self.pipeline:
                        # Open the report
                        with open('{}{}_{}.csv'.format(sample[self.analysistype].reportdir, sample.name,
                                                       self.analysistype), 'w') as report:
                            # Write the row to the report
                            report.write(row)
                dotter()
            # Create the report folder
            make_path(self.reportpath)
            # Create the report containing all the data from all samples
            if self.pipeline:
                with open('{}{}.csv'.format(self.reportpath, self.analysistype), 'w') \
                        as combinedreport:
                    # Write the results to this report
                    combinedreport.write(combinedrow)
            else:
                with open('{}{}_{:}.csv'.format(self.reportpath, self.analysistype, time.strftime("%Y.%m.%d.%H.%M.%S")),
                          'w') as combinedreport:
                    # Write the results to this report
                    combinedreport.write(combinedrow)
            # Remove the raw results csv
            [os.remove(rawresults) for rawresults in glob('{}*rawresults*'.format(self.reportpath))]
749122
def filter_for_ligands(self):
        """"""

        candidates1 = [o for o in pybel.ob.OBResidueIter(
            self.proteincomplex.OBMol) if not o.GetResidueProperty(9) and self.is_het_residue(o)]

        if config.DNARECEPTOR:  # If DNA is the receptor, don't consider DNA as a ligand
            candidates1 = [res for res in candidates1 if res.GetName() not in config.DNA+config.RNA]
        all_lignames = set([a.GetName() for a in candidates1])

        water = [o for o in pybel.ob.OBResidueIter(self.proteincomplex.OBMol) if o.GetResidueProperty(9)]
        # Filter out non-ligands
        if not config.KEEPMOD:  # Keep modified residues as ligands
            candidates2 = [a for a in candidates1 if is_lig(a.GetName()) and a.GetName() not in self.modresidues]
        else:
            candidates2 = [a for a in candidates1 if is_lig(a.GetName())]
        write_message("%i ligand(s) after first filtering step.\n" % len(candidates2), mtype='debug')

        ############################################
        # Filtering by counting and artifacts list #
        ############################################
        artifacts = []
        unique_ligs = set(a.GetName() for a in candidates2)
        for ulig in unique_ligs:
            # Discard if appearing 15 times or more and is possible artifact
            if ulig in config.biolip_list and [a.GetName() for a in candidates2].count(ulig) >= 15:
                artifacts.append(ulig)

        selected_ligands = [a for a in candidates2 if a.GetName() not in artifacts]

        return selected_ligands, all_lignames, water
1370284
def report(self, align_bam, ref_file, is_paired, bait_file, target_file,
               variant_region_file, config):
        
        dup_metrics = self._get_current_dup_metrics(align_bam)
        align_metrics = self._collect_align_metrics(align_bam, ref_file)
        # Prefer the GC metrics in FastQC instead of Picard
        # gc_graph, gc_metrics = self._gc_bias(align_bam, ref_file)
        gc_graph = None
        insert_graph, insert_metrics, hybrid_metrics = (None, None, None)
        if is_paired:
            insert_graph, insert_metrics = self._insert_sizes(align_bam)
        if bait_file and target_file:
            assert os.path.exists(bait_file), (bait_file, "does not exist!")
            assert os.path.exists(target_file), (target_file, "does not exist!")
            hybrid_metrics = self._hybrid_select_metrics(align_bam,
                                                         bait_file, target_file)
        elif (variant_region_file and
              config["algorithm"].get("coverage_interval", "").lower() in ["exome"]):
            assert os.path.exists(variant_region_file), (variant_region_file, "does not exist")
            hybrid_metrics = self._hybrid_select_metrics(
                align_bam, variant_region_file, variant_region_file)

        vrn_vals = self._variant_eval_metrics(align_bam)
        summary_info = self._parser.get_summary_metrics(align_metrics,
                dup_metrics, insert_metrics, hybrid_metrics,
                vrn_vals)
        graphs = []
        if gc_graph and os.path.exists(gc_graph):
            graphs.append((gc_graph, "Distribution of GC content across reads"))
        if insert_graph and os.path.exists(insert_graph):
            graphs.append((insert_graph, "Distribution of paired end insert sizes"))
        return summary_info, graphs
142000
def profiler(self):
        """"""
        # Initialise variables
        profiledata = defaultdict(make_dict)
        profileset = set()
        # supplementalset = ''
        genedict = {}
        # Find all the unique profiles to use with a set
        for sample in self.metadata:
            if sample[self.analysistype].profile != 'NA':
                profileset.add(sample[self.analysistype].profile[0])
                # if self.analysistype == 'rmlst':
                #     supplementalset = sample[self.analysistype].supplementalprofile
        # Extract the profiles for each set
        for sequenceprofile in profileset:
            # Clear the list of genes
            genelist = []
            for sample in self.metadata:
                if sequenceprofile == sample[self.analysistype].profile[0]:
                    # genelist = [os.path.split(x)[1].split('.')[0] for x in sample[self.analysistype].alleles]
                    genelist = sample[self.analysistype].allelenames
            try:
                # Open the sequence profile file as a dictionary
                profile = DictReader(open(sequenceprofile), dialect='excel-tab')
            # 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:
                        try:
                            profiledata[sequenceprofile][row['rST']][gene] = row[gene]
                        except KeyError:
                            raise

            # # Load the supplemental profile definitions
            # if self.analysistype == 'rmlst':
            #     supplementalprofile = DictReader(open(supplementalset), dialect='excel-tab')
            #     # Do the same with the supplemental profile
            #     for row in supplementalprofile:
            #         # Iterate through the genes
            #         for gene in genelist:
            #             # Add the sequence profile, and type, the gene name and the allele number to the dictionary
            #             profiledata[sequenceprofile][row['rST']][gene] = row[gene]
            # 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:
                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()
877526
def report(self):
        
        # Add all the genes to the header
        header = 'Sample,'
        data = str()
        with open(os.path.join(self.reportpath, '{at}.csv'.format(at=self.analysistype)), 'w') as report:
            write_header = True
            for sample in self.runmetadata:
                data += sample.name + ','
                # Iterate through all the user-supplied target names
                for target in sorted(self.genes):
                    write_results = False
                    # There was an issue with 'target' not matching 'name' due to a dash being replaced by an underscore
                    # only in 'name'. This will hopefully address this issue
                    target = target.replace('-', '_')
                    if write_header:
                        header += '{target}_match_details,{target},'.format(target=target)
                    for name, identity in sample[self.analysistype].results.items():
                        # Ensure that all dashes are replaced with underscores
                        name = name.replace('-', '_')
                        # If the current target matches the target in the header, add the data to the string
                        if name == target:
                            write_results = True
                            gene_results = '{percent_id}% ({avgdepth} +/- {stddev}),{record},'\
                                .format(percent_id=identity,
                                        avgdepth=sample[self.analysistype].avgdepth[name],
                                        stddev=sample[self.analysistype].standarddev[name],
                                        record=sample[self.analysistype].sequences[target])
                            # Populate the data string appropriately
                            data += gene_results
                    # If the target is not present, write dashes to represent the results and sequence
                    if not write_results:
                        data += '-,-,'
                data += ' \n'
                write_header = False
            header += '\n'
            # Write the strings to the report
            report.write(header)
            report.write(data)
747047
def reporter(self):
        
        logging.info('Creating {at} report'.format(at=self.analysistype))
        resistance_classes = ResistanceNotes.classes(self.targetpath)
        # Find unique gene names with the highest percent identity
        for sample in self.runmetadata.samples:
            try:
                if sample[self.analysistype].results:
                    # Initialise a dictionary to store the unique genes, and their percent identities
                    sample[self.analysistype].uniquegenes = dict()
                    for name, identity in sample[self.analysistype].results.items():
                        # Split the name of the gene from the string e.g. ARR-2_1_HQ141279 yields ARR-2
                        genename = name.split('_')[0]
                        # Set the best observed percent identity for each unique gene
                        try:
                            # Pull the previous best identity from the dictionary
                            bestidentity = sample[self.analysistype].uniquegenes[genename]
                            # If the current identity is better than the old identity, save it
                            if float(identity) > float(bestidentity):
                                sample[self.analysistype].uniquegenes[genename] = float(identity)
                        # Initialise the dictionary if necessary
                        except KeyError:
                            sample[self.analysistype].uniquegenes[genename] = float(identity)
            except AttributeError:
                pass
        # Create the path in which the reports are stored
        make_path(self.reportpath)
        # Initialise strings to store the results
        data = 'Strain,Resistance,Gene,Allele,Accession,PercentIdentity,Length,FoldCoverage\n'
        with open(os.path.join(self.reportpath, self.analysistype + '.csv'), 'w') as report:
            for sample in self.runmetadata.samples:
                # Create an attribute to store the string for the eventual pipeline report
                sample[self.analysistype].pipelineresults = list()
                data += sample.name + ','
                if sample[self.analysistype].results:
                    # If there are multiple results for a sample, don't write the name in each line of the report
                    multiple = False
                    for name, identity in sorted(sample[self.analysistype].results.items()):
                        # Extract the necessary variables from the gene name string
                        gname, genename, accession, allele = ResistanceNotes.gene_name(name)
                        # Retrieve the best identity for each gene
                        try:
                            percentid = sample[self.analysistype].uniquegenes[gname]
                        # Beta-lactamases will not have the allele and version from the gene name defined above
                        except KeyError:
                            percentid = sample[self.analysistype].uniquegenes[gname.split('-')[0]]
                        # If the percent identity of the current gene matches the best percent identity, add it to
                        # the report - there can be multiple occurrences of genes e.g.
                        # sul1,1,AY224185,100.00,840 and sul1,2,CP002151,100.00,927 are both included because they
                        # have the same 100% percent identity
                        if float(identity) == percentid:
                            try:
                                # Determine resistance phenotype of the gene
                                res = ResistanceNotes.resistance(name, resistance_classes)
                                # Treat the initial vs subsequent results for each sample slightly differently - instead
                                # of including the sample name, use an empty cell instead
                                if multiple:
                                    data += ','
                                # Populate the results
                                data += '{},{},{},{},{},{},{}\n'.format(
                                    res,
                                    genename,
                                    allele,
                                    accession,
                                    identity,
                                    len(sample[self.analysistype].sequences[name]),
                                    sample[self.analysistype].avgdepth[name])
                                sample[self.analysistype].pipelineresults.append(
                                    '{rgene} ({pid}%) {rclass}'.format(rgene=genename,
                                                                       pid=identity,
                                                                       rclass=res)
                                )
                                multiple = True
                            except KeyError:
                                pass
                else:
                    data += '\n'
            # Write the strings to the file
            report.write(data)
267254
def plot_profile(ribo_counts, transcript_name, transcript_length,
                 start_stops, read_lengths=None, read_offsets=None, rna_counts=None,
                 color_scheme='default', html_file='index.html', output_path='output'):
    
    colors = get_color_palette(scheme=color_scheme)
    gs = gridspec.GridSpec(3, 1, height_ratios=[6, 1.3, 0.5], hspace=0.35)
    font_axis = {'family': 'sans-serif', 'color': colors['color'], 'weight': 'bold', 'size': 7}

    # riboseq bar plots
    gs2 = gridspec.GridSpecFromSubplotSpec(1, 1, subplot_spec=gs[0])
    ax2 = plt.subplot(gs2[0])
    label = 'Ribo-Seq count'

    if read_lengths:
        if len(read_lengths) > 1:
            label = 'Ribo-Seq count ({}-mers)'.format(', '.join('{}'.format(item) for item in read_lengths))
        else:
            label = 'Ribo-Seq count ({}-mer)'.format('{}'.format(read_lengths[0]))

    ax2.set_ylabel(label, fontdict=font_axis, labelpad=10)

    # rna coverage if available
    ax_rna = None
    if rna_counts:
        ax_rna = ax2.twinx()
        ax_rna.set_ylabel('RNA-Seq count', fontdict=font_axis, labelpad=10)
        ax_rna.bar(rna_counts.keys(), rna_counts.values(), facecolor=colors['rna'],
                   edgecolor=colors['rna'], label='RNA')
        ax_rna.set_zorder(1)

    frame_counts = {1: {}, 2: {}, 3: {}}
    for k, v in ribo_counts.iteritems():
        for fr in (1, 2, 3):
            if v[fr] > 0:
                frame_counts[fr][k] = v[fr]
                break

    cnts = []
    [cnts.extend(item.values()) for item in frame_counts.values()]
    y_max = float(max(cnts) * 1.25)
    ax2.set_ylim(0.0, y_max)
    ax2.set_zorder(2)
    ax2.patch.set_facecolor('none')

    for frame in (1, 2, 3):
        color = colors['frames'][frame - 1]
        x_vals = frame_counts[frame].keys()
        ax2.bar(x_vals, frame_counts[frame].values(), color=color, facecolor=color, edgecolor=color)

    # ORF architecture
    gs3 = gridspec.GridSpecFromSubplotSpec(3, 1, subplot_spec=gs[1], hspace=0.1)
    if color_scheme == 'greyorfs':
        axisbg = [colors['grey'] for i in range(3)]
    else:
        axisbg = colors['frames']

    ax4 = plt.subplot(gs3[0], sharex=ax2, axisbg=axisbg[0])
    ax5 = plt.subplot(gs3[1], sharex=ax2, axisbg=axisbg[1])
    ax6 = plt.subplot(gs3[2], sharex=ax2, axisbg=axisbg[2])
    ax6.set_xlabel('Transcript length ({} nt)'.format(transcript_length), fontdict=font_axis, labelpad=6)

    # Legend
    gs4 = gridspec.GridSpecFromSubplotSpec(1, 1, subplot_spec=gs[2], hspace=0.1)
    ax7 = plt.subplot(gs4[0], axisbg=colors['background'])
    set_axis_color(ax7, colors['background'])

    ax7.text(0.02, 0.1, "AUG", size=5, ha="center", va="center", color=colors['color'],
             bbox=dict(boxstyle="square", facecolor=colors['start'], edgecolor=colors['color'], linewidth=0.3))
    ax7.text(0.06, 0.1, "STOP", size=5, ha="center", va="center", color='white',
             bbox=dict(boxstyle="square", color=colors['stop']))
    ax7.text(0.13, 0.1, "Frames", size=5, ha='center', va='center', color=colors['color'],
             fontdict={'weight': 'bold'})
    ax7.text(0.17, 0.1, "1", size=5, ha="center", va="center", color='white',
             bbox=dict(boxstyle="square", color=colors['frames'][0]))
    ax7.text(0.19, 0.1, "2", size=5, ha="center", va="center", color='white',
             bbox=dict(boxstyle="square", color=colors['frames'][1]))
    ax7.text(0.21, 0.1, "3", size=5, ha="center", va="center", color='white',
             bbox=dict(boxstyle="square", color=colors['frames'][2]))

    # No ticks or labels for ORF 1, 2 and Legend
    for axis in (ax4, ax5, ax7):
        axis.tick_params(top=False, left=False, right=False, bottom=False, labeltop=False,
                         labelleft=False, labelright=False, labelbottom=False)

    axes = [ax2]
    if ax_rna:
        axes.append(ax_rna)

    fp = FontProperties(size='5')
    for axis in axes:
        set_axis_color(axis, colors['axis'])
        axis.tick_params(colors=colors['ticks'])
        for item in (axis.get_xticklabels() + axis.get_yticklabels()):
            item.set_fontproperties(fp)
            item.set_color(colors['color'])

    for axis, frame in ((ax4, 1), (ax5, 2), (ax6, 3)):
        if color_scheme == 'greyorfs':
            color = colors['grey']
        else:
            color = colors['frames'][frame - 1]
        set_axis_color(axis, color, alpha=0.05)
        axis.patch.set_alpha(0.3)  # opacity of ORF architecture
        for item in (axis.get_xticklabels()):
            item.set_fontproperties(fp)
            item.set_color(colors['color'])
        axis.set_ylim(0, 0.2)
        axis.set_xlim(0, transcript_length)
        starts = [(item, 1) for item in start_stops[frame]['starts']]
        stops = [(item, 1) for item in start_stops[frame]['stops']]
        start_colors = [colors['start'] for item in starts]
        axis.broken_barh(starts, (0.11, 0.2), facecolors=start_colors,
                         edgecolors=start_colors, label='start', zorder=5)
        stop_colors = [colors['stop'] for item in stops]
        axis.broken_barh(stops, (0, 0.2), facecolors=stop_colors,
                         edgecolors=stop_colors, label='stop', zorder=5)
        axis.set_ylabel('{}'.format(frame),
                        fontdict={'family': 'sans-serif', 'color': colors['color'],
                                  'weight': 'normal', 'size': '6'},
                        rotation='horizontal', labelpad=10, verticalalignment='center')
        axis.tick_params(top=False, left=False, right=False, labeltop=False,
                         labelleft=False, labelright=False, direction='out', colors=colors['ticks'])
    plt.title('{}'.format(transcript_name),
              fontdict={'family': 'sans-serif', 'color': colors['color'],
                        'weight': 'bold', 'size': 8, 'y': 20})
    if not os.path.exists(output_path):
        os.mkdir(output_path)
    plt.savefig(os.path.join(output_path, 'riboplot.svg'), facecolor=colors['background'])
    plt.savefig(os.path.join(output_path, 'riboplot.png'), dpi=600, facecolor=colors['background'])

    with open(os.path.join(CONFIG.PKG_DATA_DIR, 'riboplot.html')) as g, open(os.path.join(output_path, html_file), 'w') as h:
        h.write(g.read().format(transcript_name=transcript_name))

    css_dir = os.path.join(output_path, 'css')
    if not os.path.exists(css_dir):
        os.mkdir(css_dir)

    css_data_dir = os.path.join(CONFIG.PKG_DATA_DIR, 'css')
    for fname in os.listdir(css_data_dir):
        shutil.copy(os.path.join(css_data_dir, fname), os.path.join(output_path, 'css', fname))
167089
def results( self ):
        """"""
        return self.report(self.parameters(),
                           self.metadata(),
                           self.experimentalResults())
1299051
def reporter(self):
        
        # Create the path in which the reports are stored
        make_path(self.reportpath)
        header = 'Strain,Gene,PercentIdentity,Genus,FoldCoverage\n'
        data = ''
        with open(os.path.join(self.reportpath, self.analysistype + '.csv'), 'w') as report:
            for sample in self.runmetadata.samples:
                data += sample.name + ','
                if sample[self.analysistype].results:
                    if not sample[self.analysistype].multiple:
                        for name, identity in sample[self.analysistype].results.items():
                            if name == sample[self.analysistype].besthit[0]:
                                data += '{},{},{},{}\n'.format(name, identity, sample[self.analysistype].genera[name],
                                                               sample[self.analysistype].avgdepth[name])
                    else:
                        data += '{},{},{},{}\n'.format('multiple', 'NA', ';'.join(sample[self.analysistype]
                                                                                  .classification), 'NA')
                else:
                    data += '\n'
            report.write(header)
            report.write(data)
1036526
def reporter(self):
        
        logging.info('Creating {} report'.format(self.analysistype))
        make_path(self.reportpath)
        header = 'Strain,ReferenceGenus,ReferenceFile,ReferenceGenomeMashDistance,Pvalue,NumMatchingHashes\n'
        data = ''
        for sample in self.metadata:
            try:
                data += '{},{},{},{},{},{}\n'.format(sample.name,
                                                     sample[self.analysistype].closestrefseqgenus,
                                                     sample[self.analysistype].closestrefseq,
                                                     sample[self.analysistype].mashdistance,
                                                     sample[self.analysistype].pvalue,
                                                     sample[self.analysistype].nummatches)
            except AttributeError:
                data += '{}\n'.format(sample.name)
        # Create the report file
        reportfile = os.path.join(self.reportpath, 'mash.csv')
        with open(reportfile, 'w') as report:
            report.write(header)
            report.write(data)
178275
def report(self):
        """"""
        print('\n{}{}{}'.format(c.Style.BRIGHT, c.Fore.CYAN, 'Report:'))
        # report headers
        print('{!s:<85}{!s:<20}'.format('', 'Validations'))
        print(
            '{!s:<60}{!s:<25}{!s:<10}{!s:<10}'.format(
                'Profile:', 'Execution:', 'Passed:', 'Failed:'
            )
        )
        for r in self.reports:
            d = r.data
            if not d.get('selected'):
                continue

            # execution
            execution_color = c.Fore.RED
            execution_text = 'Failed'
            if d.get('execution_success'):
                execution_color = c.Fore.GREEN
                execution_text = 'Passed'

            # pass count
            pass_count_color = c.Fore.GREEN
            pass_count = d.get('validation_pass_count', 0)
            # fail count
            fail_count = d.get('validation_fail_count', 0)
            fail_count_color = c.Fore.GREEN
            if fail_count > 0:
                fail_count_color = c.Fore.RED

            # report row
            print(
                '{!s:<60}{}{!s:<25}{}{!s:<10}{}{!s:<10}'.format(
                    d.get('name'),
                    execution_color,
                    execution_text,
                    pass_count_color,
                    pass_count,
                    fail_count_color,
                    fail_count,
                )
            )

        # write report to disk
        if self.args.report:
            with open(self.args.report, 'w') as outfile:
                outfile.write(str(self.reports))
172530
def parse_reports(self):
    

    # Set up vars
    self.gene_body_cov_hist_counts = dict()
    self.gene_body_cov_hist_percent = dict()

    # TODO - Do separate parsing step to find skewness values
    # and add these to the general stats table?

    # Go through files and parse data
    for f in self.find_log_files('rseqc/gene_body_coverage'):

        # RSeQC >= v2.4
        if f['f'].startswith('Percentile'):
            keys = []
            nrows = 0
            for l in f['f'].splitlines():
                s = l.split()
                if len(keys) == 0:
                    keys = s[1:]
                else:
                    nrows += 1
                    s_name = self.clean_s_name(s[0], f['root'])
                    if s_name in self.gene_body_cov_hist_counts:
                        log.debug("Duplicate sample name found! Overwriting: {}".format(s_name))
                    self.add_data_source(f, s_name, section='gene_body_coverage')
                    self.gene_body_cov_hist_counts[s_name] = OrderedDict()
                    for k, var in enumerate(s[1:]):
                        self.gene_body_cov_hist_counts[s_name][int(keys[k])] = float(var)
            if nrows == 0:
                log.warning("Empty geneBodyCoverage file found: {}".format(f['fn']))


        # RSeQC < v2.4
        elif f['f'].startswith('Total reads'):
            if f['s_name'].endswith('.geneBodyCoverage'):
                f['s_name'] = f['s_name'][:-17]
            if f['s_name'] in self.gene_body_cov_hist_counts:
                log.debug("Duplicate sample name found! Overwriting: {}".format(f['s_name']))
            self.add_data_source(f, section='gene_body_coverage')
            self.gene_body_cov_hist_counts[f['s_name']] = OrderedDict()
            nrows = 0
            for l in f['f'].splitlines():
                s = l.split()
                try:
                    nrows += 1
                    self.gene_body_cov_hist_counts[f['s_name']][int(s[0])] = float(s[1])
                except:
                    pass
            if nrows == 0:
                del self.gene_body_cov_hist_counts[f['s_name']]
                log.warning("Empty geneBodyCoverage file found: {}".format(f['fn']))

    # Filter to strip out ignored sample names
    self.gene_body_cov_hist_counts = self.ignore_samples(self.gene_body_cov_hist_counts)

    if len(self.gene_body_cov_hist_counts) > 0:

        # Make a normalised percentage version of the data
        for s_name in self.gene_body_cov_hist_counts:
            self.gene_body_cov_hist_percent[s_name] = OrderedDict()
            total = sum( self.gene_body_cov_hist_counts[s_name].values() )
            for k, v in self.gene_body_cov_hist_counts[s_name].items():
                self.gene_body_cov_hist_percent[s_name][k] = (v/total)*100

        # Add line graph to section
        pconfig = {
            'id': 'rseqc_gene_body_coverage_plot',
            'title': 'RSeQC: Gene Body Coverage',
            'ylab': '% Coverage',
            'xlab': "Gene Body Percentile (5' -> 3')",
            'xmin': 0,
            'xmax': 100,
            'tt_label': "<strong>{point.x}% from 5'</strong>: {point.y:.2f}",
            'data_labels': [
                {'name': 'Percentages', 'ylab': 'Percentage Coverage'},
                {'name': 'Counts', 'ylab': 'Coverage'}
            ]
        }
        self.add_section (
            name = 'Gene Body Coverage',
            anchor = 'rseqc-gene_body_coverage',
            description = '<a href="http://rseqc.sourceforge.net/#genebody-coverage-py" target="_blank">Gene Body Coverage</a>' \
                " calculates read coverage over gene bodies." \
                " This is used to check if reads coverage is uniform and" \
                " if there is any 5' or 3' bias.",
            plot = linegraph.plot([self.gene_body_cov_hist_percent, self.gene_body_cov_hist_counts], pconfig)
        )

    # Return number of samples found
    return len(self.gene_body_cov_hist_counts)