Cerebellar and frontal lobe alterations in pulmonary hypertension implicate an inflammatory lung vascular/brain axis.

Code ↔ Paper

2 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 2 matches

[1] § Results › Deficits in the L-CRBL–L-SFG axis are specific for PH. ↔ function_module.py, lines 448 to 502 · score 0.75 · central systolic blood, LV cardiac, blood pressure, ventricular rate · Report a mismatch
[2] § Results › A multicohort human PH study of brain structure and function performed across discovery and validation cohorts. ↔ function_module.py, lines 448 to 502 · score 0.53 · college degree, IQR, female, diabetes, median, BMI · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Python · 1,235 lines · 72 KB · no license · 2 matches

import sys,os
import pandas as pd
import numpy as np
import re
import seaborn as sns
import matplotlib.backends.backend_pdf
import matplotlib.pyplot as plt
import statsmodels.api as sm
import statsmodels.stats as ss
from statsmodels.formula.api import ols, logit, glm
from statsmodels.genmod import families
from statsmodels.stats.multicomp import pairwise_tukeyhsd
from statsmodels.stats.nonparametric import rank_compare_2indep
from statsmodels.stats.contingency_tables import mcnemar
from sklearn.impute import KNNImputer
from statsmodels.discrete.conditional_models import ConditionalLogit
import statsmodels.formula.api as smf
from statsmodels.tools.sm_exceptions import ConvergenceWarning
import pickle
from IPython.display import display
from statsmodels.miscmodels.ordinal_model import OrderedModel
import scipy.stats as st
from statsmodels.stats.weightstats import ztest
from statsmodels.stats import weightstats as sms
from statsmodels.graphics.regressionplots import plot_partregress_grid
from random import sample
from scipy.stats import shapiro
from scipy import stats
def filterData(df, kinship=True, onlyWhiteBritish=True):
if onlyWhiteBritish:
df = df[df['22006-0.0']==1] # only Caucasian ancestry (In_white_british_ancestry_subset)
df = df[df['31-0.0'].notna()] # sex
df = df[df['31-0.0']==df['22001-0.0']]
df = df[df['22019-0.0'].isna()]
if kinship:
df = df[df['22020-0.0']==1]
return(df)
def booleanizeEducation(df): # reference no college degree
df['6138-0.0'] = df['6138-0.0'].isin([1]).astype(int)
return(df)
def booleanizeEthnicity(df): # reference white
df['21000-0.0'] = (~df['21000-0.0'].isin([1001])).astype(int)
return(df)
def booleanizeSelfReportedDF20002(df, codes, field_name): # reference non diabetes
df[field_name] = df['20002-0.0'].isin(codes).astype(int)
return(df)
def runLogit_caseControl_formula(X2, formula, dv, returnFullResult=False, printResult=True, l1=True):
try:
res = logit(formula, X2).fit()
except np.linalg.LinAlgError as err:
if 'Singular matrix' in str(err):
return(pd.DataFrame([]))
else:
if printResult==True:
print(res.summary())
t = res.summary2().tables[1][['Coef.','Std.Err.','P>|z|','[0.025','0.975]']]
t = t.fillna(1)
t = t.drop('Intercept')
t = pd.concat([pd.DataFrame({'Coef.':[dv, res.nobs, np.sum(X2['Case']==1)], 'Std.Err.':["","",""], 'P>|z|':["","",""], '[0.025':["","",""], '0.975]':["","",""] },index=['Dep_Variable:','Sample_Count','Sample_Count_Case']),t])
return(t)
def runGLMLogit_severity_formula(X2, formula, dv, returnFullResult=False, printResult=True, l1=True, offset=None):
if isinstance(offset, type(None)):
try:
res = glm(formula, data=X2, family=families.Binomial()).fit()
except sm.tools.sm_exceptions.PerfectSeparationError as err:
if 'Perfect separation' in str(err):
return pd.DataFrame([]), pd.DataFrame([])
else:
t = res.summary2().tables[1][['Coef.','Std.Err.','P>|z|','[0.025','0.975]']]
else:
try:
res = glm(formula, offset=offset, data=X2, family=families.Binomial()).fit()
except sm.tools.sm_exceptions.PerfectSeparationError as err:
if 'Perfect separation' in str(err):
return pd.DataFrame([]), pd.DataFrame([])
else:
t = res.summary2().tables[1][['Coef.','Std.Err.','P>|z|','[0.025','0.975]']]
if printResult==True:
print(res.summary())
t = t.fillna(1)
if t.index.isin(['Intercept']).sum() >0:
t = t.drop('Intercept')
if 'TimeOffset' in X2.columns:
Coef_sx = t['Coef.']*X2.drop(['TimeOffset','IDF'],axis=1).std(axis=0)
else:
Coef_sx = t['Coef.']*X2.drop(['IDF'],axis=1).std(axis=0)
coefs = pd.concat([Coef_sx, Coef_sx / X2['IDF'].std()],axis=1)
coefs.rename({0:'Coef.sx',1:'Coef_sxsy'},axis=1,inplace=True)
t = pd.concat([t, coefs],axis=1)
t = pd.concat([pd.DataFrame({'Coef.':[dv, X2.shape[0], res.nobs], 'Std.Err.':["","",""], 'P>|z|':["","",""], '[0.025':["","",""], '0.975]':["","",""], 'Coef.sx':["","",""], 'Coef_sxsy':["","",""]},index=['Dep_Variable:','Sample_Count','Sample_Count_Case']),t])
return(res, t)
def bootstrapGLMLogitregression(X2, formula, dvmax, resT, pcol, ntrials=10000, timeOffset=False):
if timeOffset:
params = pd.DataFrame([],index=X2.columns.drop(['IDF','TimeOffset']),columns=list(range(0,ntrials)))
else:
params = pd.DataFrame([],index=X2.columns.drop(['IDF']),columns=list(range(0,ntrials)))
for trial in range(ntrials):
data_sample = X2.sample(frac=1, replace=True)
if timeOffset:
fitmodel,res = runGLMLogit_severity_formula(data_sample, formula, dvmax, offset=data_sample['TimeOffset'], printResult=False, l1=False)
else:
fitmodel,res = runGLMLogit_severity_formula(data_sample, formula, dvmax, printResult=False, l1=False)
params[trial] = res['Coef.'][3:]
point_estimate = np.mean(params,axis=1)
bootstrap_standard_error = np.std(params,ddof=1,axis=1)
z_critical = stats.norm.ppf(1 - (1 - 0.95) / 2) # confidence_level = 0.95
resT['[0.025'][3:] = point_estimate - z_critical * bootstrap_standard_error
resT['0.975]'][3:] = point_estimate + z_critical * bootstrap_standard_error
shapiroP = stats.shapiro(params.loc[caseCol,:])[1]
print('Shapiro test:'); print(shapiroP)
if pvalMethod=='ztest':
pval = ztest(params.loc[caseCol,:], value=0)
resT.loc[caseCol,pcol] = pval[1]
elif pvalMethod=='count-based':
gt = (params.loc[caseCol,:] >0).sum()
lt = (params.loc[caseCol,:] <0).sum()
resT.loc[caseCol,pcol] = np.min([gt,lt]) / (gt+lt)
return(resT)
def permutationGLMLogitregression(X2, formula, dvmax, resT, caseCol, pcol, pvalMethod, ntrials, timeOffset=False):
if timeOffset:
params = pd.DataFrame([],index=X2.columns.drop(['IDF','TimeOffset']),columns=list(range(0,ntrials)))
else:
params = pd.DataFrame([],index=X2.columns.drop(['IDF']),columns=list(range(0,ntrials)))
for trial in range(ntrials):
data_sample = X2.copy()
pred = data_sample['Predictor'].sample(frac=1)
pred.index = data_sample.index
data_sample['Predictor'] = pred
if timeOffset:
fitmodel,res = runGLMLogit_severity_formula(data_sample, formula, dvmax, offset=data_sample['TimeOffset'], printResult=False, l1=False)
else:
fitmodel,res = runGLMLogit_severity_formula(data_sample, formula, dvmax, printResult=False, l1=False)
params[trial] = res['Coef.'][3:]
params = params[params.isna().sum(axis=1)!=ntrials]
sorted_perm_stats = np.sort(np.array(params.loc[caseCol,:]))
sorted_perm_stats = sorted_perm_stats - resT.loc[caseCol,'Coef.']
resT.loc[caseCol,'[0.025'] = np.percentile(sorted_perm_stats, 2.5)
resT.loc[caseCol,'0.975]'] = np.percentile(sorted_perm_stats, 97.5)
shapiroP = stats.shapiro(params.loc[caseCol,:])[1]
print('Shapiro test:'); print(shapiroP)
if pvalMethod=='ztest':
pval = ztest(params.loc[caseCol,:], value=resT.loc[caseCol,'Coef.'])
resT.loc[caseCol,pcol] = pval[1]
elif pvalMethod=='count-based':
gt = (params.loc[caseCol,:] >resT.loc[caseCol,'Coef.']).sum()
lt = (params.loc[caseCol,:] <resT.loc[caseCol,'Coef.']).sum()
resT.loc[caseCol,pcol] = np.min([gt,lt]) / (gt+lt)
return(resT)
def bootstrapLoRegression(X2, formula, dvmax, resT, caseCol, pcol, pvalMethod, ntrials):
params = pd.DataFrame([],index=X2.columns[1:],columns=list(range(0,ntrials)))
for trial in range(ntrials):
data_sample = X2.sample(frac=1, replace=True)
res = runLogit_caseControl_formula(data_sample, formula, dvmax, printResult=False, l1=False)
if res.empty: # skip 0
continue
params[trial] = res['Coef.'][3:]
params = params[params.isna().sum(axis=1)!=ntrials]
params = params.loc[:,~params.loc['Case'].isna()] # remove NAs
point_estimate = np.mean(params,axis=1)
bootstrap_standard_error = np.std(params,ddof=1,axis=1)
z_critical = stats.norm.ppf(1 - (1 - 0.95) / 2) # confidence_level = 0.95
resT['[0.025'][3:] = point_estimate - z_critical * bootstrap_standard_error
resT['0.975]'][3:] = point_estimate + z_critical * bootstrap_standard_error
shapiroP = stats.shapiro(params.loc[caseCol,:])[1]
print('Shapiro test:'); print(shapiroP)
if pvalMethod=='ztest':
pval = ztest(params.loc[caseCol,:], value=0)
resT.loc[caseCol,pcol] = pval[1]
elif pvalMethod=='count-based':
gt = (params.loc[caseCol,:] >0).sum()
lt = (params.loc[caseCol,:] <0).sum()
resT.loc[caseCol,pcol] = np.min([gt,lt]) / (gt+lt)
return(resT)
def permutationLoRegression_ntrials(X2, formula, dvmax, resT, caseCol, pcol, pvalMethod, ntrials):
params = pd.DataFrame([],index=X2.columns[1:],columns=list(range(0,ntrials)))
for trial in range(ntrials):
data_sample = X2.copy()
data_sample[caseCol] = 0
data_sample[caseCol].iloc[sample(range(X2.shape[0]), k=X2[caseCol].sum())] = 1
res = runLogit_caseControl_formula(data_sample, formula, dvmax, printResult=False)
params[trial] = res['Coef.'][3:]
params = params[params.isna().sum(axis=1)!=ntrials]
sorted_perm_stats = np.sort(np.array(params.loc[caseCol,:]))
sorted_perm_stats = sorted_perm_stats - resT.loc[caseCol,'Coef.']
resT.loc[caseCol,'[0.025'] = np.percentile(sorted_perm_stats, 2.5)
resT.loc[caseCol,'0.975]'] = np.percentile(sorted_perm_stats, 97.5)
shapiroP = stats.shapiro(params.loc[caseCol,:])[1]
print('Shapiro test:'); print(shapiroP)
if pvalMethod=='ztest':
pval = ztest(params.loc[caseCol,:], value=resT.loc[caseCol,'Coef.'])
resT.loc[caseCol,pcol] = pval[1]
elif pvalMethod=='count-based':
gt = (params.loc[caseCol,:] >resT.loc[caseCol,'Coef.']).sum()
lt = (params.loc[caseCol,:] <resT.loc[caseCol,'Coef.']).sum()
resT.loc[caseCol,pcol] = np.min([gt,lt]) / (gt+lt)
return(resT)
def permutationGroupMeanLoRegression(X2, formula, dvmax, resT, caseCol, pcol, pvalMethod, ntrials):
resT.iloc[3:,:] = np.nan
try:
res = logit(formula, X2).fit()
except np.linalg.LinAlgError as err:
if 'Singular matrix' in str(err):
return(pd.DataFrame([]))
else:
moT = logit(formula, X2).fit()
temp = moT.predict(X2)
temp[temp >=0.5] = 1
temp[temp < 0.5] = 0
print('Residual normality Shapiro test:'); print(stats.shapiro(moT.resid)[1])
resT.loc[caseCol,'Coef.'] = temp[X2[caseCol]==1].mean() - temp[X2[caseCol]==0].mean()
resT.loc[caseCol,'Std.Err.'] = pooledSD(temp[X2[caseCol]==1], temp[X2[caseCol]==0]) # pooledSD
print('Cohens d:');
print(resT.loc[caseCol,'Coef.'] / resT.loc[caseCol,'Std.Err.']) # cohen's d
params = pd.DataFrame([],index=X2.columns[1:],columns=list(range(0,ntrials)))
for trial in range(ntrials):
data_sample = X2.copy()
data_sample[caseCol] = 0
data_sample[caseCol].iloc[sample(range(X2.shape[0]), k=X2[caseCol].sum())] = 1
try:
res = logit(formula, data_sample).fit()
except np.linalg.LinAlgError as err:
if 'Singular matrix' in str(err):
continue
else:
temp = logit(formula, data_sample).fit().predict(data_sample)
if temp.empty: # skip 0
continue
params.loc[caseCol,trial] = temp[data_sample[caseCol]==1].mean() - temp[data_sample[caseCol]==0].mean()
params = params.loc[:,~params.loc['Case'].isna()] # remove NAs
z_critical = stats.norm.ppf(1 - (1 - 0.95) / 2) # confidence_level = 0.95
resT.loc[caseCol,'[0.025'] = resT.loc[caseCol,'Coef.'] - z_critical * resT.loc[caseCol,'Std.Err.']
resT.loc[caseCol,'0.975]'] = resT.loc[caseCol,'Coef.'] + z_critical * resT.loc[caseCol,'Std.Err.']
shapiroP = stats.shapiro(params.loc[caseCol,:])[1]
print('Shapiro test:'); print(shapiroP)
if pvalMethod=='ztest':
pval = ztest(params.loc[caseCol,:], value=resT.loc[caseCol,'Coef.'])
resT.loc[caseCol,pcol] = pval[1]
elif pvalMethod=='count-based':
gt = (params.loc[caseCol,:] >resT.loc[caseCol,'Coef.']).sum()
lt = (params.loc[caseCol,:] <resT.loc[caseCol,'Coef.']).sum()
resT.loc[caseCol,pcol] = np.min([gt,lt]) / (gt+lt)
return(resT)
def runGLMOLS_severity_formula(X2, formula, dv, returnFullResult=False, printResult=True, l1=True, offset=None):
if isinstance(offset, type(None)):
try:
res = glm(formula, data=X2).fit()
except sm.tools.sm_exceptions.PerfectSeparationError as err:
if 'Perfect separation' in str(err):
return pd.DataFrame([]), pd.DataFrame([])
else:
t = res.summary2().tables[1][['Coef.','Std.Err.','P>|z|','[0.025','0.975]']]
else:
try:
res = glm(formula, offset=offset, data=X2).fit()
except sm.tools.sm_exceptions.PerfectSeparationError as err:
if 'Perfect separation' in str(err):
return pd.DataFrame([]), pd.DataFrame([])
else:
t = res.summary2().tables[1][['Coef.','Std.Err.','P>|z|','[0.025','0.975]']]
if printResult==True:
print(res.summary())
t = t.fillna(1)
if t.index.isin(['Intercept']).sum() >0:
t = t.drop('Intercept')
if 'TimeOffset' in X2.columns:
Coef_sx = t['Coef.']*X2.drop(['TimeOffset','IDF'],axis=1).std(axis=0)
else:
Coef_sx = t['Coef.']*X2.drop(['IDF'],axis=1).std(axis=0)
coefs = pd.concat([Coef_sx, Coef_sx / X2['IDF'].std()],axis=1)
coefs.rename({0:'Coef.sx',1:'Coef_sxsy'},axis=1,inplace=True)
t = pd.concat([t, coefs],axis=1)
t = pd.concat([pd.DataFrame({'Coef.':[dv, X2.shape[0], res.nobs], 'Std.Err.':["","",""], 'P>|z|':["","",""], '[0.025':["","",""], '0.975]':["","",""], 'Coef.sx':["","",""], 'Coef_sxsy':["","",""]},index=['Dep_Variable:','Sample_Count','Sample_Count_Case']),t])
return(res, t)
def bootstrapGLMOLSregression(X2, formula, dvmax, resT, caseCol, pcol, pvalMethod, ntrials, timeOffset=False):
if timeOffset:
params = pd.DataFrame([],index=X2.columns.drop(['IDF','TimeOffset']),columns=list(range(0,ntrials)))
else:
params = pd.DataFrame([],index=X2.columns.drop(['IDF']),columns=list(range(0,ntrials)))
for trial in range(ntrials):
data_sample = X2.sample(frac=1, replace=True)
if timeOffset:
fitmodel = runGLMOLS_severity_formula(data_sample, formula, dvmax, offset=data_sample['TimeOffset'], printResult=False, l1=False)
if isinstance(fitmodel, tuple):
fitmodel,res = runGLMOLS_severity_formula(data_sample, formula, dvmax, offset=data_sample['TimeOffset'], printResult=False, l1=False)
params[trial] = res['Coef.'][3:]
else:
fitmodel = runGLMOLS_severity_formula(data_sample, formula, dvmax, printResult=False, l1=False)
if isinstance(fitmodel, tuple):
fitmodel,res = runGLMOLS_severity_formula(data_sample, formula, dvmax, printResult=False, l1=False)
params[trial] = res['Coef.'][3:]
params = params.loc[:,params.isna().sum(axis=0)!=params.shape[0]]
if params.shape[1] < 10:
return(resT)
point_estimate = np.mean(params,axis=1)
bootstrap_standard_error = np.std(params,ddof=1,axis=1)
z_critical = stats.norm.ppf(1 - (1 - 0.95) / 2) # confidence_level = 0.95
resT['[0.025'][3:] = point_estimate - z_critical * bootstrap_standard_error
resT['0.975]'][3:] = point_estimate + z_critical * bootstrap_standard_error
shapiroP = stats.shapiro(params.loc[caseCol,:])[1]
print('Shapiro test:'); print(shapiroP)
if pvalMethod=='ztest':
pval = ztest(params.loc[caseCol,:], value=0)
resT.loc[caseCol,pcol] = pval[1]
elif pvalMethod=='count-based':
gt = (params.loc[caseCol,:] >0).sum()
lt = (params.loc[caseCol,:] <0).sum()
resT.loc[caseCol,pcol] = np.min([gt,lt]) / (gt+lt)
return(resT)
def permutationGLMOLSregression(X2, formula, dvmax, resT, caseCol, pcol, pvalMethod, ntrials, timeOffset=False):
if timeOffset: # timeoffset
params = pd.DataFrame([],index=X2.columns.drop(['IDF','TimeOffset']),columns=list(range(0,ntrials)))
else:
params = pd.DataFrame([],index=X2.columns.drop(['IDF']),columns=list(range(0,ntrials)))
for trial in range(ntrials):
data_sample = X2.copy()
pred = data_sample['Predictor'].sample(frac=1)
pred.index = data_sample.index
data_sample['Predictor'] = pred
if timeOffset:
fitmodel,res = runGLMOLS_severity_formula(data_sample, formula, dvmax, offset=data_sample['TimeOffset'], printResult=False, l1=False)
else:
fitmodel,res = runGLMOLS_severity_formula(data_sample, formula, dvmax, printResult=False, l1=False)
params[trial] = res['Coef.'][3:]
params = params[params.isna().sum(axis=1)!=ntrials]
sorted_perm_stats = np.sort(np.array(params.loc[caseCol,:]))
sorted_perm_stats = sorted_perm_stats - resT.loc[caseCol,'Coef.']
resT.loc[caseCol,'[0.025'] = np.percentile(sorted_perm_stats, 2.5)
resT.loc[caseCol,'0.975]'] = np.percentile(sorted_perm_stats, 97.5)
shapiroP = stats.shapiro(params.loc[caseCol,:])[1]
print('Shapiro test:'); print(shapiroP)
if pvalMethod=='ztest':
pval = ztest(params.loc[caseCol,:], value=resT.loc[caseCol,'Coef.'])
resT.loc[caseCol,pcol] = pval[1]
elif pvalMethod=='count-based':
gt = (params.loc[caseCol,:] >resT.loc[caseCol,'Coef.']).sum()
lt = (params.loc[caseCol,:] <resT.loc[caseCol,'Coef.']).sum()
resT.loc[caseCol,pcol] = np.min([gt,lt]) / (gt+lt)
return(resT)
def runOLS_formula(X2, formula, dv, returnFullResult=False, printResult=True, l1=True):
res = ols(formula, X2).fit()
if printResult==True:
print(res.summary())
t = res.summary2().tables[1][['Coef.','Std.Err.','P>|t|','[0.025','0.975]']]
t = t.fillna(1)
t = t.drop('Intercept')
t = pd.concat([pd.DataFrame({'Coef.':[dv, res.nobs, np.sum(X2['Case']==1)], 'Std.Err.':["","",""], 'P>|t|':["","",""], '[0.025':["","",""], '0.975]':["","",""]},index=['Dep_Variable:','Sample_Count','Sample_Count_Case']),t])
return(t)
def bootstrapOLSregression(X2, formula, dvmax, resT, caseCol, pcol, pvalMethod, ntrials):
params = pd.DataFrame([],index=X2.columns[1:],columns=list(range(0,ntrials)))
for trial in range(ntrials):
data_sample = X2.sample(frac=1, replace=True)
res = runOLS_formula(data_sample, formula, dvmax, printResult=False)
params[trial] = res['Coef.'][3:]
params = params[params.isna().sum(axis=1)!=ntrials]
point_estimate = np.mean(params,axis=1)
bootstrap_standard_error = np.std(params,ddof=1,axis=1)
z_critical = stats.norm.ppf(1 - (1 - 0.95) / 2) # confidence_level = 0.95
resT['[0.025'][3:] = point_estimate - z_critical * bootstrap_standard_error
resT['0.975]'][3:] = point_estimate + z_critical * bootstrap_standard_error
shapiroP = stats.shapiro(params.loc[caseCol,:])[1]
print('Shapiro test:'); print(shapiroP)
if pvalMethod=='ztest':
pval = ztest(params.loc[caseCol,:], value=0)
resT.loc[caseCol,pcol] = pval[1]
elif pvalMethod=='count-based':
gt = (params.loc[caseCol,:] >0).sum()
lt = (params.loc[caseCol,:] <0).sum()
resT.loc[caseCol,pcol] = np.min([gt,lt]) / (gt+lt)
return(resT)
def permutationOLSregression(X2, formula, dvmax, resT, caseCol, pcol, pvalMethod, ntrials):
params = pd.DataFrame([],index=X2.columns[1:],columns=list(range(0,ntrials)))
for trial in range(ntrials):
data_sample = X2.copy()
data_sample[caseCol] = 0
data_sample[caseCol].iloc[sample(range(X2.shape[0]), k=X2[caseCol].sum())] = 1
res = runOLS_formula(data_sample, formula, dvmax, printResult=False)
params[trial] = res['Coef.'][3:]
params = params[params.isna().sum(axis=1)!=ntrials]
sorted_perm_stats = np.sort(np.array(params.loc[caseCol,:]))
sorted_perm_stats = sorted_perm_stats - resT.loc[caseCol,'Coef.']
resT.loc[caseCol,'[0.025'] = np.percentile(sorted_perm_stats, 2.5)
resT.loc[caseCol,'0.975]'] = np.percentile(sorted_perm_stats, 97.5)
shapiroP = stats.shapiro(params.loc[caseCol,:])[1]
print('Shapiro test:'); print(shapiroP)
if pvalMethod=='ztest':
pval = ztest(params.loc[caseCol,:], value=resT.loc[caseCol,'Coef.'])
resT.loc[caseCol,pcol] = pval[1]
elif pvalMethod=='count-based':
gt = (params.loc[caseCol,:] >resT.loc[caseCol,'Coef.']).sum()
lt = (params.loc[caseCol,:] <resT.loc[caseCol,'Coef.']).sum()
resT.loc[caseCol,pcol] = np.min([gt,lt]) / (gt+lt)
return(resT)
def pooledSD(d1, d2):
n1, n2 = len(d1), len(d2)
s1, s2 = np.var(d1, ddof=1), np.var(d2, ddof=1)
return(np.sqrt(((n1 - 1) * s1 + (n2 - 1) * s2) / (n1 + n2 - 2)))
def permutationGroupMeanOLSregression(X2, formula, dvmax, resT, caseCol, pcol, pvalMethod, ntrials):
resT.iloc[3:,:] = np.nan
moT = ols(formula, X2).fit()
temp = moT.predict(X2)
print('Residual normality Shapiro test:'); print(stats.shapiro(moT.resid)[1])
resT.loc[caseCol,'Coef.'] = temp[X2[caseCol]==1].mean() - temp[X2[caseCol]==0].mean()
resT.loc[caseCol,'Std.Err.'] = pooledSD(temp[X2[caseCol]==1], temp[X2[caseCol]==0]) # pooledSD
print('Cohens d:')
print(resT.loc[caseCol,'Coef.'] / resT.loc[caseCol,'Std.Err.']) # cohen's d
params = pd.DataFrame([],index=X2.columns[1:],columns=list(range(0,ntrials)))
for trial in range(ntrials):
data_sample = X2.copy()
data_sample[caseCol] = 0
data_sample[caseCol].iloc[sample(range(X2.shape[0]), k=X2[caseCol].sum())] = 1
temp = ols(formula, data_sample).fit().predict(data_sample)
params.loc[caseCol,trial] = temp[data_sample[caseCol]==1].mean() - temp[data_sample[caseCol]==0].mean()
params = params[params.isna().sum(axis=1)!=ntrials]
params = params.dropna(axis=1)
z_critical = stats.norm.ppf(1 - (1 - 0.95) / 2) # confidence_level = 0.95
resT.loc[caseCol,'[0.025'] = resT.loc[caseCol,'Coef.'] - z_critical * resT.loc[caseCol,'Std.Err.']
resT.loc[caseCol,'0.975]'] = resT.loc[caseCol,'Coef.'] + z_critical * resT.loc[caseCol,'Std.Err.']
shapiroP = stats.shapiro(params.loc[caseCol,:])[1]
print('Shapiro test:'); print(shapiroP)
if pvalMethod=='ztest':
pval = ztest(params.loc[caseCol,:], value=resT.loc[caseCol,'Coef.'])
resT.loc[caseCol,pcol] = pval[1]
elif pvalMethod=='count-based':
gt = (params.loc[caseCol,:] >resT.loc[caseCol,'Coef.']).sum()
lt = (params.loc[caseCol,:] <resT.loc[caseCol,'Coef.']).sum()
resT.loc[caseCol,pcol] = np.min([gt,lt]) / (gt+lt)
return(resT)
def exportStats(X):
x1 = X[X['Case']==1]
x3 = X[X['Case']==0]
num = x1.shape[0], x3.shape[0], "" # count
temp = [num]; temp_index = ['Sample #']
if 'Age' in x1.columns:
age = str(int(np.round(x1['Age'].median(),0))) + " IQR[" + str(int(np.round(np.nanpercentile(x1['Age'],25),0))) +" "+ str(int(np.round(np.nanpercentile(x1['Age'],75),0))) +"]", str(int(np.round(x3['Age'].median(),0))) + " IQR[" + str(int(np.round(np.nanpercentile(x3['Age'],25),0))) +" "+ str(int(np.round(np.nanpercentile(x3['Age'],75),0))) +"]", (rank_compare_2indep(x1['Age'], x3['Age'])).pvalue
temp.append(age); temp_index.append('Age (years)')
if 'Age1' in x1.columns:
age = str(np.round(x1['Age1'].mean(),2)) + " \xb1 " + str(np.round(x1['Age1'].std(),2)), str(np.round(x3['Age1'].mean(),2)) + " \xb1 " + str(np.round(x3['Age1'].std(),2)), (rank_compare_2indep(x1['Age1'], x3['Age1'])).pvalue
temp.append(age); temp_index.append('Age1 (years)')
if 'Age2' in x1.columns:
age = str(np.round(x1['Age2'].mean(),2)) + " \xb1 " + str(np.round(x1['Age2'].std(),2)), str(np.round(x3['Age2'].mean(),2)) + " \xb1 " + str(np.round(x3['Age2'].std(),2)), (rank_compare_2indep(x1['Age2'], x3['Age2'])).pvalue
temp.append(age); temp_index.append('Age2 (years)')
if 'AgeDifference' in x1.columns:
age = str(np.round(x1['AgeDifference'].mean(),2)) + " \xb1 " + str(np.round(x1['AgeDifference'].std(),2)), str(np.round(x3['AgeDifference'].mean(),2)) + " \xb1 " + str(np.round(x3['AgeDifference'].std(),2)), (rank_compare_2indep(x1['AgeDifference'], x3['AgeDifference'])).pvalue
temp.append(age); temp_index.append('AgeDifference (years)')
if 'Squared_AgeDifference' in x1.columns:
age = str(np.round(x1['Squared_AgeDifference'].mean(),2)) + " \xb1 " + str(np.round(x1['Squared_AgeDifference'].std(),2)), str(np.round(x3['Squared_AgeDifference'].mean(),2)) + " \xb1 " + str(np.round(x3['Squared_AgeDifference'].std(),2)), (rank_compare_2indep(x1['Squared_AgeDifference'], x3['Squared_AgeDifference'])).pvalue
temp.append(age); temp_index.append('Squared_AgeDifference (years)')
if 'Sex' in x1.columns:
sex = str((x1['Sex']==0).sum()) + ' [' + str(int(np.round((x1['Sex']==0).mean()*100))) + '%]', str((x3['Sex']==0).sum()) + ' [' + str(int(np.round((x3['Sex']==0).mean()*100))) + '%]', (rank_compare_2indep(x1['Sex'], x3['Sex'])).pvalue
temp.append(sex); temp_index.append('Sex [female] (%)')
if 'Ethnicity' in x1.columns:
race = str((x1['Ethnicity']==0).sum()) + ' [' + str(int(np.round((x1['Ethnicity']==0).mean()*100))) + '%]', str((x3['Ethnicity']==0).sum()) + ' [' + str(int(np.round((x3['Ethnicity']==0).mean()*100))) + '%]', (rank_compare_2indep(x1['Ethnicity'], x3['Ethnicity'])).pvalue
temp.append(race); temp_index.append('Ethnicity [white] (%)')
if 'Education' in x1.columns:
edu = str((x1['Education']==1).sum()) + ' [' + str(int(np.round((x1['Education']==1).mean()*100))) + '%]', str((x3['Education']==1).sum()) + ' [' + str(int(np.round((x3['Education']==1).mean()*100))) + '%]', (rank_compare_2indep(x1['Education'], x3['Education'])).pvalue
temp.append(edu); temp_index.append('Education [college degree] (%)')
if 'Ever_smoked' in x1.columns:
smoke = str((x1['Ever_smoked']==1).sum()) + ' [' + str(int(np.round((x1['Ever_smoked']==1).mean()*100))) + '%]', str((x3['Ever_smoked']==1).sum()) + ' [' + str(int(np.round((x3['Ever_smoked']==1).mean()*100))) + '%]', (rank_compare_2indep(x1['Ever_smoked'], x3['Ever_smoked'])).pvalue
temp.append(smoke); temp_index.append('Ever smoked (%)')
if 'BMI' in x1.columns:
bmi = str(np.round(x1['BMI'].median(),2)) + " IQR[" + str(np.round(np.nanpercentile(x1['BMI'],25),2)) +" "+ str(np.round(np.nanpercentile(x1['BMI'],75),2)) +"]", str(np.round(x3['BMI'].median(),2)) + " IQR[" + str(np.round(np.nanpercentile(x3['BMI'],25),2)) +" "+ str(np.round(np.nanpercentile(x3['BMI'],75),2)) +"]", (rank_compare_2indep(x1['BMI'], x3['BMI'])).pvalue
temp.append(bmi); temp_index.append('BMI')
if '1707-0.0' in x1.columns:
handedness = str((x1['1707-0.0']==1).sum()) + ' [' + str(int(np.round((x1['1707-0.0']==1).mean()*100))) + '%]', str((x3['1707-0.0']==1).sum()) + ' [' + str(int(np.round((x3['1707-0.0']==1).mean()*100))) + '%]', (rank_compare_2indep(x1['1707-0.0'], x3['1707-0.0'])).pvalue
temp.append(handedness); temp_index.append('Handedness [right] (%)')
if 'Self_diabeticCondition' in x1.columns:
diabet = str((x1['Self_diabeticCondition']==1).sum()) + ' [' + str(float(np.round((x1['Self_diabeticCondition']==1).mean()*100,3))) + '%]', str((x3['Self_diabeticCondition']==1).sum()) + ' [' + str(float(np.round((x3['Self_diabeticCondition']==1).mean()*100,3))) + '%]', (rank_compare_2indep(x1['Self_diabeticCondition'], x3['Self_diabeticCondition'])).pvalue
temp.append(diabet); temp_index.append('Self-reported diabetic condition (%)')
if '24104-2.0' in x1.columns:
LV_cardiac_output = str(np.round(x1['24104-2.0'].median(),2)) + " IQR[" + str(np.round(np.nanpercentile(x1['24104-2.0'],25),2)) +" "+ str(np.round(np.nanpercentile(x1['24104-2.0'],75),2)) +"]", str(np.round(x3['24104-2.0'].median(),2)) + " IQR[" + str(np.round(np.nanpercentile(x3['24104-2.0'],25),2)) +" "+ str(np.round(np.nanpercentile(x3['24104-2.0'],75),2)) +"]", (rank_compare_2indep(x1['24104-2.0'], x3['24104-2.0'])).pvalue
temp.append(LV_cardiac_output); temp_index.append('LV cardiac output')
if '22424-2.0' in x1.columns:
Cardiac_output = str(np.round(x1['22424-2.0'].median(),2)) + " IQR[" + str(np.round(np.nanpercentile(x1['22424-2.0'],25),2)) +" "+ str(np.round(np.nanpercentile(x1['22424-2.0'],75),2)) +"]", str(np.round(x3['22424-2.0'].median(),2)) + " IQR[" + str(np.round(np.nanpercentile(x3['22424-2.0'],25),2)) +" "+ str(np.round(np.nanpercentile(x3['22424-2.0'],75),2)) +"]", (rank_compare_2indep(x1['22424-2.0'], x3['22424-2.0'])).pvalue
temp.append(Cardiac_output); temp_index.append('Cardiac output')
if '12336-2.0' in x1.columns:
Ventricular_rate = str(np.round(x1['12336-2.0'].median(),2)) + " IQR[" + str(np.round(np.nanpercentile(x1['12336-2.0'],25),2)) +" "+ str(np.round(np.nanpercentile(x1['12336-2.0'],75),2)) +"]", str(np.round(x3['12336-2.0'].median(),2)) + " IQR[" + str(np.round(np.nanpercentile(x3['12336-2.0'],25),2)) +" "+ str(np.round(np.nanpercentile(x3['12336-2.0'],75),2)) +"]", (rank_compare_2indep(x1['12336-2.0'], x3['12336-2.0'])).pvalue
temp.append(Ventricular_rate); temp_index.append('Ventricular rate')
if '12677-2.0' in x1.columns:
Central_systolic_blood_pressure_during_PWA = str(np.round(x1['12677-2.0'].median(),2)) + " IQR[" + str(np.round(np.nanpercentile(x1['12677-2.0'],25),2)) +" "+ str(np.round(np.nanpercentile(x1['12677-2.0'],75),2)) +"]", str(np.round(x3['12677-2.0'].median(),2)) + " IQR[" + str(np.round(np.nanpercentile(x3['12677-2.0'],25),2)) +" "+ str(np.round(np.nanpercentile(x3['12677-2.0'],75),2)) +"]", (rank_compare_2indep(x1['12677-2.0'], x3['12677-2.0'])).pvalue
temp.append(Central_systolic_blood_pressure_during_PWA); temp_index.append('Central systolic blood pressure during PWA')
T = pd.DataFrame(temp, columns=['Case','Control','P-value'], index=temp_index)
return(T)
def calcComorbidPercentage(df_idf_toE, cate, instance):
D = df_idf_toE.filter(regex=("^41270-"))
D2 = df_idf_toE.filter(regex=("^41271-"))
D3 = df_idf_toE.filter(regex=("^41200-"))
D4 = df_idf_toE.filter(regex=("^41202-"))
D5 = df_idf_toE.filter(regex=("^41203-")) # 41204 does not have time
D6 = df_idf_toE.filter(regex=("^41272-"))
D7 = df_idf_toE.filter(regex=("^41256-"))
# days
Da = df_idf_toE.filter(regex=("^41280-"))
Da2 = df_idf_toE.filter(regex=("^41281-"))
Da3 = df_idf_toE.filter(regex=("^41260-"))
Da4 = df_idf_toE.filter(regex=("^41262-"))
Da5 = df_idf_toE.filter(regex=("^41263-"))
Da6 = df_idf_toE.filter(regex=("^41282-"))
Da7 = df_idf_toE.filter(regex=("^41257-"))
T = pd.DataFrame([],index=df_idf_toE.index, columns=list(cate.keys()))
if instance=='All':
for disease, l in cate.items():
temp = D.apply(lambda x: x.isin(l).sum(),axis=1)
T[disease] = temp
else:
dataDate2 = df_idf_toE['53'+instance].astype('datetime64[ns]')
I = Da.apply(lambda x: (dataDate2 - x.astype('datetime64[ns]')).dt.days >=0,axis=0)
I.columns = D.columns
D.sub = D[I]
for disease, l in cate.items():
temp = D.sub.apply(lambda x: x.isin(l).sum(),axis=1)
T[disease] = temp
co = (T >0).sum(axis=0)
res = co.astype(str) + " [" + (np.round(co / T.shape[0]*100,1)).astype(str) + "%]"
return(res)
def calcPvalueComorbids(res, x1count, x3count):
x1 = pd.DataFrame(0, index=range(x1count), columns=res.index)
x3 = pd.DataFrame(0, index=range(x3count), columns=res.index)
P = pd.DataFrame(pd.NA,index=res.index,columns=['P-value'])
for i in range(res.shape[0]):
n = int(re.sub("^(\d+) .+$","\\1",res.iloc[i,0]))
x1.iloc[0:n,i] = 1
n = int(re.sub("^(\d+) .+$","\\1",res.iloc[i,1]))
x3.iloc[0:n,i] = 1
P.iloc[i,0] = rank_compare_2indep(x1.iloc[:,i], x3.iloc[:,i]).pvalue
P = pd.concat([res,P],axis=1)
return(P)
def value_count_of_ICDs(df_idf_toE, instance):
icd10table = pd.read_csv("UKBB_ICD10_table.csv",sep="\t",index_col=0)
icd10table = icd10table[['Meaning','Selectable']]
icd9table = pd.read_csv("UKBB_ICD9_table.csv",sep="\t",index_col=0)
icd9table = icd9table[['Meaning','Selectable']]
D = df_idf_toE.filter(regex=("^41270-"))
D2 = df_idf_toE.filter(regex=("^41271-"))
Da = df_idf_toE.filter(regex=("^41280-"))
Da2 = df_idf_toE.filter(regex=("^41281-"))
if instance=='All':
icd10 = pd.DataFrame(D[I].apply(pd.Series.value_counts).sum(axis=1).astype(int).sort_values(ascending=False),columns=['Count'])
icd10['Percent'] = icd10 / df_idf_toE.shape[0] * 100
icd10 = pd.merge(icd10, icd10table, left_index=True, right_index=True, how='left')
icd9 = pd.DataFrame(D2[I].apply(pd.Series.value_counts).sum(axis=1).astype(int).sort_values(ascending=False),columns=['Count'])
icd9['Percent'] = icd9 / df_idf_toE.shape[0] * 100
icd9 = pd.merge(icd9, icd9table, left_index=True, right_index=True, how='left')
icd9 = icd9.sort_values(by=['Count'], ascending=False)
return icd10, icd9
else:
dataDate2 = df_idf_toE['53'+instance].astype('datetime64[ns]')
# ICD10
I = Da.apply(lambda x: (dataDate2 - x.astype('datetime64[ns]')).dt.days >=0,axis=0)
I.columns = D.columns
icd10 = pd.DataFrame(D[I].apply(pd.Series.value_counts).sum(axis=1).astype(int).sort_values(ascending=False),columns=['Count'])
icd10['Percent'] = icd10 / df_idf_toE.shape[0] * 100
icd10 = pd.merge(icd10, icd10table, left_index=True, right_index=True, how='left')
# ICD9
I = Da2.apply(lambda x: (dataDate2 - x.astype('datetime64[ns]')).dt.days >=0,axis=0)
I.columns = D2.columns
D2 = D2.astype(str)
D2 = D2.apply(lambda x: x.str.replace("^(.+)\.0","\\1",regex=True), axis=0)
icd9 = pd.DataFrame(D2[I].apply(pd.Series.value_counts).sum(axis=1).astype(int).sort_values(ascending=False),columns=['Count'])
icd9['Percent'] = icd9 / df_idf_toE.shape[0] * 100
icd9 = pd.merge(icd9, icd9table, left_index=True, right_index=True, how='left')
icd9 = icd9.sort_values(by=['Count'], ascending=False)
return icd10, icd9
def snsPlotEachDF_multipage(X2, dv, yvar, xvar, yname, xname, pdf, ylabfontsize, xlabfontsize, ytickfontsize, xtickfontsize, titlefontsize):
bp = sns.boxplot(y=yvar,x=xvar,data=X2, order=[1,0])
plt.xlabel(xname,fontsize=xlabfontsize)
plt.ylabel(yname,fontsize=ylabfontsize)
plt.xticks(fontsize=xtickfontsize)
plt.yticks(fontsize=ytickfontsize)
plt.title(dv, fontsize=titlefontsize)
#bp.axhline(y=0, linewidth=2, color='black', ls=':')
fig = bp.get_figure()
pdf.savefig(fig)
fig.clf()
def missingValueImputeKNN(Before_imputation, cov):
imputer = KNNImputer(n_neighbors=2)
After_imputation = imputer.fit_transform(Before_imputation[cov])
After_imputation = pd.DataFrame(After_imputation,index=Before_imputation.index,columns=cov)
Before_imputation[cov] = After_imputation
return(Before_imputation)
def runBatchRegression_caseControl(df_idf_1, df_idf_3, maps, covariates, caseCol='Case',suffix='suffix',comorbidM=None, booleanizeComorbid=True, removeKinship=True, onlyWhiteBritish=True, instance=None, cohortStats=False, DL=None, focusedDFs=None, output_dir='.', outlierStdv=None, normalizeByTIV=None, DFtype=None, paraEst=None, postfix=None, imputeMissingValues=None):
print("paraEst: "+str(paraEst))
df_idf_1 = filterData(df_idf_1, kinship=removeKinship, onlyWhiteBritish=onlyWhiteBritish)
df_idf_3 = filterData(df_idf_3, kinship=removeKinship, onlyWhiteBritish=onlyWhiteBritish)
if focusedDFs!=None:
df_idf_1 = df_idf_1[(df_idf_1[focusedDFs].isna()).sum(axis=1)!=df_idf_1[focusedDFs].shape[1]]
df_idf_3 = df_idf_3[(df_idf_3[focusedDFs].isna()).sum(axis=1)!=df_idf_3[focusedDFs].shape[1]]
if covariates['UKBB'][0]=='None':
pe_X = pd.DataFrame(1, index=df_idf_1.index, columns=[caseCol])
ctrl_X = pd.DataFrame(0, index=df_idf_3.index, columns=[caseCol])
elif normalizeByTIV:
pe_X = df_idf_1[covariates['UKBB']+['26521-2.0']]
ctrl_X = df_idf_3[covariates['UKBB']+['26521-2.0']]
pe_X.loc[:,caseCol] = 1
ctrl_X.loc[:,caseCol] = 0
else:
pe_X = df_idf_1[covariates['UKBB']]
ctrl_X = df_idf_3[covariates['UKBB']]
pe_X.loc[:,caseCol] = 1
ctrl_X.loc[:,caseCol] = 0
display(pe_X)
display(ctrl_X)
#
X = pd.concat([pe_X, ctrl_X], axis=0)
##############################
if 'UKBB' in covariates:
if '6138-0.0' in covariates['UKBB']:
X = booleanizeEducation(X)
if '21000-0.0' in covariates['UKBB']:
X = booleanizeEthnicity(X)
if onlyWhiteBritish:
X = X[X['21000-0.0']==0]; X = X.drop(['21000-0.0'],axis=1)
if '738-0.0' in covariates['UKBB']:
mea = X.loc[X['738-0.0'] >= 0,'738-0.0'].mean()
X.loc[X['738-0.0'] < 0,'738-0.0'] = mea
X.loc[X['738-0.0'].isna(),'738-0.0'] = mea
if '20002-0.0' in covariates['UKBB']:
X = booleanizeSelfReportedDF20002(X, [1220, 1222, 1223],'Self_diabeticCondition'); X = X.drop(['20002-0.0'],axis=1)
if '20160-0.0' in covariates['UKBB']:
X['20160-0.0'].fillna(0, inplace=True) # all missing values as 0 for smoking
if '23104-0.0' in covariates['UKBB']:
X['23104-0.0'].fillna(X['23104-0.0'].mean(), inplace=True) # BMI imputation by mean
if '30760-0.0' in covariates['UKBB']:
X['30760-0.0'].fillna(X['30760-0.0'].mean(), inplace=True) # HDL cholesterol imputation
if '30780-0.0' in covariates['UKBB']:
X['30780-0.0'].fillna(X['30780-0.0'].mean(), inplace=True) # LDL direct imputation
##############################
X.rename({'31-0.0':'Sex', '6138-0.0':'Education', '21000-0.0':'Ethnicity', '23104-0.0':'BMI', '21022-0.0':'Age', '20160-0.0':'Ever_smoked', '30760-0.0':'HDL_cholesterol','30780-0.0':'LDL_direct','738-0.0':'Income','26521-2.0':'TIV'}, axis=1, inplace=True)
X = X.rename(columns=lambda x: re.sub('22009-0.','PC',x))
outfile_suffix = ''
if comorbidM is not None:
if 'comorbid' in covariates:
if booleanizeComorbid:
comorbidM = (comorbidM >0).astype(int)
X = pd.merge(X, comorbidM[covariates['comorbid']], how='inner', right_index=True, left_index=True)
if imputeMissingValues:
X = missingValueImputeKNN(X, X.columns.tolist())
X = X.dropna()
toExport = ['24486-2.0','24485-2.0','25781-2.0','27230-2.0','27199-2.0','27195-2.0','27188-2.0','27268-2.0','26557-2.0','26556-2.0','26528-2.0','26587-2.0','25895-2.0','25909-2.0','25893-2.0','25903-2.0','25897-2.0','25906-2.0','25900-2.0','25912-2.0','25915-2.0','25918-2.0','25786-2.0','24104-2.0','22424-2.0','12677-2.0','1707-0.0','12336-2.0'] # white matter, left cerebellum, SFG
if set(toExport).issubset(set(df_idf_1.columns)):
df_idf = pd.concat([df_idf_1, df_idf_3],axis=0)
g = pd.merge(X,df_idf[toExport],how='inner',left_index=True,right_index=True)
g.to_csv('DFs_WMH_white_LCerebellumSFG_'+suffix+'.csv',sep=",",quoting=2,index=True,header=True)
else:
g = X.copy()
toExport = ['Left-Cerebellum-Cortex__ctx-lh-superiorfrontal-2.0','ctx-lh-superiorfrontal__ctx-lh-superiorparietal-2.0'] # white matter, left cerebellum, SFG
if set(toExport).issubset(set(df_idf_1.columns)):
df_idf = pd.concat([df_idf_1, df_idf_3],axis=0)
g = pd.merge(X,df_idf[toExport],how='inner',left_index=True,right_index=True)
g.to_csv('DFs_eAUC_LCerebellumSFG_'+suffix+'.csv',sep=",",quoting=2,index=True,header=True)
else:
g = X.copy()
display(g)
T = exportStats(g)
if DL!=None: #
res = pd.concat([calcComorbidPercentage(df_idf_1, DL, instance), calcComorbidPercentage(df_idf_3, DL, instance)], axis=1)
res = calcPvalueComorbids(res, df_idf_1.shape[0], df_idf_3.shape[0])
res.columns = T.columns
T = pd.concat([T, res])
suffix2 = re.sub("^(.+_rfMRI_.+)_(abs|neg|pos).+?$","\\1",suffix)
T.to_csv(output_dir+'/sigStats_DFs-Case_LR'+outfile_suffix+'__'+DFtype+"_"+suffix2+'.csv',sep="\t",header=True,index=True)
icd10,icd9 = value_count_of_ICDs(df_idf_1, instance)
icd10.to_csv(output_dir+'/sigICD10_DFs-Case_LR'+outfile_suffix+'__'+DFtype+"_"+suffix2+'.csv',sep="\t",header=False,index=True)
icd9.to_csv(output_dir+'/sigICD9_DFs-Case_LR'+outfile_suffix+'__'+DFtype+"_"+suffix2+'.csv',sep="\t",header=False,index=True)
icd10,icd9 = value_count_of_ICDs(df_idf_3, instance)
icd10.to_csv(output_dir+'/sigICD10_DFs-Case_LR'+outfile_suffix+'__'+DFtype+"_"+'control'+'.csv',sep="\t",header=False,index=True)
icd9.to_csv(output_dir+'/sigICD9_DFs-Case_LR'+outfile_suffix+'__'+DFtype+"_"+'control'+'.csv',sep="\t",header=False,index=True)
U = calcPHBreakdoown(X[X['Case']==1], 2)
U.to_csv(output_dir+'/sigStats_DFs-Case_LR_PHbreakdown'+outfile_suffix+'__'+DFtype+"_"+suffix2+'.csv',sep="\t",header=True,index=True)
#
df_idf = pd.concat([df_idf_1, df_idf_3],axis=0)
df_idf = df_idf.loc[X.index,:]
pd.DataFrame(df_idf_1.index).to_csv(output_dir+'/eid_df_idf_1'+outfile_suffix+'__'+DFtype+"_"+suffix2+'.csv',sep=",",header=True,index=False,quoting=1)
pd.DataFrame(df_idf_3.index).to_csv(output_dir+'/eid_df_idf_3'+outfile_suffix+'__'+DFtype+"_"+suffix2+'.csv',sep=",",header=True,index=False,quoting=1)
#
if not df_idf.index.identical(X.index):
print('index different between df_idf and X'); exit
L = maps['data_type'].unique()
tested = []
for l in L:
print(l)
Ys = maps.loc[maps['data_type']==l,'DF']
LR = pd.DataFrame()
sub_tested = []
LRtype='DFs-Case_LR'
outpdf = output_dir+"/sigBox_"+LRtype+outfile_suffix+'__'+l+"_"+suffix+postfix+'.pdf'
pdf = matplotlib.backends.backend_pdf.PdfPages(outpdf)
for dv in Ys:
LoR = False
dvmax = str(dv)+instance
if not dvmax in df_idf_1.columns:
print(dvmax," does not exist")
else:
print(dvmax)
# remove -1.0 and -3.0 ("do not know", "prefer not to answer" )
df_idf_1.loc[(df_idf_1[dvmax]==-1) | (df_idf_1[dvmax]==-3),dvmax] = float('NaN')
df_idf_3.loc[(df_idf_3[dvmax]==-1) | (df_idf_3[dvmax]==-3),dvmax] = float('NaN')
df_idf.loc[(df_idf[dvmax]==-1) | (df_idf[dvmax]==-3),dvmax] = float('NaN')
if (len(df_idf_1[dvmax].unique())==1):
print(dvmax," is constant")
elif (len(df_idf[dvmax].unique()) < 4) and (len(df_idf[dvmax].unique()) > 2):
print(dvmax," has less than 4 unique values")
elif (pd.isna(df_idf[dvmax].unique()).sum() >0) and (len(df_idf[dvmax].unique())==2):
print(dvmax," has nan and one values")
elif (isinstance(df_idf[dvmax].unique()[1],str)):
print(dvmax," has non-numeric value")
elif (df_idf[dvmax].dtype=='O'):
print(dvmax," has object value")
else:
sub_tested.append(dv)
df_idf[caseCol] = X[caseCol]
snsPlotEachDF_multipage(df_idf, dv=dvmax, yvar=dvmax, xvar=caseCol, yname='Measurement value', xname='Case', pdf=pdf, ylabfontsize=12, xlabfontsize=12, ytickfontsize=12, xtickfontsize=12, titlefontsize=18)
if (len(df_idf[dvmax].unique())==2) and (np.sum(np.sort(df_idf[dvmax].unique())==np.sort(np.array([0,1])))==2):
LoR = True
print(dvmax," logistic regression applied")
else:
print(dvmax, " has two or more non-binary variables" )
tempX = X
if (normalizeByTIV) & (not l in ['Brain_MRI','aparc+aseg']): # exclude Abdominal, heart, DEXA
tempX = X.drop(['TIV'],axis=1)
formula = 'IDF' + " ~ " + ' + '.join([re.sub("^.+\s(.+)$","\\1",s) for s in tempX.columns.tolist()])
X2 = pd.concat([df_idf[dvmax],X],axis=1)
X2.rename({dvmax:'IDF'}, axis=1, inplace=True)
X2.columns = [re.sub("^.+\s(.+)$","\\1",s) for s in X2.columns.tolist()]
X2 = X2.dropna()
if outlierStdv!=100: ## outlier removal
high = df_idf[dvmax].mean() + outlierStdv * df_idf[dvmax].std()
low = df_idf[dvmax].mean() - outlierStdv * df_idf[dvmax].std()
I = (df_idf[dvmax] <= high) & (df_idf[dvmax] >= low)
if LoR: # logistic
resT = runLogit_caseControl_formula(X2[I], formula, dvmax, printResult=False, l1=False)
if resT.empty:
continue
resT.rename({'P>|z|':'P>|t|'},axis=1,inplace=True)
if paraEst=='permutation':
LR = pd.concat([LR,permutationLoRegression_ntrials(X2[I], formula, dvmax, resT, caseCol, pcol='P>|t|', pvalMethod='count-based', ntrials=10000)], axis=1)
elif paraEst=='bootstrapping':
LR = pd.concat([LR,bootstrapLoRegression(X2[I], formula, dvmax, resT, caseCol, pcol='P>|t|', pvalMethod='count-based', ntrials=10000)], axis=1)
elif paraEst=='meanDiffClassPermutation':
LR = pd.concat([LR, permutationGroupMeanLoRegression(X2[I], formula, dvmax, resT, caseCol, pcol='P>|t|', pvalMethod='count-based', ntrials=10000)], axis=1)
else:
LR = pd.concat([LR, resT], axis=1)
else: # OLS
resT = runOLS_formula(X2[I], formula, dvmax, printResult=True)
if paraEst=='permutation': # permutation
LR = pd.concat([LR,permutationOLSregression(X2[I], formula, dvmax, resT, caseCol, pcol='P>|t|', pvalMethod='count-based', ntrials=10000)], axis=1)
elif paraEst=='bootstrapping':
LR = pd.concat([LR,bootstrapOLSregression(X2[I], formula, dvmax, resT, caseCol, pcol='P>|t|', pvalMethod='count-based', ntrials=10000)], axis=1)
elif paraEst=='meanDiffClassPermutation':
LR = pd.concat([LR, permutationGroupMeanOLSregression(X2[I], formula, dvmax, resT, caseCol, pcol='P>|t|', pvalMethod='count-based', ntrials=10000)], axis=1)
else:
LR = pd.concat([LR,resT], axis=1)
else: ## No outlier removal
if LoR: # logistic regression
resT = runLogit_caseControl_formula(X2, formula, dvmax, printResult=False, l1=False)
if resT.empty:
continue
resT.rename({'P>|z|':'P>|t|'},axis=1,inplace=True)
if paraEst=='permutation':
LR = pd.concat([LR, permutationLoRegression_ntrials(X2, formula, dvmax, resT, caseCol, pcol='P>|t|', pvalMethod='count-based', ntrials=10000)], axis=1)
elif paraEst=='bootstrapping':
LR = pd.concat([LR, bootstrapLoRegression(X2, formula, dvmax, resT, caseCol, pcol='P>|t|', pvalMethod='count-based', ntrials=10000)], axis=1)
elif paraEst=='meanDiffClassPermutation':
LR = pd.concat([LR, permutationGroupMeanLoRegression(X2, formula, dvmax, resT, caseCol, pcol='P>|t|', pvalMethod='count-based', ntrials=10000)], axis=1)
else:
LR = pd.concat([LR, resT], axis=1)
else: # OLS
resT = runOLS_formula(X2, formula, dvmax, printResult=True)
if paraEst=='permutation': # permutation
LR = pd.concat([LR, permutationOLSregression(X2, formula, dvmax, resT, caseCol, pcol='P>|t|', pvalMethod='count-based', ntrials=10000)], axis=1)
elif paraEst=='bootstrapping':
LR = pd.concat([LR, bootstrapOLSregression(X2, formula, dvmax, resT, caseCol, pcol='P>|t|', pvalMethod='count-based', ntrials=10000)], axis=1)
elif paraEst=='meanDiffClassPermutation':
LR = pd.concat([LR, permutationGroupMeanOLSregression(X2, formula, dvmax, resT, caseCol, pcol='P>|t|', pvalMethod='count-based', ntrials=10000)], axis=1)
else:
LR = pd.concat([LR,resT], axis=1)
pdf.close()
tested.append(sub_tested)
if paraEst in ['bootstrapping','permutation']:
paraEst = "_"+paraEst
else:
paraEst=""
if LR.empty:
print(l, " has no usable variable")
continue
I = np.where(LR.loc[caseCol,'P>|t|'] <= 1.0)[0]
if(len(I)==0):
print(l, " has no significant result")
temp = pd.DataFrame(np.nan, index=maps['DF'], columns=['DF','Coefficient','P-value','Padj','CI_025','CI_975','Sample_Count','Sample_Count_Case'])
temp['DF'] = temp.index
temp.to_csv(output_dir+'/sig2_DFs-Case_LR'+outfile_suffix+'__'+l+"_"+suffix+paraEst+postfix+'.csv',sep=",",header=True,index=False)
continue
lr = LR.iloc[:,np.sort(np.concatenate([I*5,I*5+1,I*5+2,I*5+3,I*5+4]))]
lr.to_csv(output_dir+'/sig_DFs-Case_LR'+outfile_suffix+'__'+l+"_"+suffix+paraEst+postfix+'.csv',sep=",",header=True,index=True)
padj = (ss.multitest.fdrcorrection(LR.loc[caseCol,'P>|t|']))[1]
padjBonf = (ss.multitest.multipletests(LR.loc[caseCol,'P>|t|'], method='bonferroni'))[1]
target_I = int(np.where(LR.index==caseCol)[0])
b = pd.DataFrame([LR.iloc[0,np.sort(np.concatenate([I*5]))].values, LR.iloc[target_I,np.sort(np.concatenate([I*5]))].values, LR.iloc[target_I,np.sort(np.concatenate([I*5+2]))].values, padj[I], padjBonf[I], LR.iloc[target_I,np.sort(np.concatenate([I*5+3]))].values, LR.iloc[target_I,np.sort(np.concatenate([I*5+4]))].values, LR.iloc[1,np.sort(np.concatenate([I*5]))].values, LR.iloc[2,np.sort(np.concatenate([I*5]))].values]).T
b.columns = ['DF','Coefficient','P-value','Padj','Padj_Bonf','CI_025','CI_975','Sample_Count','Sample_Count_Case']
b['DF'+instance] = b['DF'].copy()
b['DF'] = [re.sub("^(.+)-.+$","\\1",i) for i in b['DF'] ]
if b['DF'][0].isnumeric():
b['DF'] = b['DF'].astype('int')
sub1 = df_idf_1.loc[:,df_idf_1.columns.isin(b['DF'+instance])]
sub3 = df_idf_3.loc[:,df_idf_3.columns.isin(b['DF'+instance])]
mean1 = sub1.mean(axis=0); std1=sub1.std(axis=0); min1=sub1.min(axis=0); q5_1=sub1.quantile(0.05, axis=0); q25_1=sub1.quantile(0.25, axis=0); med1=sub1.quantile(0.5, axis=0); q75_1=sub1.quantile(0.75, axis=0); q95_1=sub1.quantile(0.95, axis=0); max1=sub1.max(axis=0)
mean3 = sub3.mean(axis=0); std3=sub3.std(axis=0); min3=sub3.min(axis=0); q5_3=sub3.quantile(0.05, axis=0); q25_3=sub3.quantile(0.25, axis=0); med3=sub3.quantile(0.5, axis=0); q75_3=sub3.quantile(0.75, axis=0); q95_3=sub3.quantile(0.95, axis=0); max3=sub3.max(axis=0)
delta = mean1 - mean3
delta_percent = delta / mean3 *100
ms1 = pd.DataFrame(np.round(mean1,2).astype(str) + " \xb1 " + np.round(std1,2).astype(str), columns=['Case_mean[std]'])
ms3 = pd.DataFrame(np.round(mean3,2).astype(str) + " \xb1 " + np.round(std3,2).astype(str), columns=['Control_mean[std]'])
de = pd.DataFrame(delta, columns=['Case-Control']).astype(str)
de_percent = pd.DataFrame(delta_percent, columns=['Percent_Case-Control']).astype(str)
sss = pd.DataFrame([min1,q5_1,q25_1,med1,q75_1,q95_1,max1,min3,q5_3,q25_3,med3,q75_3,q95_3,max3],index=['Case_min','Case_q5','Case_q25','Case_q50','Case_q75','Case_q95','Case_max','Control_min','Control_q5','Control_q25','Control_q50','Control_q75','Control_q95','Control_max']).astype(str).T
b = pd.merge(b,ms1,how='left',left_on='DF'+instance,right_index=True)
b = pd.merge(b,ms3,how='left',left_on='DF'+instance,right_index=True)
b = pd.merge(b,de,how='left',left_on='DF'+instance,right_index=True)
b = pd.merge(b,de_percent,how='left',left_on='DF'+instance,right_index=True)
b = pd.merge(b,sss,how='left',left_on='DF'+instance,right_index=True)
temp = pd.merge(maps[['DF','description','data_type','Image_type','Sub_classification']], b, left_on='DF',right_on='DF', how='inner')
temp = temp.sort_values('P-value')
temp.to_csv(output_dir+'/sig2_DFs-Case_LR'+outfile_suffix+'__'+l+"_"+suffix+paraEst+postfix+'.csv',sep=",",header=True,index=False)
if l=="Brain_MRI":
temp['Padj_global'] = temp['Padj'].copy()
temp['Padj'] = np.float('NaN')
atlases = temp['Sub_classification'].unique()
for atla in atlases:
I = temp['Sub_classification']==atla
temp.loc[I,'Padj'] = ss.multitest.fdrcorrection(temp[I]['P-value'])[1]
temp = temp.sort_values('P-value')
temp.to_csv(output_dir+'/sig2_DFs-Case_LR'+outfile_suffix+'__'+l+"_"+suffix+paraEst+postfix+'.csv',sep=",",header=True,index=False)
temp = temp[temp['Padj'] <= 0.05]
temp = temp[temp['Sub_classification'].isin(['Freesurfer_a2009s','Freesurfer_ASEG','Freesurfer_BA_exvivo','Freesurfer_DKT','Freesurfer_subsegmentation','Regional_grey_matter_volumes_(FAST)','Subcortical_volumes_(FIRST)'])]
temp = temp.sort_values(['description','Sub_classification'], ascending=[False, True])
temp.to_csv(output_dir+'/significant_sig2_DFs-Case_LR'+outfile_suffix+'__'+l+"_"+suffix+paraEst+postfix+'.csv',sep=",",header=True,index=False)
def runBatchRegression_SNPvsPhenotype_withoutControl(df_idf_1, maps, covariates, caseCol='Case', suffix='suffix', comorbidM=None, booleanizeComorbid=True, removeKinship=False, onlyWhiteBritish=True, instance=None, cohortStats=True, DL=None, focusedDFs=None, normalizeByTIV=None, paraEst=None, output_dir='.', snps=None, postfix=None, gtypes=None, setname=None, snpmodel=None):
print("paraEst: "+str(paraEst))
df_idf_1 = filterData(df_idf_1, kinship=removeKinship, onlyWhiteBritish=onlyWhiteBritish)
if focusedDFs!=None:
df_idf_1 = df_idf_1[(df_idf_1[focusedDFs].isna()).sum(axis=1)!=df_idf_1[focusedDFs].shape[1]]
if covariates['UKBB'][0]=='None':
pass
elif normalizeByTIV:
pe_X = df_idf_1[covariates['UKBB']+['26521-2.0']]
else:
pe_X = df_idf_1[covariates['UKBB']]
display(pe_X)
X = pe_X
##############################
if 'UKBB' in covariates:
if '6138-0.0' in covariates['UKBB']:
X = booleanizeEducation(X)
if '21000-0.0' in covariates['UKBB']:
X = booleanizeEthnicity(X)
if onlyWhiteBritish:
X = X[X['21000-0.0']==0]; X = X.drop(['21000-0.0'],axis=1)
if '738-0.0' in covariates['UKBB']:
mea = X.loc[X['738-0.0'] >= 0,'738-0.0'].mean()
X.loc[X['738-0.0'] < 0,'738-0.0'] = mea
X.loc[X['738-0.0'].isna(),'738-0.0'] = mea
if '20002-0.0' in covariates['UKBB']:
X = booleanizeSelfReportedDF20002(X, [1220, 1222, 1223],'Self_diabeticCondition'); X = X.drop(['20002-0.0'],axis=1)
if '20160-0.0' in covariates['UKBB']:
X['20160-0.0'].fillna(0, inplace=True) # all missing values as 0 for smoking
if '23104-0.0' in covariates['UKBB']:
X['23104-0.0'].fillna(X['23104-0.0'].mean(), inplace=True) # BMI imputation by mean
if '30760-0.0' in covariates['UKBB']:
X['30760-0.0'].fillna(X['30760-0.0'].mean(), inplace=True) # HDL cholesterol imputation
if '30780-0.0' in covariates['UKBB']:
X['30780-0.0'].fillna(X['30780-0.0'].mean(), inplace=True) # LDL direct imputation
##############################
X.rename({'31-0.0':'Sex', '6138-0.0':'Education', '21000-0.0':'Ethnicity', '23104-0.0':'BMI', '21022-0.0':'Age', '20160-0.0':'Ever_smoked', '30760-0.0':'HDL_cholesterol','30780-0.0':'LDL_direct','738-0.0':'Income','26521-2.0':'TIV'}, axis=1, inplace=True)
X = X.rename(columns=lambda x: re.sub('22009-0.','PC',x))
outfile_suffix = ''
if comorbidM is not None:
if 'comorbid' in covariates:
if booleanizeComorbid:
comorbidM = (comorbidM >0).astype(int)
X = pd.merge(X, comorbidM[covariates['comorbid']], how='inner', right_index=True, left_index=True)
X = X.dropna()
display(X)
suffix2 = re.sub("^(.+_rfMRI_.+)_(abs|neg|pos).+?$","\\1",suffix)
df_idf = df_idf_1.copy()
df_idf = df_idf.loc[X.index,:]
df_idf_1.to_csv(output_dir+'/df_idf_1.csv',sep=",",header=True,index=True,quoting=1)
if not df_idf.index.identical(X.index):
print('index different between df_idf and X'); exit
L = maps['data_type'].unique()
tested = []
for l in L:
print(l)
Ys = maps.loc[maps['data_type']==l,'DF']
sub_tested = []
for dv in Ys:
LoR = False
dvmax = str(dv)+instance
print(dvmax)
pcol = 'P>|z|'
LR0 = pd.DataFrame()
if not dvmax in df_idf_1.columns:
print(dvmax," does not exist")
else:
# remove -1.0 and -3.0 ("do not know", "prefer not to answer" )
df_idf_1.loc[(df_idf_1[dvmax]==-1) | (df_idf_1[dvmax]==-3),dvmax] = float('NaN')
df_idf.loc[(df_idf[dvmax]==-1) | (df_idf[dvmax]==-3),dvmax] = float('NaN')
if (len(df_idf_1[dvmax].unique())==1):
print(dvmax," is constant")
elif (pd.isna(df_idf[dvmax].unique()).sum() >0) and (len(df_idf[dvmax].unique())==2):
print(dvmax," has nan and one values")
elif (isinstance(df_idf[dvmax].unique()[1],str)):
print(dvmax," has non-numeric value")
elif (df_idf[dvmax].dtype=='O'):
print(dvmax," has object value")
else:
sub_tested.append(dvmax)
if (len(df_idf[dvmax].unique())==2) and (np.sum(np.sort(df_idf[dvmax].unique())==np.sort(np.array([0,1])))==2):
print(dvmax," logistic regression applied")
elif (len(df_idf[dvmax].unique()) < 2):
print(outcome," has < 1 unique outcome variable")
continue
else:
print(dvmax, " has two or more non-binary variables" )
for eachSNP in gtypes: # regression one by one
print(eachSNP)
if not isinstance(snps, pd.DataFrame):
sys.exit('input snps not defined')
snps2 = snps[eachSNP]
X15 = pd.merge(X,snps2,left_index=True,right_index=True,how='inner')
X15 = X15[~X15[eachSNP].isna()]
X15[eachSNP] = X15[eachSNP].astype(int)
y = df_idf.loc[:,~df_idf.columns.isin(X15.columns)]
y = y[y.index.isin(X15.index)]
y2 = y[dvmax]
y2 = y2[~y2.isna()]
pcol = 'P>|z|'
if snpmodel=='dominant':
X15.loc[X15[eachSNP]==2,eachSNP] = 1
diff_X2 = pd.merge(y2, X15, right_index=True, left_index=True, how='inner')
toExport = ['rs11154337_G']
if set(toExport).issubset(set(diff_X2.columns)):
if snpmodel=='additive':
diff_X2.to_csv('DFs_regre4_SNPs_vs_'+str(dv)+"_"+suffix+'.csv',sep=",",quoting=2,index=True,header=True)
diff_X2.rename({dvmax:'IDF'},axis=1,inplace=True)
diff_X2 = diff_X2.dropna()
diff_X2.columns = [re.sub("^.+\s(.+)$","\\1",s) for s in diff_X2.columns.tolist()]
formula = 'IDF ~ '+' + '.join(diff_X2.columns.tolist())
formula = re.sub("( IDF \+| Case \+)","",formula)
if not normalizeByTIV: # exclude Abdominal, heart, DEXA
formula = formula.replace(" + TIV","")
minimumSample = len(formula.split(' + ')) + 1
if diff_X2.shape[0] < minimumSample:
print("diff_X2.shape[0] < minimumSample")
continue
if (len(diff_X2['IDF'].unique()) < 2):
print(dvmax," has < 1 unique outcome variable")
continue
if LoR: # logistic regression
fitmodel,resT = runGLMLogit_severity_formula(diff_X2, formula, dvmax, printResult=False, l1=False)
else: # OLS
fitmodel,resT = runGLMOLS_severity_formula(diff_X2, formula, dvmax, printResult=False, l1=False)
if resT.empty:
continue
if paraEst=='bootstrapping':
pass
elif paraEst=='meanDiffClassBootstrapping':
pass
else:
pass
effectAlleleNums = np.sum(diff_X2[eachSNP] >0)
totalAlele = diff_X2.shape[0]*2
effectAlleleCount = np.sum(diff_X2[eachSNP])
eaf = pd.DataFrame({'Coef.':[effectAlleleNums, effectAlleleCount / totalAlele, effectAlleleCount, totalAlele], 'Std.Err.':["","","",""], 'P>|z|':["","","",""], '[0.025':["","","",""], '0.975]':["","","",""], 'Coef.sx':["","","",""], 'Coef_sxsy':["","","",""]},index=['eaSampleCountCase','eaf','eacount','totalacount'])
resT = pd.concat([resT[0:3], eaf, resT[3:resT.shape[0]] ], axis=0)
m = pd.DataFrame(pd.concat([resT.loc[['Dep_Variable:'],'Coef.'], maps.loc[maps['DF'].isin([dv]),'description'], resT.loc[eachSNP], resT.loc[['Sample_Count','Sample_Count_Case','eaSampleCountCase','eaf','eacount','totalacount'],'Coef.']], axis=0)).T
m.columns.values[0:2] = ['DF','description']
m.index = [eachSNP]
LR0 = pd.concat([LR0,m], axis=0)
if paraEst in ['bootstrapping','meanDiffClassBootstrapping']:
paraEst2 = "_"+paraEst
else:
paraEst2 = ""
if LR0.empty:
print(l, " has no usable variable")
continue
I = np.where(LR0[pcol] <= 1.0)[0]
if(len(I)==0):
print(setname, " has no significant result")
padj = (ss.multitest.fdrcorrection(LR0[pcol]))[1]
padjBonf = (ss.multitest.multipletests(LR0[pcol], method='bonferroni'))[1]
LR0.insert(5,'Padj',padj)
LR0.insert(6,'padjBonf',padj)
LR0 = LR0.sort_values(pcol)
LR0.index.name='SNP'
LR0.to_csv(output_dir+'/sig2_SNPvsPhenotype'+outfile_suffix+'__'+suffix+paraEst2+postfix+"_"+setname+'_'+snpmodel+"__"+dvmax+'.csv',sep=",",header=True,index=True)
def runBatchSeverityRegression_inputBothOutcomeAndExposureVectors(df_idf_1, maps, covariates, caseCol='Predictor',suffix='suffix',comorbidM=None, booleanizeComorbid=True, removeKinship=True, onlyWhiteBritish=True, instance=None, LRtype=None, cohortStats=False, DL=None, focusedDFs=None, outcomeDescriptions=None, output_dir='.', outlierStdv=None, normalizeByTIV=False, DFtype=None, paraEst=None, postfix=None, timeOffset=False, ntrials=None):
print("paraEst: "+str(paraEst))
df_idf_1 = filterData(df_idf_1, kinship=removeKinship, onlyWhiteBritish=onlyWhiteBritish)
if focusedDFs!=None:
df_idf_1 = df_idf_1[(df_idf_1[focusedDFs].isna()).sum(axis=1)!=df_idf_1[focusedDFs].shape[1]]
if timeOffset:
if not 'TimeOffset' in covariates['UKBB']:
covariates['UKBB'] = covariates['UKBB']+['TimeOffset']
if covariates['UKBB'][0]=='None':
pe_X = pd.DataFrame(1, index=df_idf_1.index, columns=[caseCol])
pass
elif normalizeByTIV:
pe_X = df_idf_1[covariates['UKBB']+['26521-2.0']]
else:
pe_X = df_idf_1[covariates['UKBB']]
display(pe_X)
X = pe_X
#
##############################
if 'UKBB' in covariates:
if '6138-0.0' in covariates['UKBB']:
X = booleanizeEducation(X)
if '21000-0.0' in covariates['UKBB']:
X = booleanizeEthnicity(X)
if onlyWhiteBritish:
X = X[X['21000-0.0']==0]; X = X.drop(['21000-0.0'],axis=1)
if '738-0.0' in covariates['UKBB']:
mea = X.loc[X['738-0.0'] >= 0,'738-0.0'].mean()
X.loc[X['738-0.0'] < 0,'738-0.0'] = mea
X.loc[X['738-0.0'].isna(),'738-0.0'] = mea
if '20002-0.0' in covariates['UKBB']:
X = booleanizeSelfReportedDF20002(X, [1220, 1222, 1223],'Self_diabeticCondition'); X = X.drop(['20002-0.0'],axis=1)
if '20160-0.0' in covariates['UKBB']:
X['20160-0.0'].fillna(0, inplace=True) # all missing values as 0 for smoking
if '23104-0.0' in covariates['UKBB']:
X['23104-0.0'].fillna(X['23104-0.0'].mean(), inplace=True) # BMI imputation by mean
if '30760-0.0' in covariates['UKBB']:
X['30760-0.0'].fillna(X['30760-0.0'].mean(), inplace=True) # HDL cholesterol imputation
if '30780-0.0' in covariates['UKBB']:
X['30780-0.0'].fillna(X['30780-0.0'].mean(), inplace=True) # LDL direct imputation
##############################
X.rename({'31-0.0':'Sex', '6138-0.0':'Education', '21000-0.0':'Ethnicity', '23104-0.0':'BMI', '21022-0.0':'Age', '20160-0.0':'Ever_smoked', '30760-0.0':'HDL_cholesterol','30780-0.0':'LDL_direct','738-0.0':'Income','26521-2.0':'TIV'}, axis=1, inplace=True)
X = X.rename(columns=lambda x: re.sub('22009-0.','PC',x))
outfile_suffix = ''
if comorbidM is not None:
if 'comorbid' in covariates:
if booleanizeComorbid:
comorbidM = (comorbidM >0).astype(int)
X = pd.merge(X, comorbidM[covariates['comorbid']], how='inner', right_index=True, left_index=True)
X = X.dropna()
toExport = ['30870-0.0','30860-0.0','30840-0.0','30820-0.0','30810-0.0','30800-0.0','30790-0.0','30780-0.0','30760-0.0','30740-0.0','30730-0.0','30660-0.0','30690-0.0','30680-0.0','30650-0.0','30640-0.0','30630-0.0','30610-0.0','30620-0.0','30890-0.0','30670-0.0','30700-0.0','30720-0.0','30880-0.0','30750-0.0','30770-0.0','30600-0.0','30830-0.0','30850-0.0','30710-0.0','27230-2.0','27199-2.0','26557-2.0','26556-2.0','26521-2.0'] # white matter, left cerebellum, SFG, TIV
if set(toExport).issubset(set(df_idf_1.columns)):
g = pd.merge(X,df_idf_1[toExport],how='inner',left_index=True,right_index=True)
g.to_csv('DFs_Blood_biochemistry_LCerebellumSFG_'+suffix+'.csv',sep=",",quoting=2,index=True,header=True)
toExport = ['Max_wsumed_cardsupdays','Summed_episode_length','Summed_wsumed_cardsupdays','Summed_All_supdays','Max_wsumed_card_res_supdays','Summed_wsumed_card_res_supdays','Max_All_supdays','Max_orgsupmax','Summed_wsumed_ressupdays','Count_gisupcares','Max_wsumed_ressupdays','LCerebellum__LSFG_z0.675_lower','L_Cerebellum_WM_and_LCerebellum__LSFG_AND_L_SFG_Th_and_LCerebellum__LSFG_z0.675_lower'] # white matter, left cerebellum, SFG
toExport = [s+instance for s in toExport]
if set(toExport).issubset(set(df_idf_1.columns)):
g = pd.merge(X,df_idf_1[toExport],how='inner',left_index=True,right_index=True)
g.to_csv('DFs_regre13_ClinicalOutcomesvsMRI_RI_LCerebellumWMSFG_'+suffix+'.csv',sep=",",quoting=2,index=True,header=True)
display(X)
suffix2 = re.sub("^(.+_rfMRI_.+)_(abs|neg|pos).+?$","\\1",suffix)
df_idf = df_idf_1.copy()
df_idf = df_idf.loc[X.index,:]
pd.DataFrame(df_idf_1.index).to_csv(output_dir+'/eid_df_idf_1'+outfile_suffix+'__'+DFtype+"_"+suffix2+'.csv',sep=",",header=True,index=False,quoting=1)
if not df_idf.index.identical(X.index):
print('index different between df_idf and X'); exit
L = maps['data_type'].unique()
tested = []
for outcome in focusedDFs: # export output file for each outcome
print(outcome)
for l in L:
print(l)
Ys = maps.loc[maps['data_type']==l,'DF']
LR = pd.DataFrame()
sub_tested = []
outpdf = output_dir+"/regressExog_"+LRtype+outfile_suffix+'__'+l+"_"+suffix+postfix+'.pdf'
for dv in Ys:
LoR = False
if instance=='NA':
dvmax = str(dv)
else:
dvmax = str(dv)+instance
if dvmax==outcome:
continue
if not dvmax in df_idf_1.columns:
print(dvmax," does not exist")
else:
print(dvmax)
# remove -1.0 and -3.0 ("do not know", "prefer not to answer" )
df_idf_1.loc[(df_idf_1[dvmax]==-1) | (df_idf_1[dvmax]==-3),dvmax] = float('NaN')
df_idf.loc[(df_idf[dvmax]==-1) | (df_idf[dvmax]==-3),dvmax] = float('NaN')
if (len(df_idf_1[dvmax].unique())==1):
print(dvmax," is constant")
elif (pd.isna(df_idf[dvmax].unique()).sum() >0) and (len(df_idf[dvmax].unique())==2):
print(dvmax," has nan and one values")
elif (isinstance(df_idf[dvmax].unique()[1],str)):
print(dvmax," has non-numeric value")
elif (df_idf[dvmax].dtype=='O'):
print(dvmax," has object value")
else:
sub_tested.append(dv)
X[caseCol] = df_idf[dvmax]
if (len(df_idf[outcome].unique())==2) and (np.sum(np.sort(df_idf[outcome].unique())==np.sort(np.array([0,1])))==2):
LoR = True
print(outcome," logistic regression applied")
elif (len(df_idf[outcome].unique()) < 2):
print(outcome," has < 1 unique outcome variable")
continue
else:
print(outcome, " has two or more non-binary variables" )
X2 = X.copy()
X2.columns = [re.sub("^.+\s(.+)$","\\1",s) for s in X2.columns.tolist()]
X2['IDF'] = df_idf[outcome]
X2 = X2.dropna()
if X2.shape[0] <10:
print('too few input')
continue
if timeOffset:
ax = X2['TimeOffset'].hist() # s is an instance of Series
print(shapiro(X2['TimeOffset']))
X2['TimeOffset'] = stats.zscore(X2['TimeOffset'])
X2['Predictor'] = stats.zscore(X2['Predictor'])
if normalizeByTIV:
X2['TIV'] = stats.zscore(X2['TIV'])
formula = 'IDF' + " ~ " + ' + '.join([re.sub("^.+\s(.+)$","\\1",s) for s in X.drop(['TimeOffset'],axis=1).columns.tolist()]) + " - 1"
X3 = X2.copy()
X3['Predictor_over_TimeOffset'] = X2['Predictor'] / X2['TimeOffset']
X3['Predictor_minus_TimeOffset'] = X2['Predictor'] - X2['TimeOffset']
else:
formula = 'IDF' + " ~ " + ' + '.join([re.sub("^.+\s(.+)$","\\1",s) for s in X.columns.tolist()])
if not normalizeByTIV: # exclude Abdominal, heart, DEXA
formula = formula.replace(" + TIV","")
if outlierStdv!=100: ## outlier removal
high = df_idf[dvmax].mean() + outlierStdv * df_idf[dvmax].std()
low = df_idf[dvmax].mean() - outlierStdv * df_idf[dvmax].std()
I = (df_idf[dvmax] <= high) & (df_idf[dvmax] >= low)
X2 = X2[I]
if LoR: # logistic regression
if timeOffset:
fitmodel,resT = runGLMLogit_severity_formula(X2, formula, dvmax, offset=X2['TimeOffset'], printResult=False, l1=False)
else:
fitmodel,resT = runGLMLogit_severity_formula(X2, formula, dvmax, printResult=False, l1=False)
if resT.empty:
continue
if paraEst=='permutation':
if timeOffset:
LR = pd.concat([LR, permutationGLMLogitregression(X2, formula, dvmax, resT, caseCol, pcol='P>|z|', pvalMethod='count-based', ntrials=ntrials, timeOffset=timeOffset)], axis=1)
else:
LR = pd.concat([LR, permutationGLMLogitregression(X2, formula, dvmax, resT, caseCol, pcol='P>|z|', pvalMethod='count-based', ntrials=ntrials, timeOffset=False)], axis=1)
elif paraEst=='bootstrapping':
if timeOffset:
LR = pd.concat([LR, bootstrapGLMLogitregression(X2, formula, dvmax, resT, caseCol,pcol='P>|z|', pvalMethod='count-based', ntrials=ntrials, timeOffset=timeOffset)], axis=1)
else:
LR = pd.concat([LR, bootstrapGLMLogitregression(X2, formula, dvmax, resT, caseCol,pcol='P>|z|', pvalMethod='count-based', ntrials=ntrials, timeOffset=False)], axis=1)
else:
LR = pd.concat([LR, resT], axis=1)
else: # OLS
if timeOffset:
fitmodel,resT = runGLMOLS_severity_formula(X2, formula, dvmax, offset=X2['TimeOffset'], printResult=False, l1=False)
else:
fitmodel,resT = runGLMOLS_severity_formula(X2, formula, dvmax, printResult=False, l1=False)
if paraEst=='permutation':
if timeOffset:
LR = pd.concat([LR, permutationGLMOLSregression(X2, formula, dvmax, resT, caseCol, pcol='P>|z|', pvalMethod='count-based', ntrials=ntrials, timeOffset=timeOffset)], axis=1)
else:
LR = pd.concat([LR, permutationGLMOLSregression(X2, formula, dvmax, resT, caseCol, pcol='P>|z|', pvalMethod='count-based', ntrials=ntrials, timeOffset=False)], axis=1)
elif paraEst=='bootstrapping':
if timeOffset:
LR = pd.concat([LR, bootstrapGLMOLSregression(X2, formula, dvmax, resT, caseCol, pcol='P>|z|', pvalMethod='count-based', ntrials=ntrials, timeOffset=timeOffset)], axis=1)
else:
LR = pd.concat([LR, bootstrapGLMOLSregression(X2, formula, dvmax, resT, caseCol, pcol='P>|z|', pvalMethod='count-based', ntrials=ntrials, timeOffset=False)], axis=1)
else:
LR = pd.concat([LR,resT], axis=1)
tested.append(sub_tested)
if paraEst in ['bootstrapping','permutation']:
paraEst2 = "_"+paraEst
else:
paraEst2 = ""
if LR.empty:
print(l, " has no usable variable")
continue
I = np.where(LR.loc[caseCol,'P>|z|'] <= 1.0)[0]
if(len(I)==0):
print(l, " has no significant result")
continue
lr = LR.iloc[:,np.sort(np.concatenate([I*7,I*7+1,I*7+2,I*7+3,I*7+4,I*7+5,I*7+6]))]
padj = (ss.multitest.fdrcorrection(LR.loc[caseCol,'P>|z|']))[1]
padjBonf = (ss.multitest.multipletests(LR.loc[caseCol,'P>|z|'], method='bonferroni'))[1]
target_I = int(np.where(LR.index==caseCol)[0])
b = pd.DataFrame([LR.iloc[0,np.sort(np.concatenate([I*7]))].values, LR.iloc[target_I,np.sort(np.concatenate([I*7]))].values, LR.iloc[target_I,np.sort(np.concatenate([I*7+2]))].values, padj[I], padjBonf[I], LR.iloc[target_I,np.sort(np.concatenate([I*7+3]))].values, LR.iloc[target_I,np.sort(np.concatenate([I*7+4]))].values, LR.iloc[target_I,np.sort(np.concatenate([I*7+5]))].values, LR.iloc[target_I,np.sort(np.concatenate([I*7+6]))].values, LR.iloc[1,np.sort(np.concatenate([I*7]))].values, LR.iloc[2,np.sort(np.concatenate([I*7]))].values]).T
b.columns = ['DF','Coefficient','P-value','Padj','Padj_Bonf','CI_025','CI_975','Coef_sx','Coef_sxsy','Sample_Count','Sample_Count_Case']
b['DF'+instance] = b['DF'].copy()
b['DF'] = [re.sub("^(.+)-\d+\.\d+$","\\1",i) for i in b['DF'] ]
b['DF'] = [re.sub("^(.+)-All$","\\1",i) for i in b['DF'] ]
if b['DF'].str.isnumeric().sum()==b.shape[0]:
b['DF'] = b['DF'].astype('int')
else:
b['DF'] = b['DF'].astype(str)
maps['DF'] = maps['DF'].astype(str)
sub1 = df_idf_1.loc[:,df_idf_1.columns.isin(b['DF'+instance])]
mean1 = sub1.mean(axis=0); std1=sub1.std(axis=0); min1=sub1.min(axis=0); q5_1=sub1.quantile(0.05, axis=0); q25_1=sub1.quantile(0.25, axis=0); med1=sub1.quantile(0.5, axis=0); q75_1=sub1.quantile(0.75, axis=0); q95_1=sub1.quantile(0.95, axis=0); max1=sub1.max(axis=0)
ms1 = pd.DataFrame(np.round(mean1,2).astype(str) + " \xb1 " + np.round(std1,2).astype(str), columns=['Case_mean[std]'])
sss = pd.DataFrame([min1,q5_1,q25_1,med1,q75_1,q95_1,max1],index=['Case_min','Case_q5','Case_q25','Case_q50','Case_q75','Case_q95','Case_max']).astype(str).T
b = pd.merge(b,ms1,how='left',left_on='DF'+instance,right_index=True)
b = pd.merge(b,sss,how='left',left_on='DF'+instance,right_index=True)
temp = pd.merge(maps[['DF','description','data_type','Image_type','Sub_classification']], b, left_on='DF',right_on='DF', how='inner')
temp = temp.sort_values('P-value')
if isinstance(outcomeDescriptions, pd.DataFrame):
temp.insert(0,'Outcome_description',outcomeDescriptions[outcomeDescriptions['DF']==outcome]['description'].tolist()[0])
else:
temp.insert(0,'Outcome_description','')
temp.insert(0,'Outcome',outcome)
temp.to_csv(output_dir+'/sig2_'+LRtype+outfile_suffix+'__'+l+"_"+suffix+paraEst2+postfix+"__"+outcome+"_"+instance+'_.csv',sep=",",header=True,index=False)

