PCA for gene count comparison between NCBI Crassostrea gigas Annotation Release 101 and GCF_963853765.1-RS_2024_06 reference genomes

import pandas as pd
import matplotlib.pyplot as plt
from matplotlib_venn import venn2
import matplotlib
import numpy as np

# Use non-interactive backend for headless environments
matplotlib.use('Agg')

# Function to load gene count files (ignoring the second column)
def load_gene_counts(file_path):
    df = pd.read_csv(file_path, sep='\t')
    df = df.drop(columns=df.columns[1])  # Drop the second column (gene_name)
    return df

# Function to prepare the comparison table
def compare_gene_counts(df1, df2):
    # Merge data on gene_id
    # The suffixes modifies the column names for the respective dataframes
    merged = pd.merge(df1, df2, on='gene_id', how='outer', suffixes=('_run1', '_run2'))

    # Mark presence in each run
    # : is the row indexer; it means select all rows
    # 1 selects the second column
    # len(df1.columns selects the first column in df2 data
    # The 'present_in' columns are appended to each row
    # Only one column of each original dataframe is checked because the outer
    # join assigns NaN to all colums if the joined-on key isn't present.
    merged['present_in_run1'] = ~merged.iloc[:, 1].isna()
    merged['present_in_run2'] = ~merged.iloc[:, len(df1.columns)].isna()

    return merged

# Function to generate visualizations
def generate_visualizations(merged_df, output_prefix, df1_columns):
    # Venn Diagram
    only_in_run1 = merged_df['present_in_run1'] & ~merged_df['present_in_run2']
    only_in_run2 = merged_df['present_in_run2'] & ~merged_df['present_in_run1']
    in_both_runs = merged_df['present_in_run1'] & merged_df['present_in_run2']

    venn2(
        subsets=(only_in_run1.sum(), only_in_run2.sum(), in_both_runs.sum()),
        set_labels=('Run 1', 'Run 2')
    )
    plt.title("Gene ID Overlap Between Runs")
    plt.savefig(f"{output_prefix}_venn.png")
    plt.clf()

    # Scatter Plot for common genes
    common_genes = merged_df[merged_df['present_in_run1'] & merged_df['present_in_run2']]

    plt.xlabel("Log2 Run 1 Counts")
    plt.ylabel("Log2 Run 2 Counts")
    plt.title(f"Scatter Plot of Log2 Gene Counts for Common Genes")
    plt.grid(True)
    for col_idx in range(1, len(df1_columns)):
        column_name = df1_columns[col_idx]
        print(f"Processing sample: {column_name}")
        sample_1_run1 = column_name + '_run1'
        sample_1_run2 = column_name + '_run2'

        # Check if columns exist before plotting
        if sample_1_run1 in common_genes.columns and sample_1_run2 in common_genes.columns:
            x = np.log2(common_genes[sample_1_run1] + 1)
            y = np.log2(common_genes[sample_1_run2] + 1)

            plt.scatter(x, y, alpha=0.3, s=5)  # Reduced dot size (s=5)
        else:
            print(f"Error: Columns {sample_1_run1} or {sample_1_run2} not found in the merged dataset.")
    plt.savefig(f"{output_prefix}_scatter.png")

# Main Function
def main(file1, file2, output_prefix):
    # Load data
    df1 = load_gene_counts(file1)
    df2 = load_gene_counts(file2)

    # Perform comparison
    merged_df = compare_gene_counts(df1, df2)

    # Save tabular comparison
    merged_df.to_csv(f"{output_prefix}_comparison.csv", index=False)

    # Generate visualizations
    generate_visualizations(merged_df, output_prefix, df1.columns)

# Example usage:
# main("run1_counts.tsv", "run2_counts.tsv", "gene_comparison")


# Example usage:
main("ncbi_dataset_GCA_000297895.1/ncbi_dataset/data/GCF_000297895.1/salmon.merged.gene_counts.tsv", "seqera_roberto/salmon.merged.gene_counts_seqera.tsv", "gene_comparison")

