prs_comparisons/plot_feature.py at main · jmillerlab/prs_comparisons · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import sys
# Load your CSV file
####df = pd.read_csv("percentiles_race_round.csv")
#df = pd.read_csv("pop_projections/percentiles/EUR")
df = pd.read_csv(sys.argv[1])


'''
race='Black'
#race='White'
#race='Asian'
# Filter for Race == 'White'
race_df = df[df['Race'] == f'{race}']
'''
#### Take out this line and remove comment above to work on a specific reported race
race_df = df
####


ad_count = (race_df['Class'] == 'AD').sum()
baseline_percentage = (ad_count / len(race_df)) *100
ax=plt.figure(figsize=(12, 6))
plt.axhline(y=baseline_percentage, color='blue', linestyle='--',alpha=0.2)
#plt.axhline(y=baseline_percentage, color='black', linestyle='--', label='Baseline')
#plt.text(100, baseline_percentage, 'Baseline', color='black', va='bottom', ha='right')
colorscheme = {'e33' : 'purple', 'e34' : 'green', 'e44' : 'red' }
titlescheme = {'e33' : 'APOE ε3/ε3', 'e34' : 'APOE ε3/ε4', 'e44' : 'APOE ε4/ε4' }
#print(len(race_df))
#sys.exit()
for apoe in ['e33','e34','e44']:
    subset = race_df[race_df['APOE_reported'] == apoe]
    if len(subset) > 0:
        ad_count = (subset['Class'] == 'AD').sum()
        percentage = 100*(ad_count / len(subset))
        color = colorscheme[apoe]
        title = titlescheme[apoe]
        plt.axhline(y=percentage, color=color, linestyle='--', alpha=0.2)
        #plt.axhline(y=percentage, color=color, linestyle='--', label=title)
        #plt.text(97, percentage, title, color='black', va='bottom', ha='left')
        print(apoe,percentage,len(subset), ad_count, len(subset)-ad_count)
    else:
        percentages.append(np.nan)


#sys.exit()

# Sort unique values of the Kunkle column
#column='Bellenguez_s0.2_lambda0.0233572146909012'
#column='Jansen_1e.04_0.07_sparse'
column='Kunkle_0.1_0.07_sparse'
thresholds = sorted(race_df[column].unique())

# Calculate percentage of 'AD' for each threshold
percentages = []
apoe_percentages={'e33':[],'e34':[],'e44':[]} #e33,e34,e44
for t in thresholds:
    subset = race_df[race_df[column] >= t]
    if len(subset) > 0:
        ad_count = (subset['Class'] == 'AD').sum()
        percentage = ad_count / len(subset)
        percentages.append(percentage*100)
        for genotype in apoe_percentages:
            subset2= subset[subset['APOE_reported'] == genotype]
            percent = 100* ((subset2['Class'] == 'AD').sum() / len(subset2))
            apoe_percentages[genotype].append(percent)


        if thresholds[-1] == t:
            subset2= subset[subset['Class'] == 'AD']
            print("AD Cases",len(subset2), (subset2['Sex'] == 'M').sum(), (subset2['Sex'] == 'F').sum())
            subset2= subset[subset['Class'] == 'Control']
            print("Controls",len(subset2), (subset2['Sex'] == 'M').sum(), (subset2['Sex'] == 'F').sum())
            subset2= subset[subset['APOE_reported'] == 'e24']
            print("e24",len(subset2), (subset2['Class'] == 'AD').sum(), (subset2['Class'] == 'Control').sum())
            subset2= subset[subset['APOE_reported'] == 'e23']
            print("e23",len(subset2), (subset2['Class'] == 'AD').sum(), (subset2['Class'] == 'Control').sum())
            subset2= subset[subset['APOE_reported'] == 'e33']
            print("e33",len(subset2), (subset2['Class'] == 'AD').sum(), (subset2['Class'] == 'Control').sum())
            subset2= subset[subset['APOE_reported'] == 'e34']
            print("e34",len(subset2), (subset2['Class'] == 'AD').sum(), (subset2['Class'] == 'Control').sum())
            subset2= subset[subset['APOE_reported'] == 'e44']
            print("e44",len(subset2), (subset2['Class'] == 'AD').sum(), (subset2['Class'] == 'Control').sum())
    else:
        percentages.append(np.nan)

# Plotting
#plt.figure(figsize=(10, 6))
plt.plot(thresholds, percentages, label='Average', color='blue')
percent=round(percentages[-1])
#plt.text(0,percentages[0], f'PRS={percent}' + "%", color='blue',va="bottom",ha='left')
plt.text(0,percentages[0]+1, f'Average', color='blue',va="bottom",ha='left')

for apoe in ['e33','e34','e44']:
    plt.plot(thresholds, apoe_percentages[apoe], label=titlescheme[apoe], color=colorscheme[apoe])
    plt.text(0, apoe_percentages[apoe][0]+1, titlescheme[apoe], color=colorscheme[apoe], va='bottom', ha='left')

#plt.axhline(y=32.37, color='black', linestyle='--', label='Baseline')
#plt.text(thresholds[-1], 32.37, 'Baseline', color='black', va='bottom', ha='right')
'''
plt.axhline(y=26.93739425, color='orange', linestyle='--', label='APOE ε3/ε3')
plt.text(thresholds[0], 26.93739425, 'APOE ε3/ε3', color='black', va='bottom', ha='left')

plt.axhline(y=42.71660991, color='green', linestyle='--', label='APOE ε3/ε4')
plt.text(thresholds[0], 42.71660991, 'APOE ε3/ε4', color='black', va='bottom', ha='left')

plt.axhline(y=55.56872038, color='red', linestyle='--', label='APOE ε4/ε4')
plt.text(thresholds[0], 55.56872038, 'APOE ε4/ε4', color='black', va='bottom', ha='left')
'''
##plt.xlim(0,100.1)
plt.ylim(0,100.1)
plt.xlabel(f'{column} Precentile Risk')
plt.ylabel(f'Probability of AD Diagnosis')
#######plt.title(f'AD Risk vs PRS Percentile (Race: {race})')
plt.title(f'AD Risk vs PRS Percentile')
plt.legend(loc='upper left')
#plt.grid(True)
plt.tight_layout()
plt.show()