function_module.py at commit 3ae4b25, no license · at the source

Overview

Authors: Satoshi Okawa^1,2,3, Yinan Shi¹, Somanshu Banerjee⁴, Yi Yin Tai¹, Jamal Nabhanizadeh⁵, Stefan Guenther⁶, Anisha N. Shah¹, Siyi Jiang^1,7, Timothy Kim⁸, Yassmin Al Aaraj¹, Caroline Pidro⁹, Metin Aytekin¹, Omolara Fatukasi¹⁰, Nishita Kalepalli¹, Joseph Mettenburg¹¹, Timothy D. Girard⁹, James T. Becker¹⁰, Kathiravan Kaliyappan^12,13, Shanshan Song^12,13, Sandra Breuils-Bonnet¹⁴

and 10 other authors

Steeve Provencher¹⁴, Sébastien Bonnet^5,14, Theodore Huppert¹⁵, Lloyd D. Harvey¹, Seyed Mehdi Nouraie¹⁶, Rajkumar Savai^5,6,14, Dandan Sun^12,13, Soban Umar⁴, Soni S. Pullamsetti^5,6,14, Stephen Y. Chan^1,8

ORCID iDs: Yi Yin Tai, Steeve Provencher, Lloyd D. Harvey, Rajkumar Savai, Stephen Y. Chan