# Import required libraries
import pandas as pd
import numpy as np
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import plotly.graph_objects as go

matplotlib.use('Agg')

destdir = '/mnt/c/Users/inter/Downloads/GMGI/ChatGPT/'
# Load the files
data1_file = f'{destdir}salmon.merged.gene_counts_2023_seqera.tsv'
data2_file = f'{destdir}salmon.merged.gene_counts_roberto_seqera.tsv'
metadata_file = f'{destdir}deseq_metadata2.csv'

# Load the data1 and data2, ignoring the second column
data1 = pd.read_csv(data1_file, sep='\t', index_col=0).iloc[:, 1:]
data2 = pd.read_csv(data2_file, sep='\t', index_col=0).iloc[:, 1:]

# Load the metadata
metadata = pd.read_csv(metadata_file)

# Find common genes between data1 and data2
common_genes = data1.index.intersection(data2.index)

# Subset both datasets to include only these common genes
data1_common = data1.loc[common_genes]
data2_common = data2.loc[common_genes]

# Compute variance for each gene in the common data
variances = data1_common.var(axis=1) + data2_common.var(axis=1)

# Select the top 500 most variable genes
top_genes = variances.nlargest(500).index
data1_top = data1_common.loc[top_genes]
data2_top = data2_common.loc[top_genes]

# Standardize the gene counts
scaler = StandardScaler()
data1_standardized = pd.DataFrame(
    scaler.fit_transform(data1_top.T).T, index=data1_top.index, columns=data1_top.columns
)
data2_standardized = pd.DataFrame(
    scaler.fit_transform(data2_top.T).T, index=data2_top.index, columns=data2_top.columns
)

# Perform PCA to reduce to 3 components
pca1 = PCA(n_components=3)
pca2 = PCA(n_components=3)
pca1_results = pd.DataFrame(
    pca1.fit_transform(data1_standardized.T), 
    index=data1_standardized.columns, 
    columns=["PC1", "PC2", "PC3"]
)
pca2_results = pd.DataFrame(
    pca2.fit_transform(data2_standardized.T), 
    index=data2_standardized.columns, 
    columns=["PC1", "PC2", "PC3"]
)

# Merge PCA results with metadata
metadata_subset = metadata.set_index("sample")
pca1_metadata = pca1_results.join(metadata_subset, how="inner")
pca2_metadata = pca2_results.join(metadata_subset, how="inner")

# Combine PCA results for both datasets
pca1_metadata['Dataset'] = 'Data1'
pca2_metadata['Dataset'] = 'Data2'
pca_combined = pd.concat([pca1_metadata, pca2_metadata], axis=0)

# Marker mapping
pca_combined['Marker'] = pca_combined['thermal.tolerance'].map(
    {'resistant': 'o', 'susceptible': 'x'}
)
# pca_combined['Marker Size'] = pca_combined['day'].map({0: 5, 30: 15})
pca_combined['Marker Size'] = pca_combined['day'].map({0: 3, 30: 8})


# Marker size adjustment is necessary for the 3D plot
def marker_size(row):
    if row['day'] == 30:
        return row['Marker Size'] * (1.5 if row['Marker'] == 'o' else 0.6)
    return row['Marker Size'] * (2.2 if row['Marker'] == 'o' else 1.4)


pca_combined['Scaled Marker Size'] = pca_combined.apply(
    lambda row: marker_size(row), axis=1
)

# Interactive 3D PCA plot with updated marker styles
fig = go.Figure()

# Add Data1 points with markers
data1_points = pca_combined[pca_combined['Dataset'] == 'Data1']
fig.add_trace(go.Scatter3d(
    x=data1_points['PC1'],
    y=data1_points['PC2'],
    z=data1_points['PC3'],
    mode='markers+text',
    marker=dict(
        size=data1_points['Scaled Marker Size'],
        color='blue',
        symbol=data1_points['Marker'].map({'o': 'circle', 'x': 'x'})
    ),
    name='Data1',
    text=data1_points.index.str[-2:],  # Sample ID as last two digits
    textposition='top center'
))

