Positive selection on brain cis-regulatory elements in the human lineage drives changes in gene expression and susceptibility to neuropsychiatric disorders.

1. import pandas as pd
2. import numpy as np
3. from scipy import stats
4. import warnings
5. warnings.filterwarnings('ignore')
7. # Install sklearn if not available
8. try:
9. from sklearn.linear_model import LogisticRegression
10. except ImportError:
11. import subprocess
12. import sys
13. print("Installing scikit-learn...")
14. subprocess.check_call([sys.executable, "-m", "pip", "install", "scikit-learn"])
15. from sklearn.linear_model import LogisticRegression
17. def benjamini_hochberg_fdr(p_values, alpha=0.05):
18. """Apply Benjamini-Hochberg FDR correction"""
19. p_values = np.array(p_values)
20. n = len(p_values)
22. # Sort p-values and keep track of original indices
23. sorted_indices = np.argsort(p_values)
24. sorted_p = p_values[sorted_indices]
26. # Calculate adjusted p-values
27. adjusted_p = np.zeros(n)
28. for i in range(n-1, -1, -1):
29. if i == n-1:
30. adjusted_p[sorted_indices[i]] = sorted_p[i]
31. else:
32. adjusted_p[sorted_indices[i]] = min(
33. sorted_p[i] * n / (i + 1),
34. adjusted_p[sorted_indices[i+1]]
35. )
37. return np.clip(adjusted_p, 0, 1)
39. def run_wgcna_enrichment(df, gene_list, gene_type_name):
40. """Run WGCNA enrichment analysis for a specific gene list"""
42. print(f"\n{'='*70}")
43. print(f"ANALYZING {gene_type_name.upper()} GENES")
44. print(f"{'='*70}")
46. # Create HAR status for all genes in WGCNA background
47. df['HAR_status'] = df['ID'].isin(gene_list).astype(int)
48. genes_in_wgcna = df['HAR_status'].sum()
49. print(f"Total {gene_type_name} genes in list: {len(gene_list)}")
50. print(f"{gene_type_name} genes found in WGCNA background: {genes_in_wgcna}")
52. if genes_in_wgcna == 0:
53. print(f"WARNING: No {gene_type_name} genes found in WGCNA background!")
54. return None
56. # Test enrichment for each module using logistic linear regression
57. print(f"\nTesting enrichment using logistic linear regression across {df['Modules'].nunique()} modules...")
59. modules = sorted(df['Modules'].unique())
60. results = []
62. for module in modules:
63. # Create binary predictor: gene is in this module (1) or not (0)
64. in_module = (df['Modules'] == module).astype(int)
66. # Binary outcome: gene is HAR-associated (1) or not (0)
67. har_status = df['HAR_status']
69. # Skip modules with no genes or if no HAR genes exist
70. if in_module.sum() == 0 or har_status.sum() == 0:
71. continue
73. # Logistic linear regression as described in the paper
74. X = in_module.values.reshape(-1, 1)
75. y = har_status.values
77. try:
78. # Fit logistic regression
79. lr = LogisticRegression(fit_intercept=True, max_iter=1000)
80. lr.fit(X, y)
82. # Get coefficient (log odds ratio)
83. coef = lr.coef_[0][0]
85. # Calculate p-value using likelihood ratio test
86. lr_null = LogisticRegression(fit_intercept=True, max_iter=1000)
87. lr_null.fit(np.zeros((len(y), 1)), y)
89. # Calculate log-likelihoods
90. y_pred_full = lr.predict_proba(X)[:, 1]
91. y_pred_null = lr_null.predict_proba(np.zeros((len(y), 1)))[:, 1]
93. # Avoid numerical issues
94. epsilon = 1e-15
95. y_pred_full = np.clip(y_pred_full, epsilon, 1-epsilon)
96. y_pred_null = np.clip(y_pred_null, epsilon, 1-epsilon)
98. ll_full = np.sum(y * np.log(y_pred_full) + (1-y) * np.log(1-y_pred_full))
99. ll_null = np.sum(y * np.log(y_pred_null) + (1-y) * np.log(1-y_pred_null))
101. # Likelihood ratio test statistic
102. lr_stat = 2 * (ll_full - ll_null)
103. p_value = 1 - stats.chi2.cdf(lr_stat, df=1)
105. except Exception as e:
106. print(f"Logistic regression failed for module {module}: {e}")
107. continue
109. # Calculate descriptive statistics
110. har_in_module = ((har_status == 1) & (in_module == 1)).sum()
111. total_in_module = in_module.sum()
112. total_har = har_status.sum()
114. # Expected number of HAR genes in module under null hypothesis
115. expected_har_in_module = (total_har * total_in_module) / len(df)
116. fold_enrichment = har_in_module / max(expected_har_in_module, 1)
118. # Odds ratio
119. odds_ratio = np.exp(coef)
121. results.append({
122. 'Module': module,
123. 'Total_genes_in_module': total_in_module,
124. 'HAR_genes_in_module': har_in_module,
125. 'Total_HAR_genes': total_har,
126. 'Expected_HAR_in_module': expected_har_in_module,
127. 'Fold_enrichment': fold_enrichment,
128. 'Log_odds_ratio': coef,
129. 'Odds_ratio': odds_ratio,
130. 'P_value': p_value
131. })
133. # Convert results to DataFrame
134. results_df = pd.DataFrame(results)
136. # Multiple testing correction
137. results_df['P_value_bonferroni'] = results_df['P_value'] * len(results_df)
138. results_df['P_value_bonferroni'] = results_df['P_value_bonferroni'].clip(upper=1.0)
140. # FDR correction
141. results_df['P_value_fdr'] = benjamini_hochberg_fdr(results_df['P_value'], alpha=0.05)
143. # Sort by p-value
144. results_df = results_df.sort_values('P_value')
146. # Display summary
147. print(f"\nLogistic linear regression analysis complete!")
148. print(f"Modules tested: {len(results_df)}")
149. print(f"Significant modules (p < 0.05): {(results_df['P_value'] < 0.05).sum()}")
150. print(f"Significant modules (Bonferroni p < 0.05): {(results_df['P_value_bonferroni'] < 0.05).sum()}")
151. print(f"Significant modules (FDR p < 0.05): {(results_df['P_value_fdr'] < 0.05).sum()}")
153. # Show top enriched modules
154. print(f"\nTop 10 most significantly enriched modules:")
155. display_cols = ['Module', 'HAR_genes_in_module', 'Total_genes_in_module',
156. 'Fold_enrichment', 'Odds_ratio', 'P_value', 'P_value_bonferroni', 'P_value_fdr']
157. print(results_df[display_cols].head(10))
159. return results_df
162. # Read the Excel file with WGCNA modules
163. print("Reading WGCNA modules data...")
164. df = pd.read_excel('Gene Modules.xlsx', header=2)
166. # Clean the data
167. df = df.dropna(subset=['ID'])
168. df['ID'] = df['ID'].astype(str).str.strip()
170. # Handle merged cells - forward fill module numbers
171. df['Modules'] = df['Modules'].fillna(method='ffill')
172. df['Modules'] = df['Modules'].astype(int)
174. print(f"Total genes in WGCNA: {len(df)}")
175. print(f"Total modules: {df['Modules'].nunique()}")
177. # Define gene list files and their names
178. gene_lists = {
179. 'shared': 'shared_HPS-bCREs.txt',
180. 'fetal': 'fetal_HPS-bCREs.txt',
181. 'adult': 'adult_all_HPS-bCREs.txt'
182. }
184. # Process each gene list
185. for gene_type, filename in gene_lists.items():
186. try:
187. # Read gene list
188. print(f"\n\nReading {gene_type} genes list from {filename}...")
189. with open(filename, 'r') as f:
190. gene_list = [line.strip() for line in f if line.strip()]
192. # Run enrichment analysis
193. results_df = run_wgcna_enrichment(df.copy(), gene_list, gene_type)
195. if results_df is not None:
196. # Save results
197. output_filename = f'HAR_WGCNA_enrichment_results_{gene_type}.csv'
198. results_df.to_csv(output_filename, index=False)
199. print(f"\nResults saved to '{output_filename}'")
201. except FileNotFoundError:
202. print(f"\nWARNING: Could not find file '{filename}' - skipping {gene_type} analysis")
203. except Exception as e:
204. print(f"\nERROR processing {gene_type} genes: {e}")
206. print(f"\n\n{'='*70}")
207. print("ALL ANALYSES COMPLETE")
208. print(f"{'='*70}")
209. print("\nMethod used (as per paper):")
210. print("- Logistic linear regression")
211. print(f"- Background: {len(df)} genes included in WGCNA analysis")
212. print("- No regression of exome or gene length (promoter-based interactions)")
213. print("- Multiple testing: Bonferroni and FDR corrections applied")

WGCNA.py at commit 5a3e90a, no license · at the source