16 affiliations

Center for Pulmonary Vascular Biology and Medicine, Pittsburgh Heart, Lung, and Blood Vascular Medicine Institute, Division of Cardiology, Department of Medicine
Department of Computational and Systems Biology; and
McGowan Institute for Regenerative Medicine, University of Pittsburgh School of Medicine, Pittsburgh, Pennsylvania, USA
Department of Anesthesiology and Perioperative Medicine, Division of Molecular Medicine, University of California, Los Angeles, Los Angeles, California, USA
Institute for Lung Health, Department of Internal Medicine, German Center for Lung Research (DZL), Cardio-pulmonary Institute (CPI), Justus Liebig University, Giessen, Germany
Max Planck Institute for Heart and Lung Research, DZL, CPI, Bad Nauheim, Germany
Cellular and Molecular Pathology Graduate Program, University of Pittsburgh, Pittsburgh, Pennsylvania, USA
Heart and Vascular Institute, University of Pittsburgh Medical Center (UPMC), Pittsburgh, Pennsylvania, USA
Department of Critical Care Medicine
Department of Psychiatry
Department of Radiology, and
Department of Neurology, University of Pittsburgh School of Medicine, Pittsburgh, Pennsylvania, USA
Pittsburgh Veteran’s Affairs Research Center, Pittsburgh, Pennsylvania, USA
PAH Research Group, Institut Universitaire de Cardiologie et Pneumologie de Québec Research Center, Department of Medicine, Université Laval, Quebec City, Québec, Canada
Department of Electrical and Computer Engineering, Swanson School of Engineering, University of Pittsburgh, Pittsburgh, Pennsylvania, USA
Division of Pulmonary, Allergy, and Critical Care Medicine, Department of Medicine, University of Pittsburgh School of Medicine, Pittsburgh, Pennsylvania, USA

