Day 9: Visualization & Publication - Making Data Tell Stories

from pycirclize import Circos
import pandas as pd

# Initialize
circos = Circos(sectors={"chr1": 10, "chr2": 15, "chr3": 12})

# Add tracks
for sector in circos.sectors:
    track = sector.add_track((80, 100))
    track.axis()

    # GC content
    track.heatmap(gc_data)

    # Gene positions
    track.xticks_by_interval(1000000)

# Add links between genomes
circos.link(("chr1", 1000, 2000), ("chr2", 5000, 6000))

# Save
circos.savefig("genome_comparison.pdf")

import igraph as ig
import matplotlib.pyplot as plt
import pandas as pd

# Load co-occurrence network
edges = pd.read_csv('cooccurrence_network.csv')

# Create graph
g = ig.Graph.TupleList(edges.itertuples(index=False),
                       directed=False,
                       edge_attrs=['weight'])

# Calculate metrics
g.vs['degree'] = g.degree()
g.vs['betweenness'] = g.betweenness()

# Community detection
communities = g.community_multilevel()
g.vs['community'] = communities.membership

# Visual style
visual_style = {
    'vertex_size': [d*3 for d in g.vs['degree']],
    'vertex_color': g.vs['community'],
    'vertex_label': g.vs['name'],
    'edge_width': [w*2 for w in g.es['weight']],
    'layout': g.layout_fruchterman_reingold(),
    'bbox': (800, 800),
    'margin': 50
}

# Plot
ig.plot(g, 'network.pdf', **visual_style)

import networkx as nx
import matplotlib.pyplot as plt
import pandas as pd

# Create network
G = nx.Graph()

# Add edges from co-occurrence data
edges = pd.read_csv('cooccurrence_network.csv')
for _, row in edges.iterrows():
    G.add_edge(row['MAG1'], row['MAG2'], weight=row['correlation'])

# Remove weak connections
G_filtered = G.copy()
edges_to_remove = [(u, v) for u, v, d in G_filtered.edges(data=True)
                   if abs(d['weight']) < 0.3]
G_filtered.remove_edges_from(edges_to_remove)

# Layout
pos = nx.spring_layout(G_filtered, k=0.5, iterations=50)

# Node sizes by degree
node_sizes = [G_filtered.degree(node) * 100 for node in G_filtered.nodes()]

# Node colors by community
communities = nx.community.greedy_modularity_communities(G_filtered)
node_colors = []
for node in G_filtered.nodes():
    for i, comm in enumerate(communities):
        if node in comm:
            node_colors.append(i)
            break

# Plot
plt.figure(figsize=(14, 14))

# Draw edges
edges = G_filtered.edges()
weights = [G_filtered[u][v]['weight'] for u, v in edges]
nx.draw_networkx_edges(G_filtered, pos, alpha=0.3, width=weights)

# Draw nodes
nx.draw_networkx_nodes(G_filtered, pos,
                       node_size=node_sizes,
                       node_color=node_colors,
                       cmap=plt.cm.Set3,
                       alpha=0.8)

# Draw labels (only for large nodes)
large_nodes = [node for node in G_filtered.nodes()
               if G_filtered.degree(node) > 5]
labels = {node: node for node in large_nodes}
nx.draw_networkx_labels(G_filtered, pos, labels, font_size=8)

plt.title('MAG Co-occurrence Network', fontsize=16, fontweight='bold')
plt.axis('off')
plt.tight_layout()
plt.savefig('cooccurrence_network.pdf', dpi=300, bbox_inches='tight')
plt.show()

# Export for Cytoscape
import pandas as pd

# Nodes file
nodes = pd.DataFrame({
    'id': list(G.nodes()),
    'degree': [G.degree(n) for n in G.nodes()],
    'betweenness': [nx.betweenness_centrality(G)[n] for n in G.nodes()]
})
nodes.to_csv('network_nodes.csv', index=False)

# Edges file
edges_export = []
for u, v, d in G.edges(data=True):
    edges_export.append({
        'source': u,
        'target': v,
        'weight': d['weight']
    })
edges_df = pd.DataFrame(edges_export)
edges_df.to_csv('network_edges.csv', index=False)

print("Import these files into Cytoscape for advanced styling")

import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np
from scipy.cluster import hierarchy
from scipy.spatial.distance import pdist

# Load data
abundance = pd.read_csv('mag_abundance.csv', index_col=0)
metadata = pd.read_csv('metadata.csv', index_col=0)

# Log transform
abundance_log = np.log10(abundance + 1)

# Hierarchical clustering
row_linkage = hierarchy.linkage(pdist(abundance_log), method='ward')
col_linkage = hierarchy.linkage(pdist(abundance_log.T), method='ward')

# Create figure with GridSpec
fig = plt.figure(figsize=(14, 12))
gs = fig.add_gridspec(3, 3,
                      height_ratios=[0.3, 3, 0.2],
                      width_ratios=[0.2, 3, 0.3],
                      hspace=0.05, wspace=0.05)

# Top annotation (metadata)
ax_top = fig.add_subplot(gs[0, 1])
treatment_colors = {'Control': '#3498db', 'Treatment': '#e74c3c'}
treatment_bars = [treatment_colors[metadata.loc[s, 'Treatment']]
                  for s in abundance.columns]