# Add Data2 points with markers
data2_points = pca_combined[pca_combined['Dataset'] == 'Data2']
fig.add_trace(go.Scatter3d(
    x=data2_points['PC1'],
    y=data2_points['PC2'],
    z=data2_points['PC3'],
    mode='markers+text',
    marker=dict(
        size=data2_points['Scaled Marker Size'],
        color='red',
        symbol=data2_points['Marker'].map({'o': 'circle', 'x': 'x'})
    ),
    name='Data2',
    text=data2_points.index.str[-2:],  # Sample ID as last two digits
    textposition='top center'
))

# Update layout
fig.update_layout(
    title="Interactive 3D PCA Plot",
    scene=dict(
        xaxis_title=f'PC1 ({pca1.explained_variance_ratio_[0]*100:.2f}%)',
        yaxis_title=f'PC2 ({pca1.explained_variance_ratio_[1]*100:.2f}%)',
        zaxis_title=f'PC3 ({pca1.explained_variance_ratio_[2]*100:.2f}%)',
    ),
    scene_camera=dict(
        up=dict(x=0, y=1, z=0),
        eye=dict(x=0, y=0, z=2.5)  # Adjust perspective to make PC1 horizontal and PC2 vertical
    ),
    legend=dict(title="Dataset")
)

# Save interactive 3D plot as HTML
html_path = f'{destdir}/Interactive_3D_PCA_Plot.html'
fig.write_html(html_path)


# Create static 2D PCA plots with markers
def plot_2d_pca(data, x, y, title, xlabel, ylabel, filename):
    fig, ax = plt.subplots(figsize=(8, 6))
    for (thermal_tolerance, marker), color in [
        (('resistant', 'o'), 'blue'), (('susceptible', 'x'), 'red')
    ]:
        subset = data[(data['thermal.tolerance'] == thermal_tolerance) & (data['Dataset'] == 'Data1')]
        ax.scatter(
            subset[x], subset[y], c='blue', alpha=0.7, label=f'{thermal_tolerance} Data1',
            marker=marker, s=subset['Marker Size'] * 10
        )
        for i, txt in enumerate(subset.index):
            ax.annotate(txt[-2:], (subset[x].iloc[i], subset[y].iloc[i]), fontsize=8)

        subset = data[(data['thermal.tolerance'] == thermal_tolerance) & (data['Dataset'] == 'Data2')]
        ax.scatter(
            subset[x], subset[y], c='red', alpha=0.7, label=f'{thermal_tolerance} Data2',
            marker=marker, s=subset['Marker Size'] * 10
        )
        for i, txt in enumerate(subset.index):
            ax.annotate(txt[-2:], (subset[x].iloc[i], subset[y].iloc[i]), fontsize=8)

    ax.set_xlabel(xlabel)
    ax.set_ylabel(ylabel)
    plt.title(title)
    handles, _labels = ax.get_legend_handles_labels()
    handles.extend([
        plt.Line2D([0], [0], marker='o', color='w', markerfacecolor='gray', markersize=6, label='Day 0'),
        plt.Line2D([0], [0], marker='o', color='w', markerfacecolor='gray', markersize=10, label='Day 30')
    ])
    ax.legend(handles=handles, title="Legend")
    plt.savefig(f'{destdir}{filename}.png')
    plt.close()

# Regenerate 2D PCA projections
plot_2d_pca(
    pca_combined, 'PC1', 'PC2',
    "2D PCA Projection (PC1 vs PC2)", 
    f'PC1 ({pca1.explained_variance_ratio_[0]*100:.2f}%)',
    f'PC2 ({pca1.explained_variance_ratio_[1]*100:.2f}%)',
    "2D_PCA_PC1_PC2"
)
plot_2d_pca(
    pca_combined, 'PC1', 'PC3',
    "2D PCA Projection (PC1 vs PC3)", 
    f'PC1 ({pca1.explained_variance_ratio_[0]*100:.2f}%)',
    f'PC3 ({pca1.explained_variance_ratio_[2]*100:.2f}%)',
    "2D_PCA_PC1_PC3"
)