ReportTools.py

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from matplotlib_venn import venn3
import pandas as pd
import seaborn as sns
from scipy.stats import ttest_ind, mannwhitneyu, pearsonr
import os.path
import math

from ParseResults import parse_regions, get_stats, parse_gene_file
from ensemblQueries import gene_name_to_ID, gene_ID_to_name

def make_report(gistic_results, rubic_gain_results, rubic_loss_results,
                  census_genes, known_genes, ref_genome, biomart_file,
                  file_tools, file_regions, file_venn, file_swarmplot,
                  file_genes_both, file_genes_GISTIC, file_genes_RUBIC,
                  bed_overlap):
    """Make a report, tables, figures and gene files from the analyses with GISTIC2 and RUBIC."""
    parsed, stats = [], []
    for tool in 'GISTIC', 'RUBIC', 'Overlap':
        if tool == 'GISTIC':
            results = gistic_results
        elif tool == 'RUBIC':
            results = [rubic_gain_results, rubic_loss_results]
        else:
            results = bed_overlap
        parsed_results = parse_regions(results, known_genes, census_genes, tool, ref_genome)
        parsed.append(parsed_results)
        stats_results = get_stats(parsed_results, tool)
        stats.append(stats_results)
    make_table_regions(parsed[0:2], file_regions, ref_genome)                                                                    #make table with all recurrent regions
    make_tool_report(file_tools, stats)                                                                         #make tool report
    extract_gene_lists(parsed, file_genes_GISTIC, file_genes_RUBIC, file_genes_both, known_genes, file_venn, ref_genome)    #make gene files and venn overlap plot
    plot_histogram_sizes(parsed, file_swarmplot)                                                                #plot histogram with the sizes of the regions

def make_tool_report(file_tools, stats_tools):
    """Make a report of the results from GISTIC and RUBIC."""
    row_names = ["Tool", "Type", "No. regions", "Avg. size (Kb)", "Total size (Mb)", "No. genes", "Average driver density", "No. focal regions", "No. known regions",
                 "No. census regions", "No. known focal regions"]
    with open(file_tools, 'w') as out:
        for i in range(len(row_names)):
            list_out = [row_names[i]]
            for j in range(len(stats_tools)):
                list_out += [stats_tools[j][0][i], stats_tools[j][1][i]]
            list_out = list(map(str, list_out))
            out.write("\t".join(list_out) + "\n")

def make_table_regions(parsed_tools, file_regions, ref_genome):
    """Make a table with all recurent regions detected by both tools."""
    with open(file_regions, 'w') as out:
        header = ["Method", "Chr", "Start", "End", "Type", "Negative log10 q-value", "Known genes", "Census genes", "Total number of genes"]
        out.write("\t".join(header) + "\n")
        for tool in parsed_tools:
            tool_name = "GISTIC" if tool == parsed_tools[0] else "RUBIC"
            for cnv in tool:
                out.write(tool_name + "\t")
                cnv_type = "Amplification" if cnv[3] == "amp" else "Deletion"
                q_vals = ", ".join(cnv[4]) if type(cnv[4]) == tuple else cnv[4]
                begin_row = cnv[0:3] + [cnv_type, str(q_vals)]
                out.write("\t".join(begin_row))
                for list_genes in range(5,7):
                    out_genes = gene_ID_to_name(cnv[list_genes], ref_genome)
                    out.write("\t" + ", ".join(out_genes))
                out.write("\t" + str(len(cnv[7])))
                #for list_genes in range(5,8):
                    #out.write("\t" + ", ".join(cnv[list_genes]))
                out.write("\n")

def extract_gene_lists(parsed_results, out_GISTIC, out_RUBIC, out_both, known_genes, out_venn, ref_genome):
    """Extract lists of genes from the parsed results and a list of known genes to make files with gene lists and a venn diagram on the overlap."""
    genes_all = []
    genes_venn = []
    for tool_results in parsed_results:     #tool
        genes_tool = []
        for cnv_type in tool_results:       #type CNV
            for cnv in tool_results:
                genes = cnv[7]
                genes_tool += genes
        unique_genes = list(set(genes_tool))
        genes_all.append(unique_genes)
        out_file = out_GISTIC if tool_results == parsed_results[0] else out_RUBIC
        make_gene_file_tool(genes_tool, out_file)
    make_gene_file_overlap(genes_all, out_both)
    venn3_overlap(genes_all, known_genes, out_venn, ref_genome)

def make_gene_file_tool(gene_names, out_file):
    """Make file with list of genes reported by tool (one gene per line)."""
    with open(out_file, 'w') as out:
        #gene_IDs = gene_name_to_ID(gene_names)       #Convert gene names to IDs in order to do GO enrichment analysis
        #out.write("\n".join(gene_IDs))
        out.write("\n".join(gene_names))

def make_gene_file_overlap(gene_lists, out_both):
    """Make file with list of all genes reported by both tools (one gene per line)."""
    overlapping_genes = set(gene_lists[0]).intersection(set(gene_lists[1]))
    #overlapping_IDs = gene_name_to_ID(list(overlapping_genes))
    with open(out_both, 'w') as overlap_file:
        #overlap_file.write("\n".join(overlapping_IDs))
        overlap_file.write("\n".join(overlapping_genes))