ax_top.imshow([treatment_bars], aspect='auto')
ax_top.set_xticks([])
ax_top.set_yticks([0])
ax_top.set_yticklabels(['Treatment'])
ax_top.spines['top'].set_visible(False)
ax_top.spines['right'].set_visible(False)

# Left dendrogram
ax_left = fig.add_subplot(gs[1, 0])
hierarchy.dendrogram(row_linkage, orientation='left', ax=ax_left)
ax_left.set_xticks([])
ax_left.set_yticks([])
ax_left.spines['top'].set_visible(False)
ax_left.spines['right'].set_visible(False)
ax_left.spines['bottom'].set_visible(False)
ax_left.spines['left'].set_visible(False)

# Main heatmap
ax_main = fig.add_subplot(gs[1, 1])
im = ax_main.imshow(abundance_log, aspect='auto', cmap='viridis')
ax_main.set_xticks(range(len(abundance.columns)))
ax_main.set_xticklabels(abundance.columns, rotation=90)
ax_main.set_yticks(range(len(abundance.index)))
ax_main.set_yticklabels(abundance.index)

# Colorbar
ax_cbar = fig.add_subplot(gs[1, 2])
plt.colorbar(im, cax=ax_cbar, label='log10(Abundance + 1)')

plt.savefig('complex_heatmap.pdf', dpi=300, bbox_inches='tight')
plt.show()

from bokeh.plotting import figure, show, output_file
from bokeh.models import HoverTool, ColumnDataSource
from bokeh.layouts import column, row
import pandas as pd

# Load data
abundance = pd.read_csv('mag_abundance.csv')

# Create ColumnDataSource
source = ColumnDataSource(data=dict(
    x=abundance['Sample'],
    y=abundance['MAG_1'],
    desc=abundance['Taxonomy']
))

# Create figure
p = figure(x_range=abundance['Sample'].unique(),
           height=400, width=800,
           title="MAG Abundance Across Samples",
           toolbar_location="above")

# Add bars
p.vbar(x='x', top='y', width=0.8, source=source,
       color='navy', alpha=0.8)

# Add hover tool
hover = HoverTool(tooltips=[
    ("Sample", "@x"),
    ("Abundance", "@y"),
    ("Taxonomy", "@desc")
])
p.add_tools(hover)

# Styling
p.xaxis.major_label_orientation = 45
p.xgrid.grid_line_color = None

# Save
output_file("interactive_abundance.html")
show(p)

# Colorblind-friendly palettes
import matplotlib.pyplot as plt

# Qualitative (categorical data)
cb_friendly_qual = ['#E69F00', '#56B4E9', '#009E73', '#F0E442',
                    '#0072B2', '#D55E00', '#CC79A7']

# Sequential (continuous data)
cb_friendly_seq = plt.cm.viridis  # or 'plasma', 'cividis'

# Diverging (centered data)
cb_friendly_div = plt.cm.RdBu_r

# Use in plots
plt.scatter(x, y, c=categories, cmap=cb_friendly_qual)

import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec

# Create figure
fig = plt.figure(figsize=(16, 10))
gs = GridSpec(2, 3, figure=fig, hspace=0.3, wspace=0.3)

# Panel A: Taxonomy barplot
ax1 = fig.add_subplot(gs[0, :2])
# ... plot code ...
ax1.text(-0.1, 1.05, 'A', transform=ax1.transAxes,
         fontsize=16, fontweight='bold', va='top')

# Panel B: Alpha diversity
ax2 = fig.add_subplot(gs[0, 2])
# ... plot code ...
ax2.text(-0.15, 1.05, 'B', transform=ax2.transAxes,
         fontsize=16, fontweight='bold', va='top')

# Panel C: Ordination
ax3 = fig.add_subplot(gs[1, 0])
# ... plot code ...
ax3.text(-0.15, 1.05, 'C', transform=ax3.transAxes,
         fontsize=16, fontweight='bold', va='top')

# Panel D: Heatmap
ax4 = fig.add_subplot(gs[1, 1:])
# ... plot code ...
ax4.text(-0.08, 1.05, 'D', transform=ax4.transAxes,
         fontsize=16, fontweight='bold', va='top')

plt.savefig('figure1_composite.pdf', dpi=300, bbox_inches='tight')
plt.show()

# Always clean data first
def prepare_data_for_plotting(df):
    # Remove zeros
    df = df[(df > 0).any(axis=1)]

    # Log transform
    df_log = np.log10(df + 1)

    # Normalize if needed
    df_norm = df.div(df.sum(axis=0), axis=1) * 100

    return df_log, df_norm

# Set global matplotlib style
import matplotlib.pyplot as plt

plt.rcParams.update({
    'font.size': 12,
    'font.family': 'sans-serif',
    'font.sans-serif': ['Arial'],
    'axes.linewidth': 1.5,
    'axes.labelsize': 12,
    'axes.titlesize': 14,
    'xtick.labelsize': 10,
    'ytick.labelsize': 10,
    'legend.fontsize': 10,
    'figure.dpi': 100,
    'savefig.dpi': 300,
    'savefig.bbox': 'tight'
})

def save_figure(fig, name):
    """Save figure in multiple formats"""
    fig.savefig(f'{name}.pdf', dpi=300, bbox_inches='tight')
    fig.savefig(f'{name}.png', dpi=300, bbox_inches='tight')
    fig.savefig(f'{name}.svg', bbox_inches='tight')  # For editing
    print(f"✓ Saved {name} in PDF, PNG, and SVG formats")