Institutions: University of Pittsburgh (United States); McGowan Institute for Regenerative Medicine (United States); University of California, Los Angeles (United States); Justus-Liebig-Universität Gießen (Germany); German Center for Lung Research (Germany); Cardio-Pulmonary Institute (Germany); Max Planck Institute for Heart and Lung Research (Germany); UPMC Health System (United States); University of Pittsburgh Medical Center (United States); Institut universitaire de cardiologie et de pneumologie de Québec (Canada); Université Laval (Canada)

Journal: The Journal of clinical investigation, volume 136, issue 21, article e195110

Dates: received 9 May 2025; accepted 3 September 2026; published online 2 November 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1172/jci195110 · PMID 42969225 · PMCID PMC13722401 · OpenAlex W7225064535

Open access: gold, a free copy (OpenAlex)

Status: code verified (how the code is found and verified)

Categories: genetics / omics (modality), human (organism), mouse (organism), other condition (population)

Methods: Spectral & time-frequency, Statistics

Keywords: Inflammation, Neuroscience, Pulmonology, Genetic risk factors, Neuroimaging, Transcriptomics

MeSH: Cerebellum*, Frontal Lobe*, Hypertension, Pulmonary*, Lung*, Aged, Animals, Female, Humans, Male, Mice, Mice, Knockout, Middle Aged, White Matter (* major topic)