def venn3_overlap(gene_lists, known_genes_file, out_venn, ref_genome):
    """Plot the overlap between the gene lists per tool and the list of known genes."""
    known_genes = parse_gene_file(known_genes_file, np.inf, ref_genome)
    overlap_tools = set(gene_lists[0]).intersection(set(gene_lists[1]))
    overlap_all = list(overlap_tools.intersection(set(known_genes)))
    overlap_GISTIC_known = list(set(gene_lists[0]).intersection(set(known_genes)))
    overlap_RUBIC_known = list(set(gene_lists[1]).intersection(set(known_genes)))
    print('Overlap RUBIC, GISTIC & known:\t' + ", ".join(gene_ID_to_name(overlap_all, ref_genome)))
    print('Overlap GISTIC & known:\t' + ", ".join(gene_ID_to_name(overlap_GISTIC_known, ref_genome)))
    print('Overlap RUBIC & known:\t' + ", ".join(gene_ID_to_name(overlap_RUBIC_known, ref_genome)))
    plt.subplots(figsize=(7, 7))
    venn_sets = [set(gene_lists[0]), set(gene_lists[1]), set(known_genes)]
    c = venn3(venn_sets, ('GISTIC', 'RUBIC', '                 Known genes'))
    try:
        c.get_patch_by_id('100').set_color('#5975A4'), c.get_patch_by_id('100').set_alpha(1)
        c.get_patch_by_id('010').set_color('#5F9E6E'), c.get_patch_by_id('010').set_alpha(1)
        c.get_patch_by_id('110').set_color('#5C888B'), c.get_patch_by_id('110').set_alpha(1)
        c.get_patch_by_id('001').set_color('#857AAA'), c.get_patch_by_id('001').set_alpha(1)
        c.get_patch_by_id('101').set_color('#6D77A7'), c.get_patch_by_id('101').set_alpha(1)
        c.get_patch_by_id('011').set_color('#748A8F'), c.get_patch_by_id('011').set_alpha(1)
        c.get_patch_by_id('111').set_color('#B55D60'), c.get_patch_by_id('111').set_alpha(1)
    except:
        pass    #too small overlap from three groups to show the overlapping region
    plt.savefig(out_venn, dpi=300)

def plot_histogram_sizes(parsed_results, plot_file):
    """Make boxplot with swarmplot on top showing sizes of recurrent regions detected by GISTIC2.0 and RUBIC."""
    size_list, size_list_known = parse_sizes_for_plot(parsed_results)
    sns.set_style("white")
    f, ax = plt.subplots(figsize=(7, 7))
    ax.set(yscale="log")
    g = sns.boxplot(x="Tool", y="Size", hue="CNV type", data=size_list, whis=np.inf, palette={"Amplification": "#71AEC0", "Deletion": "#B55D60"})
    g = sns.swarmplot(x="Tool", y="Size", dodge=True, hue="CNV type", data=size_list, size=5, palette={"Amplification": "#527F8C", "Deletion": "#844446"})
    g = sns.swarmplot(x="Tool", y="Size", dodge=True, hue="CNV type", data=size_list_known, size=5, palette={"Amplification": "black", "Deletion": "black"})
    sns.despine()
    g.set_xlabel("Tool",fontsize=16)
    g.set_ylabel("Log size of recurrent regions",fontsize=16)
    g.tick_params(labelsize=12)
    handles, labels = g.get_legend_handles_labels()         #Fix legend to prevent duplicates
    plt.legend(handles[0:2], labels[:2], fontsize=12)
    plt.savefig(plot_file, dpi=300)

def parse_sizes_for_plot(parsed_results):
    """Extract the sizes of regions to use as input for plot of sizes of recurrent regions."""
    sizes, sizes_known = [], []
    for tool_results in parsed_results:
        tool = 'GISTIC' if tool_results == parsed_results[0] else 'RUBIC'
        amp_sizes, del_sizes = [], []
        for cnv in tool_results:
            size = int(cnv[2]) - int(cnv[1])
            cnv_label = 'Amplification' if cnv[3] == 'amp' else 'Deletion'
            sizes += [(size, tool, cnv_label)]
            if cnv[5] != []:
                sizes_known += [(size, tool, cnv_label)]
    labels = ['Size', 'Tool', 'CNV type']
    df_sizes = pd.DataFrame(sizes, columns=labels)
    df_sizes_known = pd.DataFrame(sizes_known, columns=labels)
    return df_sizes, df_sizes_known

def make_bed_genes(known_file, biomart_file, out_bed, ref_genome):
    """Make bed file with location of known genes."""
    known_genes = parse_gene_file(known_file, np.inf, ref_genome)
    bed_lines = get_location_genes(known_genes, biomart_file)
    with open(out_bed, 'w') as out:
        out.write("track name=Known genes color=0,0,255\n")
        for line in bed_lines:
            out.write("\t".join(line) + "\n")

def get_location_genes(known_genes, biomart_file):
    """Extract locations of a list of genes from a biomart file with all genes to use for a bed file."""
    bed_lines = []
    with open(biomart_file, 'r') as biomart:
        biomart.readline()
        for line in biomart:
            gene = line.split("\t")[0]
            if gene in known_genes:
                gene_ID, gene_name, chrom, start, end = line.strip().split("\t")
                bed_lines.append(["chr" + chrom, start, end, gene_name])
    return bed_lines