ubuntu2404
=>PYTHONTEX#py#default#default#0#code#####62#
# Import core bioinformatics and data analysis libraries
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use('Agg')
import seaborn as sns
import pandas as pd
from scipy import stats
from scipy.cluster.hierarchy import dendrogram, linkage
from scipy.spatial.distance import squareform
# BioPython imports for sequence analysis
try:
from Bio import SeqIO, Align, Phylo
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.SeqUtils import GC, molecular_weight
from Bio.SeqUtils.ProtParam import ProteinAnalysis
from Bio.Blast import NCBIXML
from Bio import Entrez
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor
from Bio.Align import MultipleSeqAlignment
biopython_available = True
except ImportError:
print("BioPython not available - using simulated data")
biopython_available = False
# Set visualization parameters for genomic data
plt.style.use('seaborn-v0_8-whitegrid')
np.random.seed(42)
sns.set_palette("Set2")
# Genomics-optimized figure settings
plt.rcParams['figure.figsize'] = (10, 6)
plt.rcParams['figure.dpi'] = 150
plt.rcParams['savefig.bbox'] = 'tight'
plt.rcParams['savefig.pad_inches'] = 0.1
plt.rcParams['font.family'] = 'sans-serif'
# Configure pandas display to avoid overfull boxes
pd.set_option('display.max_columns', 8)
pd.set_option('display.width', 80)
pd.set_option('display.precision', 3)
# Define nucleotide and amino acid color schemes
nucleotide_colors = {'A': '#FF6B6B', 'T': '#4ECDC4', 'G': '#45B7D1', 'C': '#96CEB4'}
amino_acid_colors = {
'A': '#FF6B6B', 'R': '#4ECDC4', 'N': '#45B7D1', 'D': '#96CEB4',
'C': '#FECA57', 'Q': '#48CAE4', 'E': '#F38BA8', 'G': '#A8DADC',
'H': '#FFB3BA', 'I': '#BFEFFF', 'L': '#FFDFBA', 'K': '#FFFFBA',
'M': '#BAE1FF', 'F': '#FFBAE1', 'P': '#E1BAFF', 'S': '#D4E2FC',
'T': '#FCE4D4', 'W': '#E4FCD4', 'Y': '#D4FCE4', 'V': '#FCD4E4'
}
=>PYTHONTEX#py#default#default#1#code#####192#
# Generate synthetic E. coli genome data for demonstration
# Real analysis would use: Entrez.efetch() from NCBI
def generate_synthetic_sequence(length, gc_content=0.51):
"""Generate synthetic DNA sequence with specified GC content"""
np.random.seed(42)
# Calculate nucleotide probabilities based on GC content
gc_prob = gc_content / 2 # Equal probability for G and C
at_prob = (1 - gc_content) / 2 # Equal probability for A and T
nucleotides = ['A', 'T', 'G', 'C']
probabilities = [at_prob, at_prob, gc_prob, gc_prob]
sequence = np.random.choice(nucleotides, size=length, p=probabilities)
return ''.join(sequence)
# Generate synthetic E. coli strain sequences
strain_names = [
'EcoliK12', 'EcoliO157H7', 'EcoliCFT073', 'EcoliUTI89',
'EcoliEDL933', 'EcoliMG1655', 'EcoliDH10B', 'EcoliBL21'
]
# Simulate genomic regions (e.g., 16S rRNA gene sequences)
sequence_length = 1500 # Typical 16S rRNA length
genomic_sequences = {}
for strain in strain_names:
# Vary GC content slightly between strains (E. coli ~50-52%)
gc_content = 0.50 + np.random.normal(0, 0.01)
gc_content = max(0.48, min(0.54, gc_content)) # Bound within realistic range
sequence = generate_synthetic_sequence(sequence_length, gc_content)
genomic_sequences[strain] = sequence
print(f"Generated synthetic genomic sequences for {len(strain_names)} E. coli strains")
print(f"Sequence length: {sequence_length} bp")
# Calculate basic sequence statistics
sequence_stats = {}
for strain, sequence in genomic_sequences.items():
gc_content = (sequence.count('G') + sequence.count('C')) / len(sequence)
at_content = (sequence.count('A') + sequence.count('T')) / len(sequence)
sequence_stats[strain] = {
'Length': len(sequence),
'GCcontent': gc_content,
'ATcontent': at_content,
'Acount': sequence.count('A'),
'Tcount': sequence.count('T'),
'Gcount': sequence.count('G'),
'Ccount': sequence.count('C')
}
# Convert to DataFrame for analysis
stats_df = pd.DataFrame.from_dict(sequence_stats, orient='index')
print("\nSequence composition statistics:")
# Print in a compact format without problematic characters
print(stats_df.round(4).to_string(max_cols=6, max_colwidth=12))
=>PYTHONTEX#py#default#default#2#code#####258#
# Calculate pairwise sequence distances
def hamming_distance(seq1, seq2):
"""Calculate Hamming distance between two sequences"""
if len(seq1) != len(seq2):
raise ValueError("Sequences must be of equal length")
return sum(c1 != c2 for c1, c2 in zip(seq1, seq2))
# Create distance matrix
strains = list(genomic_sequences.keys())
n_strains = len(strains)
distance_matrix = np.zeros((n_strains, n_strains))
for i, strain1 in enumerate(strains):
for j, strain2 in enumerate(strains):
if i != j:
hamming_dist = hamming_distance(genomic_sequences[strain1],
genomic_sequences[strain2])
# Convert to evolutionary distance (proportion of differences)
evolutionary_distance = hamming_dist / sequence_length
distance_matrix[i, j] = evolutionary_distance
# Convert to DataFrame for better presentation
distance_df = pd.DataFrame(distance_matrix, index=strains, columns=strains)
print("Pairwise evolutionary distances (proportion of differences):")
# Print in a compact format to avoid overfull boxes
print(distance_df.round(4).to_string(max_cols=8, max_colwidth=10))
# Calculate mean distances for each strain
mean_distances = distance_df.mean(axis=1)
print(f"\nMean evolutionary distances:")
for strain, dist in mean_distances.items():
print(f"{strain}: {dist:.4f}")
=>PYTHONTEX#py#default#default#3#code#####299#
# Phylogenetic reconstruction using hierarchical clustering
from scipy.cluster.hierarchy import dendrogram, linkage
from scipy.spatial.distance import squareform
# Convert distance matrix to condensed form for clustering
condensed_distances = squareform(distance_matrix)
# Perform hierarchical clustering (UPGMA method)
linkage_matrix = linkage(condensed_distances, method='average')
print("Phylogenetic reconstruction completed using UPGMA method")
print(f"Linkage matrix shape: {linkage_matrix.shape}")
# Extract clustering information
cluster_info = []
for i, merge in enumerate(linkage_matrix):
cluster_info.append({
'step': i + 1,
'cluster1': int(merge[0]) if merge[0] < len(strains) else f"Cluster{int(merge[0])}",
'cluster2': int(merge[1]) if merge[1] < len(strains) else f"Cluster{int(merge[1])}",
'distance': merge[2],
'size': int(merge[3])
})
cluster_df = pd.DataFrame(cluster_info)
print("\nClustering steps:")
print(cluster_df.round(4).to_string(max_cols=6, max_colwidth=12))
=>PYTHONTEX#py#default#default#4#code#####340#
import os
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
# Ensure the figures directory exists
os.makedirs('figures', exist_ok=True)
# Create comprehensive sequence composition visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.suptitle('E. coli Strain Genomic Sequence Composition Analysis',
fontsize=16, fontweight='bold')
# GC content distribution
ax1 = axes[0, 0]
gc_values = stats_df['GCcontent'].values
ax1.hist(gc_values, bins=8, alpha=0.7, color='skyblue', edgecolor='black')
ax1.axvline(gc_values.mean(), color='red', linestyle='--', linewidth=2,
label=f'Mean: {gc_values.mean():.3f}')
ax1.set_xlabel('GC Content')
ax1.set_ylabel('Number of Strains')
ax1.set_title('GC Content Distribution')
ax1.legend()
ax1.grid(True, alpha=0.3)
# Nucleotide composition heatmap
ax2 = axes[0, 1]
nucleotide_data = stats_df[['Acount', 'Tcount', 'Gcount', 'Ccount']].T
sns.heatmap(nucleotide_data, annot=True, fmt='d', cmap='YlOrRd', ax=ax2,
xticklabels=[s.replace('Ecoli', '') for s in strains],
yticklabels=['A', 'T', 'G', 'C'])
ax2.set_title('Nucleotide Counts by Strain')
ax2.set_xlabel('E. coli Strains')
# GC vs AT content scatter plot
ax3 = axes[1, 0]
scatter = ax3.scatter(stats_df['GCcontent'], stats_df['ATcontent'],
s=80, alpha=0.7, c=range(len(strains)), cmap='viridis')
ax3.set_xlabel('GC Content')
ax3.set_ylabel('AT Content')
ax3.set_title('GC vs AT Content Relationship')
# Add strain labels
for i, strain in enumerate(strains):
ax3.annotate(strain.replace('Ecoli', ''),
(stats_df.loc[strain, 'GCcontent'],
stats_df.loc[strain, 'ATcontent']),
xytext=(5, 5), textcoords='offset points', fontsize=8)
ax3.grid(True, alpha=0.3)
# Sequence diversity analysis
ax4 = axes[1, 1]
diversity_values = [distance_df.loc[strain].sum() for strain in strains]
bars = ax4.bar(range(len(strains)), diversity_values,
color=plt.cm.Set3(np.linspace(0, 1, len(strains))))
ax4.set_xlabel('E. coli Strains')
ax4.set_ylabel('Total Evolutionary Distance')
ax4.set_title('Genomic Diversity Index')
ax4.set_xticks(range(len(strains)))
ax4.set_xticklabels([s.replace('Ecoli', '') for s in strains], rotation=45)
ax4.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('figures/sequence_composition.pdf', dpi=300, bbox_inches='tight')
plt.close()
print("Figure saved to figures/sequence composition.pdf")
=>PYTHONTEX#py#default#default#5#code#####425#
# Create phylogenetic analysis visualization
fig, axes = plt.subplots(2, 2, figsize=(15, 12))
fig.suptitle('Phylogenetic Analysis of E. coli Strain Relationships',
fontsize=16, fontweight='bold')
# Dendrogram (phylogenetic tree)
ax1 = axes[0, 0]
dendrogram_result = dendrogram(linkage_matrix,
labels=[s.replace('Ecoli', '') for s in strains],
ax=ax1, leaf_rotation=90, leaf_font_size=10)
ax1.set_title('UPGMA Phylogenetic Tree')
ax1.set_xlabel('E. coli Strains')
ax1.set_ylabel('Evolutionary Distance')
# Distance matrix heatmap
ax2 = axes[0, 1]
mask = np.triu(np.ones_like(distance_matrix, dtype=bool)) # Mask upper triangle
sns.heatmap(distance_df, mask=mask, annot=True, fmt='.3f',
cmap='RdYlBu_r', ax=ax2,
xticklabels=[s.replace('Ecoli', '') for s in strains],
yticklabels=[s.replace('E_coli_', '') for s in strains])
ax2.set_title('Pairwise Evolutionary Distance Matrix')
# Clustering dendrogram (horizontal)
ax3 = axes[1, 0]
dendrogram(linkage_matrix,
labels=[s.replace('Ecoli', '') for s in strains],
ax=ax3, orientation='left', leaf_font_size=10)
ax3.set_title('Horizontal Dendrogram View')
ax3.set_xlabel('Evolutionary Distance')
# Distance distribution analysis
ax4 = axes[1, 1]
# Extract all pairwise distances (excluding diagonal)
all_distances = distance_matrix[np.triu_indices_from(distance_matrix, k=1)]
ax4.hist(all_distances, bins=15, alpha=0.7, color='lightcoral',
edgecolor='black', density=True)
ax4.axvline(all_distances.mean(), color='red', linestyle='--', linewidth=2,
label=f'Mean: {all_distances.mean():.4f}')
ax4.axvline(np.median(all_distances), color='blue', linestyle='--', linewidth=2,
label=f'Median: {np.median(all_distances):.4f}')
ax4.set_xlabel('Evolutionary Distance')
ax4.set_ylabel('Density')
ax4.set_title('Distribution of Pairwise Distances')
ax4.legend()
ax4.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('figures/phylogenetic_analysis.pdf', dpi=300, bbox_inches='tight')
plt.close()
# Print summary statistics
print(f"\nPhylogenetic Analysis Summary:")
print(f"Number of strains analyzed: {len(strains)}")
print(f"Mean pairwise distance: {all_distances.mean():.4f}")
print(f"Standard deviation: {all_distances.std():.4f}")
print(f"Minimum distance: {all_distances.min():.4f}")
print(f"Maximum distance: {all_distances.max():.4f}")
=>PYTHONTEX#py#default#default#6#code#####498#
# Simulate gene content and functional analysis
np.random.seed(42)
# Generate synthetic gene presence/absence data
gene_families = [
'CoreMetabolism', 'DNArepair', 'CellWall', 'Transport', 'Regulation',
'Pathogenicity', 'AntibioticResistance', 'MobileElements', 'StressResponse',
'SecretionSystems', 'Chemotaxis', 'FlagellarBiosynthesis'
]
# Simulate gene presence (1) or absence (0) for each strain and gene family
gene_matrix = np.random.binomial(1, 0.7, size=(len(strains), len(gene_families)))
# Ensure core genes are present in all strains
core_genes = ['CoreMetabolism', 'DNArepair', 'CellWall']
for i, gene_family in enumerate(gene_families):
if gene_family in core_genes:
gene_matrix[:, i] = 1
# Create gene content DataFrame
gene_content_df = pd.DataFrame(gene_matrix,
index=strains,
columns=gene_families)
print("Gene family presence/absence matrix:")
print(gene_content_df.to_string(max_cols=8, max_colwidth=15))
# Calculate gene family conservation
conservation_scores = gene_content_df.sum(axis=0) / len(strains)
print(f"\nGene family conservation scores:")
for gene_family, score in conservation_scores.items():
print(f"{gene_family}: {score:.2f}")
# Calculate strain-specific gene diversity
strain_diversity = gene_content_df.sum(axis=1)
print(f"\nStrain gene family diversity:")
for strain, diversity in strain_diversity.items():
print(f"{strain}: {diversity} gene families")
=>PYTHONTEX#py#default#default#7#code#####543#
# Create gene content and functional analysis visualization
fig, axes = plt.subplots(2, 2, figsize=(15, 10))
fig.suptitle('Gene Content and Functional Diversity Analysis',
fontsize=16, fontweight='bold')
# Gene presence/absence heatmap
ax1 = axes[0, 0]
sns.heatmap(gene_content_df.T, cmap='RdYlGn', cbar_kws={'label': 'Gene Presence'},
xticklabels=[s.replace('Ecoli', '') for s in strains],
yticklabels=gene_families, ax=ax1)
ax1.set_title('Gene Family Presence/Absence Matrix')
ax1.set_xlabel('E. coli Strains')
ax1.set_ylabel('Gene Families')
# Gene family conservation bar plot
ax2 = axes[0, 1]
bars = ax2.bar(range(len(gene_families)), conservation_scores.values,
color=plt.cm.viridis(conservation_scores.values))
ax2.set_xlabel('Gene Families')
ax2.set_ylabel('Conservation Score')
ax2.set_title('Gene Family Conservation Across Strains')
ax2.set_xticks(range(len(gene_families)))
ax2.set_xticklabels(gene_families, rotation=45, ha='right')
ax2.grid(True, alpha=0.3)
# Add conservation threshold line
ax2.axhline(y=0.8, color='red', linestyle='--', alpha=0.7,
label='High conservation threshold')
ax2.legend()
# Strain gene diversity
ax3 = axes[1, 0]
diversity_colors = plt.cm.Set3(np.linspace(0, 1, len(strains)))
bars = ax3.bar(range(len(strains)), strain_diversity.values, color=diversity_colors)
ax3.set_xlabel('E. coli Strains')
ax3.set_ylabel('Number of Gene Families')
ax3.set_title('Gene Family Diversity by Strain')
ax3.set_xticks(range(len(strains)))
ax3.set_xticklabels([s.replace('E_coli_', '') for s in strains], rotation=45)
ax3.grid(True, alpha=0.3)
# Gene content clustering
ax4 = axes[1, 1]
# Calculate Jaccard distances for gene content
from scipy.spatial.distance import pdist, squareform
from scipy.cluster.hierarchy import dendrogram, linkage
gene_distances = pdist(gene_matrix, metric='jaccard')
gene_linkage = linkage(gene_distances, method='average')
dendrogram(gene_linkage,
labels=[s.replace('Ecoli', '') for s in strains],
ax=ax4, leaf_rotation=90, leaf_font_size=10)
ax4.set_title('Gene Content-Based Clustering')
ax4.set_xlabel('E. coli Strains')
ax4.set_ylabel('Jaccard Distance')
plt.tight_layout()
plt.savefig('figures/gene_content_analysis.pdf', dpi=300, bbox_inches='tight')
plt.close()
=>PYTHONTEX#py#default#default#8#i#####621#
f"{all_distances.min():.4f}"
=>PYTHONTEX#py#default#default#9#i#####621#
f"{all_distances.max():.4f}"
=>PYTHONTEX#py#default#default#10#i#####634#
f"{strain_diversity.min()}"
=>PYTHONTEX#py#default#default#11#i#####634#
f"{strain_diversity.max()}"
=>PYTHONTEX:SETTINGS#
version=0.18
outputdir=pythontex-files-main
workingdir=.
workingdirset=false
gobble=none
rerun=default
hashdependencies=default
makestderr=false
stderrfilename=full
keeptemps=none
pyfuture=default
pyconfuture=none
pygments=true
pygglobal=:GLOBAL||
fvextfile=-1
pyconbanner=none
pyconfilename=stdin
depythontex=false
pygfamily=py|python3|
pygfamily=pycon|pycon|
pygfamily=sympy|python3|
pygfamily=sympycon|pycon|
pygfamily=pylab|python3|
pygfamily=pylabcon|pycon|