Topic: Pulmonary Hypertension Research and Treatments (Pulmonary and Respiratory Medicine, Medicine), according to OpenAlex

Funding: WoodNext Foundation (Research Grant); Veterans Administration Medical Center (IK6 BX005647-01,I01BX006894-01); National Institutes of Health (R01HL161038, R35GM164224); American Heart Association (24CDA1274310, 26SFRNPVHD1609644, 24SFRNCCN1276089,24SFRNPCN1280228,25ISFRNCP1501443); National Aeronautics and Space Administration (80NSSC26K0720); Deutsche Forschungsgemeinschaft (SA 1923/7-1, Collaborative Research Center 1213 (A01/A05/A10N) (268555672)); European Research Council (866051)

Citations: not cited yet (Europe PMC); 82 references in the paper

Abstract

Connections between vascular and neurocognitive diseases are emerging but remain poorly defined. Cognitive impairment is increasingly recognized in pulmonary hypertension (PH), a disease traditionally viewed as limited to the heart and lungs. Given that inflammation is shared by PH and neurodegenerative diseases, we investigated its role in driving brain and cognitive deficits in PH. Using structural and functional MRI data from UK Biobank (N = 248 PH) and a prospective validation cohort, we identified reduced cerebellar white matter volume, thinner frontal gyri, and weakened cerebellar-frontal connectivity in patients with PH but not in other heart or lung diseases. These changes were accompanied by worsened executive function and processing speed. Olink proteomics and targeted genomics linked inflammatory plasma proteins and a PH-associated nuclear receptor coactivator 7 (NCOA7) single nucleotide polymorphism, rs11154337, to these deficits. PH cerebellum showed increased microglial activation and white matter derangement. In PH mice, Ncoa7 deficiency worsened cerebellar inflammation and atrophy. Microglial functional assays, cerebellar single-nucleus RNA sequencing, and spatial transcriptomic analyses further implicated Ncoa7-dependent lysosomal, inflammatory, and cell-specific mechanisms of cerebellar pathobiology. These findings identify NCOA7-linked inflammatory dysfunction in the cerebellar-frontal axis in PH, redefining this disease with neurologic consequences and thus offering neural targets for diagnostic, prognostic, and therapeutic development.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 2 matches between paragraphs and lines of code.

