snputils.visualization.admixture_manhattan_plot

 1import os
 2import numpy as np
 3import pandas as pd
 4from typing import Optional, Tuple, Dict
 5import matplotlib.pyplot as plt
 6
 7def manhattan_plot(
 8    input_file: str, 
 9    colors: list,
10    significance_threshold: float = 0.05,
11    figsize: Optional[Tuple[float, float]] = None,
12    title: Optional[str] = None,
13    fontsize: Optional[Dict[str, float]] = None,
14    save: Optional[bool] = None,
15    output_filename: Optional[str] = None,
16):
17    """
18    Generates a Manhattan plot from an input file the results of an admixture mapping association. 
19    The plot is highly customizable, allowing users to specify colors for each chromosome and apply bonferroni correction for p-values.
20
21    Args:
22        input_file: Path to the input file containing columns '#CHROM', 'POS', and 'P'. default plik2 output files are supported.
23        colors: List of colors to apply to each chromosome. the chromosome number is used as an index to select the color.
24        significance_threshold: Significance threshold for p-values. Default is 0.05.
25        figsize: Optional tuple to specify figure dimensions (width, height).
26        title: Plot title. If None, no title is shown.
27        fontsize: Dictionary with font sizes for title, labels, and legend.
28        save: If True, saves the plot to a file. If None, the plot is not saved.
29        output_filename: Filename for saving the plot (PNG).
30    """
31    # Read the input file
32    df = pd.read_csv(input_file, sep='\t')
33
34    # Calculate the maximum distance within each chromosome to scale absolute positions
35    max_distance = 0
36    for chrom, chrom_data in df.groupby('#CHROM'):
37        chrom_max_pos = chrom_data['POS'].max()
38        if chrom_max_pos > max_distance:
39            max_distance = chrom_max_pos
40
41    # Calculate absolute positions for each SNP
42    df['ABS_POS'] = df['POS'] + max_distance * df['#CHROM']
43
44    # Bonferroni threshold
45    bonferroni_threshold = significance_threshold / len(df)
46
47    # Create the plot
48    plt.figure(figsize=figsize)
49    chrom_offsets = {chrom: max_distance * (chrom - 1) for chrom in range(1, 23)}
50
51    # Display Manhattan plot points for each chromosome
52    for i, (chrom, chrom_data) in enumerate(df.groupby('#CHROM')):
53        chrom_data['ABS_POS'] = chrom_data['POS'] + chrom_offsets[chrom]
54        plt.scatter(chrom_data['ABS_POS'], -np.log10(chrom_data['P']), 
55                    color=colors[int(chrom+1) % len(colors)])
56
57    # X-axis settings
58    plt.xlim(0, 22 * max_distance)
59    chrom_labels = [str(c) for c in range(1, 23)]
60    chrom_positions = [chrom_offsets[c] + max_distance / 2 for c in range(1, 23)]
61    plt.xticks(chrom_positions, chrom_labels)
62
63    # Significance thresholds
64    plt.axhline(y=-np.log10(bonferroni_threshold), color='r', linestyle='--', label='Bonferroni')
65
66    # Labels and title
67    if title:
68        plt.title(title, fontsize=fontsize.get('title', 20))
69    plt.xlabel('Chromosomes', fontsize=fontsize.get('xlabel', 15))
70    plt.ylabel('-log10(p-value)', fontsize=fontsize.get('ylabel', 15))
71    plt.legend(fontsize=fontsize.get('legend', 15))
72
73    # Save the plot
74    plt.tight_layout()
75    if save:
76        plt.savefig(output_filename)
77    plt.show()