Turn into right format from camber cloud jupyter notebook format
This commit is contained in:
-318
File diff suppressed because one or more lines are too long
@@ -0,0 +1,221 @@
|
|||||||
|
import numpy as np
|
||||||
|
import matplotlib.pyplot as plt
|
||||||
|
import pandas as pd
|
||||||
|
from scipy import ndimage
|
||||||
|
|
||||||
|
matrix = np.loadtxt('chimpanzee_v_human.txt', skiprows=2)
|
||||||
|
print(f"Matrix shape: {matrix.shape}")
|
||||||
|
|
||||||
|
def apply_screening_window(matrix, window_size=5, threshold=0):
|
||||||
|
|
||||||
|
# Option 1: Moving average filter
|
||||||
|
filtered = ndimage.uniform_filter(matrix, size=window_size)
|
||||||
|
|
||||||
|
# Option 2: Apply threshold to remove low values
|
||||||
|
filtered[filtered < threshold] = 0
|
||||||
|
|
||||||
|
# Option 3: Downsample to reduce size
|
||||||
|
# Take every Nth point
|
||||||
|
downsampled = filtered[::2, ::2] # Takes every 2nd point
|
||||||
|
|
||||||
|
return filtered
|
||||||
|
|
||||||
|
# Apply filter
|
||||||
|
filtered_matrix = apply_screening_window(matrix, window_size=5, threshold=0.1)
|
||||||
|
|
||||||
|
def get_chromosome_boundaries(chr_sizes, matrix_dimension, genome_total_size):
|
||||||
|
"""
|
||||||
|
Convert chromosome sizes to matrix positions
|
||||||
|
|
||||||
|
Parameters:
|
||||||
|
-----------
|
||||||
|
chr_sizes : dict
|
||||||
|
Dictionary of chromosome names and their sizes in base pairs
|
||||||
|
matrix_dimension : int
|
||||||
|
The size of the matrix dimension (e.g., 2000 for x-axis or y-axis)
|
||||||
|
genome_total_size : int
|
||||||
|
Total genome size in base pairs (sum of all chromosome sizes)
|
||||||
|
|
||||||
|
Returns:
|
||||||
|
--------
|
||||||
|
boundaries : list
|
||||||
|
Matrix positions where chromosomes start
|
||||||
|
labels : list
|
||||||
|
Chromosome names
|
||||||
|
midpoints : list
|
||||||
|
Midpoints for placing chromosome labels
|
||||||
|
"""
|
||||||
|
cumulative_bp = 0
|
||||||
|
boundaries = [0] # Start at 0
|
||||||
|
labels = []
|
||||||
|
midpoints = []
|
||||||
|
|
||||||
|
for chr_name, size_bp in chr_sizes.items():
|
||||||
|
# Calculate the proportion of the genome this chromosome represents
|
||||||
|
proportion = size_bp / genome_total_size
|
||||||
|
|
||||||
|
# Convert to matrix position
|
||||||
|
matrix_size_for_chr = proportion * matrix_dimension
|
||||||
|
|
||||||
|
# Update cumulative position
|
||||||
|
cumulative_bp += size_bp
|
||||||
|
cumulative_matrix_pos = int((cumulative_bp / genome_total_size) * matrix_dimension)
|
||||||
|
|
||||||
|
# Calculate midpoint for label
|
||||||
|
midpoint = (boundaries[-1] + cumulative_matrix_pos) / 2
|
||||||
|
midpoints.append(midpoint)
|
||||||
|
|
||||||
|
# Store
|
||||||
|
boundaries.append(cumulative_matrix_pos)
|
||||||
|
labels.append(chr_name)
|
||||||
|
|
||||||
|
# Remove the last boundary (it's just the end of the matrix)
|
||||||
|
boundaries = boundaries[:-1]
|
||||||
|
|
||||||
|
return boundaries, labels, midpoints
|
||||||
|
|
||||||
|
# Example usage for human and chimpanzee
|
||||||
|
# Human chromosome sizes (hg38)
|
||||||
|
chr_sizes_human = {
|
||||||
|
'chr1': 248956422,
|
||||||
|
'chr2': 242193529,
|
||||||
|
'chr3': 198295559,
|
||||||
|
'chr4': 190214555,
|
||||||
|
'chr5': 181538259,
|
||||||
|
'chr6': 170805979,
|
||||||
|
'chr7': 159345973,
|
||||||
|
'chr8': 145138636,
|
||||||
|
'chr9': 138394717,
|
||||||
|
'chr10': 133797422,
|
||||||
|
'chr11': 135086622,
|
||||||
|
'chr12': 133275309,
|
||||||
|
'chr13': 114364328,
|
||||||
|
'chr14': 107043718,
|
||||||
|
'chr15': 101991189,
|
||||||
|
'chr16': 90338345,
|
||||||
|
'chr17': 83257441,
|
||||||
|
'chr18': 80373285,
|
||||||
|
'chr19': 58617616,
|
||||||
|
'chr20': 64444167,
|
||||||
|
'chr21': 46709983,
|
||||||
|
'chr22': 50818468,
|
||||||
|
'chrX': 156040895,
|
||||||
|
'chrY': 57227415,
|
||||||
|
}
|
||||||
|
|
||||||
|
# Chimpanzee chromosome sizes (panTro6)
|
||||||
|
chr_sizes_chimp = {
|
||||||
|
'chr1': 231299448,
|
||||||
|
'chr2': 201579334, # Note: chimp has 2A and 2B instead of fused chr2
|
||||||
|
'chr3': 191670129,
|
||||||
|
'chr4': 178964043,
|
||||||
|
'chr5': 183320697,
|
||||||
|
'chr6': 174648847,
|
||||||
|
'chr7': 155691389,
|
||||||
|
'chr8': 146288486,
|
||||||
|
'chr9': 136991367,
|
||||||
|
'chr10': 140160335,
|
||||||
|
'chr11': 118863153,
|
||||||
|
'chr12': 122853238,
|
||||||
|
'chr13': 144476839,
|
||||||
|
'chr14': 120149342,
|
||||||
|
'chr15': 104701121,
|
||||||
|
'chr16': 93846347,
|
||||||
|
'chr17': 90695042,
|
||||||
|
'chr18': 94528409,
|
||||||
|
'chr19': 99595324,
|
||||||
|
'chr20': 68087067,
|
||||||
|
'chr21': 78115444,
|
||||||
|
'chr22': 51350302,
|
||||||
|
'chr23': 59520411,
|
||||||
|
'chrX': 154259566,
|
||||||
|
'chrY': 26350395,
|
||||||
|
}
|
||||||
|
|
||||||
|
# Calculate total genome sizes
|
||||||
|
human_total_size = sum(chr_sizes_human.values())
|
||||||
|
chimp_total_size = sum(chr_sizes_chimp.values())
|
||||||
|
|
||||||
|
print(f"Human genome size: {human_total_size:,} bp")
|
||||||
|
print(f"Chimp genome size: {chimp_total_size:,} bp")
|
||||||
|
|
||||||
|
# Get boundaries for the 2000x2000 matrix
|
||||||
|
# X-axis: human (reference), Y-axis: chimp (query)
|
||||||
|
human_boundaries, human_labels, human_midpoints = get_chromosome_boundaries(
|
||||||
|
chr_sizes_human,
|
||||||
|
matrix_dimension=2000, # x-axis dimension
|
||||||
|
genome_total_size=human_total_size
|
||||||
|
)
|
||||||
|
|
||||||
|
chimp_boundaries, chimp_labels, chimp_midpoints = get_chromosome_boundaries(
|
||||||
|
chr_sizes_chimp,
|
||||||
|
matrix_dimension=2000, # y-axis dimension
|
||||||
|
genome_total_size=chimp_total_size
|
||||||
|
)
|
||||||
|
|
||||||
|
def apply_threshold(n):
|
||||||
|
if n > 30:
|
||||||
|
return 1
|
||||||
|
else:
|
||||||
|
return 0
|
||||||
|
|
||||||
|
|
||||||
|
def plot_dotplot_with_chromosomes(matrix,
|
||||||
|
ref_boundaries, ref_labels, ref_midpoints,
|
||||||
|
query_boundaries, query_labels, query_midpoints,
|
||||||
|
figsize=(14, 14), cmap='binary'):
|
||||||
|
apply_threshold_vectorized = np.vectorize(apply_threshold)
|
||||||
|
"""
|
||||||
|
Create dot plot with chromosome boundary lines using actual genome proportions
|
||||||
|
"""
|
||||||
|
fig, ax = plt.subplots(figsize=figsize)
|
||||||
|
|
||||||
|
# Plot the matrix
|
||||||
|
im = ax.imshow(apply_threshold_vectorized(matrix), cmap='binary', aspect='auto',
|
||||||
|
origin='lower', interpolation='nearest')
|
||||||
|
|
||||||
|
# Add colorbar
|
||||||
|
plt.colorbar(im, ax=ax, label='Similarity Score', shrink=0.8)
|
||||||
|
|
||||||
|
plt.style.use("default")
|
||||||
|
|
||||||
|
# Add vertical lines for reference genome chromosomes (x-axis)
|
||||||
|
for boundary in ref_boundaries[1:]: # Skip first boundary (0)
|
||||||
|
ax.axvline(x=boundary, color='cyan', linestyle='--',
|
||||||
|
linewidth=0.8, alpha=0.6)
|
||||||
|
|
||||||
|
# Add horizontal lines for query genome chromosomes (y-axis)
|
||||||
|
for boundary in query_boundaries[1:]: # Skip first boundary (0)
|
||||||
|
ax.axhline(y=boundary, color='cyan', linestyle='--',
|
||||||
|
linewidth=0.8, alpha=0.6)
|
||||||
|
|
||||||
|
# Add chromosome labels at midpoints
|
||||||
|
ax.set_xticks(ref_midpoints)
|
||||||
|
ax.set_xticklabels(ref_labels, rotation=45, ha='right', fontsize=9)
|
||||||
|
|
||||||
|
ax.set_yticks(query_midpoints)
|
||||||
|
ax.set_yticklabels(query_labels, fontsize=9)
|
||||||
|
|
||||||
|
ax.set_xlabel('Human Genome (Reference)', fontsize=12, fontweight='bold')
|
||||||
|
ax.set_ylabel('Chimpanzee Genome (Query)', fontsize=12, fontweight='bold')
|
||||||
|
ax.set_title('Human vs Chimpanzee Genome-wide Dot Plot',
|
||||||
|
fontsize=14, fontweight='bold', pad=20)
|
||||||
|
|
||||||
|
# Add grid for better readability
|
||||||
|
ax.grid(False)
|
||||||
|
|
||||||
|
plt.tight_layout()
|
||||||
|
return fig, ax
|
||||||
|
|
||||||
|
# Create the plot
|
||||||
|
fig, ax = plot_dotplot_with_chromosomes(
|
||||||
|
matrix,
|
||||||
|
human_boundaries, human_labels, human_midpoints,
|
||||||
|
chimp_boundaries, chimp_labels, chimp_midpoints,
|
||||||
|
figsize=(16, 16),
|
||||||
|
cmap='hot' # or 'RdYlBu_r', 'viridis', 'plasma'
|
||||||
|
)
|
||||||
|
|
||||||
|
plt.savefig('human_chimp_dotplot.png', dpi=300, bbox_inches='tight')
|
||||||
|
plt.show()
|
||||||
|
|
||||||
Reference in New Issue
Block a user