PittVMIOSB/Okawa_et_al

License: none: the authors keep all their rights

State: the link answers, verified on 3 November 2026

Evidence: files inventoried

Commit: 3ae4b2540f1491eca8a5292851ea8d5ac437de4a, 22 April 2025

Languages: Python (4)

Size: 14 files, 4 scripts

Software Heritage: not archived

Found in: “Data availability.”

Holds: README

Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: NumPy (4 files), pandas (4 files), statsmodels (4 files), Matplotlib (1 file), scikit-learn (1 file), SciPy (1 file), seaborn (1 file)

Availability: 1 check, the latest on 3 November 2026: the link answers

3 November 2026: the link answers

5 files

function_module.py, Python, 1,235 lines, 2 matches
mouseMRI_analysis.py, Python, 140 lines
ukbb_MRI_analysis.py, Python, 227 lines
upmc_MRI_analysis.py, Python, 150 lines
README.md, Text, 1 line

At the source: github.com/PittVMIOSB/Okawa_et_al

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
4 scripts, each with its path and the digest of its content;
2 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

ukbiobank.ac.uk/enable-your-research, at UK Biobank; found in “Data availability.”

Data availability

All materials and data from this study are available to researchers upon reasonable request and after issuance of relevant institutional data sharing agreements by contacting the corresponding authors. Human data with identifiers are restricted owing to patient confidentiality. The UKBB data can be accessed through the UKBB Research Analysis Portal (https://www.ukbiobank.ac.uk/enable-your-research). Raw and processed snRNA-seq and spatial transcriptomics data have been deposited in Gene Expression Omnibus under accession numbers GSE334773 and GSE334959, respectively. Analysis code is available with open access at the following GitHub repository: https://github.com/PittVMIOSB/Okawa_et_al (commit ID 3ae4b25). Any additional information required to reanalyze the data reported in this paper is available from the corresponding authors upon request. All raw data values are included in the Supporting Data Values file.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 3 November 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 30 authors, 6 keywords, 13 MeSH terms, 7 funders, 82 references.

Cite

This paper

Okawa, S., Shi, Y., Banerjee, S., Tai, Y. Y., Nabhanizadeh, J., Guenther, S., Shah, A. N., Jiang, S., Kim, T., Al Aaraj, Y., Pidro, C., Aytekin, M., Fatukasi, O., Kalepalli, N., Mettenburg, J., Girard, T. D., Becker, J. T., Kaliyappan, K., Song, S., . . . Chan, S. Y. (2026). Cerebellar and frontal lobe alterations in pulmonary hypertension implicate an inflammatory lung vascular/brain axis. The Journal of clinical investigation, 136(21), e195110. https://doi.org/10.1172/jci195110

BibTeX

@article{okawa2026cerebellar, author = {Okawa, Satoshi and Shi, Yinan and Banerjee, Somanshu and Tai, Yi Yin and Nabhanizadeh, Jamal and Guenther, Stefan and Shah, Anisha N. and Jiang, Siyi and Kim, Timothy and Al Aaraj, Yassmin and Pidro, Caroline and Aytekin, Metin and Fatukasi, Omolara and Kalepalli, Nishita and Mettenburg, Joseph and Girard, Timothy D. and Becker, James T. and Kaliyappan, Kathiravan and Song, Shanshan and Breuils-Bonnet, Sandra and Provencher, Steeve and Bonnet, Sébastien and Huppert, Theodore and Harvey, Lloyd D. and Nouraie, Seyed Mehdi and Savai, Rajkumar and Sun, Dandan and Umar, Soban and Pullamsetti, Soni S. and Chan, Stephen Y.}, title = {{Cerebellar and frontal lobe alterations in pulmonary hypertension implicate an inflammatory lung vascular/brain axis}}, journal = {The Journal of clinical investigation}, year = {2026}, month = nov, volume = {136}, number = {21}, pages = {e195110}, publisher = {American Society for Clinical Investigation}, issn = {0021-9738}, doi = {10.1172/jci195110}, url = {https://doi.org/10.1172/jci195110}, pmid = {42969225}, pmcid = {PMC13722401} }

RIS

TY - JOUR AU - Okawa, Satoshi AU - Shi, Yinan AU - Banerjee, Somanshu AU - Tai, Yi Yin AU - Nabhanizadeh, Jamal AU - Guenther, Stefan AU - Shah, Anisha N. AU - Jiang, Siyi AU - Kim, Timothy AU - Al Aaraj, Yassmin AU - Pidro, Caroline AU - Aytekin, Metin AU - Fatukasi, Omolara AU - Kalepalli, Nishita AU - Mettenburg, Joseph AU - Girard, Timothy D. AU - Becker, James T. AU - Kaliyappan, Kathiravan AU - Song, Shanshan AU - Breuils-Bonnet, Sandra AU - Provencher, Steeve AU - Bonnet, Sébastien AU - Huppert, Theodore AU - Harvey, Lloyd D. AU - Nouraie, Seyed Mehdi AU - Savai, Rajkumar AU - Sun, Dandan AU - Umar, Soban AU - Pullamsetti, Soni S. AU - Chan, Stephen Y. TI - Cerebellar and frontal lobe alterations in pulmonary hypertension implicate an inflammatory lung vascular/brain axis T2 - The Journal of clinical investigation J2 - J Clin Invest PY - 2026 DA - 2026/11/02 VL - 136 IS - 21 SP - e195110 SN - 0021-9738 PB - American Society for Clinical Investigation DO - 10.1172/jci195110 UR - https://doi.org/10.1172/jci195110 LA - en ER -

CSL-JSON

{ "id": "10.1172/jci195110", "type": "article-journal", "title": "Cerebellar and frontal lobe alterations in pulmonary hypertension implicate an inflammatory lung vascular/brain axis", "container-title": "The Journal of clinical investigation", "author": [ { "family": "Okawa", "given": "Satoshi" }, { "family": "Shi", "given": "Yinan" }, { "family": "Banerjee", "given": "Somanshu" }, { "family": "Tai", "given": "Yi Yin" }, { "family": "Nabhanizadeh", "given": "Jamal" }, { "family": "Guenther", "given": "Stefan" }, { "family": "Shah", "given": "Anisha N." }, { "family": "Jiang", "given": "Siyi" }, { "family": "Kim", "given": "Timothy" }, { "family": "Al Aaraj", "given": "Yassmin" }, { "family": "Pidro", "given": "Caroline" }, { "family": "Aytekin", "given": "Metin" }, { "family": "Fatukasi", "given": "Omolara" }, { "family": "Kalepalli", "given": "Nishita" }, { "family": "Mettenburg", "given": "Joseph" }, { "family": "Girard", "given": "Timothy D." }, { "family": "Becker", "given": "James T." }, { "family": "Kaliyappan", "given": "Kathiravan" }, { "family": "Song", "given": "Shanshan" }, { "family": "Breuils-Bonnet", "given": "Sandra" }, { "family": "Provencher", "given": "Steeve" }, { "family": "Bonnet", "given": "Sébastien" }, { "family": "Huppert", "given": "Theodore" }, { "family": "Harvey", "given": "Lloyd D." }, { "family": "Nouraie", "given": "Seyed Mehdi" }, { "family": "Savai", "given": "Rajkumar" }, { "family": "Sun", "given": "Dandan" }, { "family": "Umar", "given": "Soban" }, { "family": "Pullamsetti", "given": "Soni S." }, { "family": "Chan", "given": "Stephen Y." } ], "container-title-short": "J Clin Invest", "volume": "136", "issue": "21", "page": "e195110", "DOI": "10.1172/jci195110", "PMID": "42969225", "PMCID": "PMC13722401", "ISSN": "0021-9738", "publisher": "American Society for Clinical Investigation", "URL": "https://doi.org/10.1172/jci195110", "language": "en", "issued": { "date-parts": [ [ 2026, 11, 2 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41591-024-03426-4 [code]: Plasma proteomic evidence for increased β-amyloid pathology after SARS-CoV-2 infection.
Journal: Nature medicine
In common: statsmodels, seaborn, scikit-learn, 4 other tools, ukbiobank.ac.uk/enable-your-research, genetics / omics, other condition
[2] doi:10.1002/hbm.26267 [code]: Optic radiations representing different eccentricities age differently.
Journal: Human brain mapping
In common: statsmodels, seaborn, pandas, 3 other tools, ukbiobank.ac.uk/enable-your-research
[3] doi:10.1038/s41586-026-10345-6 [code]: Population-scale repeat expansions elucidate disease risk and brain atrophy.
Journal: Nature
In common: statsmodels, pandas, SciPy, 2 other tools, ukbiobank.ac.uk/enable-your-research, genetics / omics, other condition
[4] doi:10.1038/s41467-026-72875-x [code]: Lifespan normative modeling of brain microstructure.
Journal: Nature communications
In common: seaborn, scikit-learn, pandas, 3 other tools, ukbiobank.ac.uk/enable-your-research
[5] doi:10.1038/s41398-025-03290-x [code]: Neurofind: using deep learning to make individualised inferences in brain-based disorders.
Journal: Translational psychiatry
In common: seaborn, scikit-learn, pandas, 3 other tools, ukbiobank.ac.uk/enable-your-research
[6] doi:10.1038/s43856-024-00496-w [code]: Convolutional neural network-based classification of glaucoma using optic radiation tissue properties.
Journal: Communications medicine
In common: statsmodels, scikit-learn, pandas, 2 other tools, ukbiobank.ac.uk/enable-your-research, other condition
[7] doi:10.1016/j.celrep.2025.116624 [code]: Spatial and single-cell transcriptomics reveal the reorganization of cerebellar microglia with aging.
Journal: Cell reports
In common: statsmodels, seaborn, scikit-learn, 4 other tools, genetics / omics, mouse, 2 references
[8] doi:10.1038/s41467-025-67578-8 [code]: Genetic Insights into Head-to-Body Ratios Via Deep Learning-Based Image Segmentation and Implications for Common Diseases.
Journal: Nature communications
In common: seaborn, scikit-learn, pandas, 2 other tools, ukbiobank.ac.uk/enable-your-research, genetics / omics
[9] doi:10.1038/s41514-026-00456-9 [code]: Exploring the link between body physiology and cognition: the role of the brain and aging.
Journal: npj aging
In common: statsmodels, seaborn, scikit-learn, 4 other tools, 2 references
[10] doi:10.3389/fnimg.2022.1031253 [code]: Automatic measure and normalization of spinal cord cross-sectional area using the pontomedullary junction.
Journal: Frontiers in neuroimaging
In common: seaborn, pandas, SciPy, 2 other tools, ukbiobank.ac.uk/enable-your-research

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Add the badge to the code's README

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown

<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.1172/jci195110.svg)](https://openscicode.org/badge/doi/10.1172/jci195110)
<!-- oscr-badge:end -->

reStructuredText

.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.1172/jci195110.svg
   :target: https://openscicode.org/badge/doi/10.1172/jci195110
   :alt: OSCR

.. oscr-badge:end

HTML

<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.1172/jci195110"><img src="https://openscicode.org/badge/doi/10.1172/jci195110.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

PittVMIOSB/Okawa_et_al: README.md on GitHub

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.