Running LDSC on raw sumstats

The below code runs ldsc() on a collection of traits used throughout the paper. However, we have only provided the summary statistics from the novel GWAS conducted in the course of this study. For summary statistics from GWAS conducted by other researchers, to the citations found in Tables S2, S7 and S8.

Raw summary statistics must be wrangled/munged prior to submitting them to ldsc(). Guidance on this process may be found here: https://github.com/GenomicSEM/GenomicSEM/wiki/2.1-Calculating-Sum-of-Effective-Sample-Size-and-Preparing-GWAS-Summary-Statistics

LD-reference files used below may be found here: https://github.com/GenomicSEM/GenomicSEM/wiki/3.-Genome%E2%80%90wide-Models

traits<-c("regenie_cut1.sumstats.gz",
          "regenie_cut2.sumstats.gz",
          "regenie_cut3.sumstats.gz",
          "regenie_cut4.sumstats.gz",
          "regenie_cut5.sumstats.gz",
          "regenie_cut6.sumstats.gz",
          "regenie_cut7.sumstats.gz",
          "regenie_cut8.sumstats.gz",
          "regenie_cut9.sumstats.gz",
          "regenie_cut10.sumstats.gz",
          "regenie_cut11.sumstats.gz",
          "regenie_cut12.sumstats.gz",
          "MDD.sumstats.gz",
          "Wellbeing.sumstats.gz",
          "Neuroticism.sumstats.gz",
          "cut1_diffactor.sumstats.gz",
          "cut11_diffactor.sumstats.gz",
          "PTSD.sumstats.gz",
          "SCZ.2n.sumstats.gz",
          "ANX.sumstats.gz",
          "BIP.2n.sumstats.gz",
          "ADHD.sumstats.gz",
          "ASD.sumstats.gz",
          "OCD.sumstats.gz",
          "AN.sumstats.gz",
          "COMP.sumstats.gz",
          "INT.sumstats.gz",
          "NEURODEV.sumstats.gz",
          "THOUGHT.sumstats.gz",
          "ADDICT.sumstats.gz",
          "OUD.sumstats.gz",
          "ALCH.sumstats.gz",
          "TS.sumstats.gz",
          "CUD.sumstats.gz"
          )
 
#entering sample prevalences
sample.prev<-c(0.85556773,
               0.75867410,
               0.64767354,
               0.53828939,
               0.43558773,
               0.34194621,
               0.25818629,
               0.19038580,
               0.13159727,
               0.08571119,
               0.04947786,
               0.02166147,
               0.5,
               NA,
               NA,
               0.5,
               0.5,
               0.5,
               0.5,
               0.3,
               0.5,
               0.5,
               0.5,
               0.5,
               0.5,
               0.5,
               0.5,
               0.5,
               0.5,
               0.5,
               0.22,
               0.5, 
               0.5,
               0.5
               )

###vector of population prevalences (it is assumed the sample prevs match population
###prevalences for each neuroticism cut point)
population.prev<-c(0.85556773,
                 0.75867410,
                 0.64767354,
                 0.53828939,
                 0.43558773,
                 0.34194621,
                 0.25818629,
                 0.19038580,
                 0.13159727,
                 0.08571119,
                 0.04947786,
                 0.02166147, 
                 0.15, 
                 NA,
                 NA,
                 0.85556773, 
                 0.04947786,
                 0.15,
                 0.01,
                 0.20,
                 0.02,
                 0.06,
                 0.01,
                 0.015,
                 0.01,
                 NA,
                 NA,
                 NA,
                 NA,
                 NA,
                 0.01,
                 0.01,
                 0.006,
                 0.02
                  )

###the folder of LD scores
  ld<-"eur_w_ld_chr/"

###the folder of LD weights [typically the same as folder of LD scores]
  wld<-"eur_w_ld_chr/"

###name the traits
trait.names<-c("cut1",
               "cut2",
               "cut3",
               "cut4",
               "cut5",
               "cut6",
               "cut7",
               "cut8",
               "cut9",
               "cut10",
               "cut11",
               "cut12",
               "MDD",
               "Wellbeing",
               "Neuroticism",
               "cut1_diffactor",
               "cut11_diffactor",
               "ptsd",
               "scz",
               "anxiety",
               "bipolar",
               "adhd",
               "autism",
               "ocd",
               "an",
               "compulsive",
               "internalizing",
               "neurodev",
               "thought",
               "addict",
               "oud",
               "alch",
               "ts",
               "cud"
               )

###run LDSC
  LDSCoutput<-ldsc(traits=traits,
                   sample.prev=sample.prev,
                   population.prev=population.prev,
                   ld=ld,
                   wld=wld,
                   trait.names=trait.names,
                   stand = TRUE)
                 #n.blocks = 1000)

summary(LDSCoutput)

GTACCC of Negative Emotionality in Relation to Major Depressive Disorder

The below code reproduces Figure 1 from the paper, and outputs the parameter estimates of the linear GLS regression of MDD (REGENIE) found in Extened Data Table 2.

#Loading LDSC object 
load("LDSCoutput_Regenie.RData")

###Predictors for GLS
PREDICTORS <- -1*qnorm(c(0.85556773,0.75867410,0.64767354,
                         0.53828939,0.43558773,0.34194621,
                         0.25818629,0.19038580,0.13159727, 
                         0.08571119,0.04947786,0.02166147))

summaryGLSbands(OBJECT = subSV(LDSC_OBJECT = LDSCoutput,
                INDEXVALS = indexS(MATRIX = LDSCoutput$S, R = F)[1:12,13], TYPE = "S_Stand"),
                PREDICTORS = PREDICTORS,
                INTERVALS = 40
                ,
                XLAB = "Neuroticism Severity (Z-Distribution)",
                YLAB = "Genetic Correlation w/ MDD",
                QUAD = F,
                XCOORDS = c(-1,2),
                YCOORDS = c(0,1))
##         betas        pvals         SE        Z
## b0 0.52648636 1.055411e-18 0.05963095 8.829079
## b1 0.09049874 1.584398e-03 0.02865005 3.158764

For parameter estimates from the quadratic GLS set the argument QUAD = T. Parameter estimates for GLS regression of the remaining phenotypes in Extended Data Table 2, and Tables S4-S6,S9 may be found by specifying the columns associated with those phenotypes with the appropriate LDSC object.

#Loading LDSC object 
load("LDSCoutput_Regenie.RData")

###Predictors for GLS
PREDICTORS <- -1*qnorm(c(0.85556773,0.75867410,0.64767354,
                         0.53828939,0.43558773,0.34194621,
                         0.25818629,0.19038580,0.13159727, 
                         0.08571119,0.04947786,0.02166147))

summaryGLSbands(OBJECT = subSV(LDSC_OBJECT = LDSCoutput,
                INDEXVALS = indexS(MATRIX = LDSCoutput$S, R = F)[1:12,13], TYPE = "S_Stand"),
                PREDICTORS = PREDICTORS,
                INTERVALS = 40
                ,
                XLAB = "Neuroticism Severity (Z-Distribution)",
                YLAB = "Genetic Correlation w/ MDD",
                QUAD = T,
                XCOORDS = c(-1,2),
                YCOORDS = c(0,1))
##         betas        pvals         SE        Z
## b0 0.50936285 4.495604e-17 0.06064427 8.399192
## b1 0.08047254 6.145083e-03 0.02937021 2.739938
## b2 0.03375832 1.208844e-01 0.02176451 1.551072

Parameter estimates for GLS regression of the remaining phenotypes in Extended Data Table 2, and Tables S4-S6 may be found by specifying the columns associated with those prototypes.

colnames(LDSCoutput$S_Stand)
##  [1] "cut1"        "cut2"        "cut3"        "cut4"        "cut5"       
##  [6] "cut6"        "cut7"        "cut8"        "cut9"        "cut10"      
## [11] "cut11"       "cut12"       "MDD"         "Wellbeing"   "Neuroticism"

Modelling of Severity Factors using Genomic SEM

Below is code to reproduce the models presented in the main text in a section titled “Modelling of Severity Factors using Genomic SEM,” as well as in Figure 2.

One-Factor Model

load("LDSCoutput_Regenie.RData")

one_severity_factors_model <-'Severity =~ NA*cut1 + cut2 + cut3 + cut4 + 
                                      cut5 + cut6 + cut7 + cut8 + cut9 + cut10 + cut11 + cut12
                              MDD ~~ Severity
                              Severity ~~ 1*Severity'

###Fitting the model
one_severity_factor <-usermodel(LDSCoutput,
                                     model = one_severity_factors_model)

One-factor Model Fit Indices

chisq df p_chisq AIC CFI SRMR
1604.257 65 0 1656.257 0.999 0.053

One-factor Model Parameter Estimates

lhs op rhs Unstand_Est Unstand_SE STD_Genotype STD_Genotype_SE STD_All p_value
14 Severity =~ cut1 0.297 0.00838889668556512 0.933 0.026305011517975 0.933 1.98624135526816e-275
18 Severity =~ cut2 0.329 0.00779549585042523 0.956 0.022682421349244 0.956 < 5e-300
19 Severity =~ cut3 0.346 0.00729731521832829 0.980 0.0206650261455009 0.980 < 5e-300
20 Severity =~ cut4 0.351 0.00775926135178155 0.991 0.0219097866159061 0.991 < 5e-300
21 Severity =~ cut5 0.354 0.00741269693815855 0.993 0.0208112800777277 0.993 < 5e-300
22 Severity =~ cut6 0.359 0.00730251254186436 0.996 0.0202704526869198 0.996 < 5e-300
23 Severity =~ cut7 0.364 0.00780566478218322 0.995 0.0213088434336923 0.995 < 5e-300
24 Severity =~ cut8 0.361 0.0078864381630413 0.995 0.0217202107769853 0.995 < 5e-300
25 Severity =~ cut9 0.351 0.00851795202167433 0.974 0.0236543344182371 0.974 < 5e-300
15 Severity =~ cut10 0.337 0.00913890651220113 0.939 0.0254572160668675 0.939 5.08511326967504e-298
16 Severity =~ cut11 0.312 0.00997702465852634 0.926 0.0296078010200161 0.926 7.31038860156998e-215
17 Severity =~ cut12 0.291 0.0119559534445321 0.935 0.0383910669233824 0.935 4.56364956060229e-131
26 Severity ~~ MDD 0.176 0.0189751263618886 0.594 0.0639780055994698 0.594 1.74530148103812e-20
1 cut1 ~~ cut1 0.013 0.00311301188212919 0.130 0.0306089257939989 0.130 2.19822019510362e-05
5 cut2 ~~ cut2 0.010 0.00195588179665866 0.086 0.0165590001129238 0.086 1.84820762378177e-07
6 cut3 ~~ cut3 0.005 0.00150641383111434 0.039 0.0120806410523472 0.039 0.00109677221042213
7 cut4 ~~ cut4 0.002 0.001123558289508 0.018 0.00895841259712292 0.018 0.050264853416805
8 cut5 ~~ cut5 0.002 0.000995104215645592 0.014 0.00784356832534762 0.014 0.0689034499678154
9 cut6 ~~ cut6 0.001 0.0010695813370285 0.007 0.00824130702538061 0.007 0.376356562109321
10 cut7 ~~ cut7 0.001 0.0011857300146954 0.010 0.00883661208799619 0.010 0.251747360272973
11 cut8 ~~ cut8 0.001 0.00130942799743943 0.011 0.00993225410088934 0.011 0.285252166066245
12 cut9 ~~ cut9 0.007 0.001893610320315 0.051 0.0146029877239283 0.051 0.000453526983030212
2 cut10 ~~ cut10 0.015 0.00287755642848173 0.118 0.0223283826951658 0.118 1.35269527703377e-07
3 cut11 ~~ cut11 0.016 0.00449643867453152 0.142 0.0395985196458468 0.142 0.000336825227366282
4 cut12 ~~ cut12 0.012 0.00850015904070892 0.125 0.0876433532457637 0.125 0.152478422458314
13 MDD ~~ MDD 0.088 0.0118036802126993 1.000 0.13418670710626 1.000 9.17255016432065e-14
27 Severity ~~ Severity 1.000 1.000 1.000 NA

Two-factor Model

two_severity_factors_model <-'Flow =~ NA*cut1 + cut2 + cut3 + cut4 + 
                                      cut5 + cut6 + cut7 + cut8 + cut9 + cut10 + cut11
                              Fhigh =~ NA*cut2 + cut3 + cut4 + cut5 +
                                       cut6 + cut7 + cut8 + cut9 + cut10 + cut11 + cut12
                              Flow ~~ Fhigh
                              MDD ~~ Flow + Fhigh
                            
                              Flow ~~ 1*Flow
                              Fhigh ~~ 1*Fhigh'

###Fitting the model
two_severity_factors <-usermodel(LDSCoutput,
                                     model = two_severity_factors_model)

Two-factor Model Fit Indices

chisq df p_chisq AIC CFI SRMR
194.173 53 0 270.173 1 0.012

Two-factor Model Parameter Estimates

lhs op rhs Unstand_Est Unstand_SE STD_Genotype STD_Genotype_SE STD_All p_value
26 Flow =~ cut1 0.312 0.00911560118289078 0.979 0.0285837339428487 0.979 0.000
29 Flow =~ cut2 0.311 0.0185001134993025 0.906 0.0538294790206733 0.906 0.000
30 Flow =~ cut3 0.299 0.0225116779934765 0.847 0.0637500869741013 0.847 0.000
31 Flow =~ cut4 0.277 0.0211657043903759 0.782 0.0597654943534013 0.782 0.000
32 Flow =~ cut5 0.266 0.0237175145036054 0.746 0.0665873580533125 0.746 0.000
33 Flow =~ cut6 0.238 0.0248905996647439 0.662 0.0690918080643291 0.662 0.000
34 Flow =~ cut7 0.209 0.0269539940375222 0.571 0.0735822644811752 0.571 0.000
35 Flow =~ cut8 0.159 0.0293975431059885 0.437 0.0809644066794419 0.437 0.000
36 Flow =~ cut9 0.097 0.0342365224355613 0.269 0.0950747574225929 0.269 0.005
27 Flow =~ cut10 0.029 0.0405794161094462 0.080 0.113037540999848 0.080 0.479
28 Flow =~ cut11 -0.001 0.0355603063537817 -0.004 0.105528704800155 -0.004 0.973
16 Fhigh =~ cut2 0.033 0.0198183341598964 0.097 0.0576650855169325 0.097 0.094
17 Fhigh =~ cut3 0.064 0.024322978071973 0.180 0.0688794529229774 0.180 0.009
18 Fhigh =~ cut4 0.090 0.0226199481495772 0.253 0.0638718372212661 0.253 0.000
19 Fhigh =~ cut5 0.103 0.0249751528013201 0.291 0.0701182011895148 0.291 0.000
20 Fhigh =~ cut6 0.136 0.0254211932738472 0.376 0.0705646410099901 0.376 0.000
21 Fhigh =~ cut7 0.171 0.0270822946223777 0.467 0.0739325216481835 0.467 0.000
22 Fhigh =~ cut8 0.220 0.027465887107612 0.605 0.0756443970894861 0.605 0.000
23 Fhigh =~ cut9 0.273 0.0311776757839206 0.759 0.0865803500528233 0.759 0.000
13 Fhigh =~ cut10 0.331 0.0369675703627583 0.923 0.10297642613481 0.923 0.000
14 Fhigh =~ cut11 0.336 0.033121920834356 0.997 0.0982925626605277 0.997 0.000
15 Fhigh =~ cut12 0.312 0.0141425373608029 1.001 0.0454122801649043 1.001 0.000
37 Flow ~~ Fhigh 0.821 0.0431464462961825 0.821 0.0431464457040539 0.821 0.000
39 Flow ~~ MDD 0.133 0.0201611806788112 0.449 0.0679769988499139 0.449 0.000
25 Fhigh ~~ MDD 0.214 0.0249209011514986 0.721 0.0840252393687385 0.721 0.000
1 cut1 ~~ cut1 0.004 0.00249206280859702 0.042 0.024503395286935 0.042 0.087
5 cut2 ~~ cut2 0.003 0.0013818411660784 0.025 0.0116990220690767 0.025 0.030
6 cut3 ~~ cut3 0.000 0.000862365861376442 0.001 0.00691571765404036 0.001 0.896
7 cut4 ~~ cut4 0.000 0.000738657543684859 -0.001 0.00588950210655267 -0.001 0.825
8 cut5 ~~ cut5 0.000 0.000921896435972384 0.003 0.00726653229127145 0.003 0.636
9 cut6 ~~ cut6 0.001 0.000982954635951843 0.011 0.00757383328347134 0.011 0.150
10 cut7 ~~ cut7 0.002 0.00117304756070602 0.019 0.00874209639911103 0.019 0.033
11 cut8 ~~ cut8 0.001 0.00126513164449796 0.010 0.00959625773415476 0.010 0.275
12 cut9 ~~ cut9 0.002 0.00151332598243703 0.016 0.0116703411018504 0.016 0.161
2 cut10 ~~ cut10 0.003 0.00235667072683269 0.020 0.0182865821028072 0.020 0.280
3 cut11 ~~ cut11 0.001 0.0030114485432539 0.011 0.0265207386376074 0.011 0.681
4 cut12 ~~ cut12 0.000 0.00690451296355195 -0.003 0.0711909808174056 -0.003 0.971
40 MDD ~~ MDD 0.088 0.0118036802126993 1.000 0.13418670710626 1.000 0.000
38 Flow ~~ Flow 1.000 1.000 1.000 NA
24 Fhigh ~~ Fhigh 1.000 1.000 1.000 NA

Three-factor Model

three_severity_factors_model <- '
  Flow =~ NA*cut1 + cut2 + cut3 + cut4 + cut5 
  Fmid =~ NA*cut2 + cut3 + cut4 + cut5 + cut6 + cut7 + cut8 + cut9 + cut10 + cut11  
  Fhigh =~ NA*cut8 + cut9 + cut10 + cut11 + cut12

  Flow ~~ Fhigh
  Flow ~~ Fmid
  Fmid ~~ Fhigh

  MDD ~~ Flow + Fhigh + Fmid 

  Flow ~~ 1*Flow
  Fhigh ~~ 1*Fhigh
  Fmid ~~ 1*Fmid
'

three_severity_factors_fit <- usermodel(
  LDSCoutput,
  model = three_severity_factors_model
)

Three-factor Model Fit Indices

chisq df p_chisq AIC CFI SRMR
182.112 52 0 260.112 1 0.012

Three-factor Model Parameter Estimates

lhs op rhs Unstand_Est Unstand_SE STD_Genotype STD_Genotype_SE STD_All p_value
20 Flow =~ cut1 0.315 0.00917486072887723 0.987 0.0287695546469164 0.987 6.65101180598834e-258
21 Flow =~ cut2 0.264 0.0434484779674912 0.768 0.126421361324284 0.768 1.25688357024297e-09
22 Flow =~ cut3 0.198 0.0442819873179708 0.561 0.125400811817729 0.561 7.68521868199764e-06
23 Flow =~ cut4 0.143 0.0355950647028522 0.404 0.100509696669515 0.404 5.7153275738797e-05
24 Flow =~ cut5 0.110 0.0312395782313007 0.309 0.0877056974948271 0.309 0.00041760667366969
31 Fmid =~ cut2 0.079 0.0435027537857497 0.231 0.12657928096463 0.231 0.0677877806525361
32 Fmid =~ cut3 0.159 0.0437573048920025 0.450 0.123914981530592 0.450 0.000282256313020956
33 Fmid =~ cut4 0.215 0.0350669671866059 0.608 0.09901851438389 0.608 8.40958255966121e-10
34 Fmid =~ cut5 0.249 0.0312302291834818 0.699 0.0876794434499092 0.699 1.54647218831524e-15
35 Fmid =~ cut6 0.359 0.00731054013371482 0.997 0.0202927353432421 0.997 < 5e-300
36 Fmid =~ cut7 0.365 0.00772519097444024 0.995 0.0210891553264821 0.995 < 5e-300
37 Fmid =~ cut8 0.258 0.0299684988450407 0.710 0.0825368655146348 0.710 8.12158194464441e-18
38 Fmid =~ cut9 0.160 0.0439935663659197 0.445 0.12217010194426 0.445 0.000272511632911141
29 Fmid =~ cut10 0.050 0.0622745925128981 0.138 0.173471162853926 0.138 0.426143607208316
30 Fmid =~ cut11 -0.004 0.0583496642681591 -0.012 0.17315809606216 -0.012 0.944887916006604
16 Fhigh =~ cut8 0.110 0.0293587547169786 0.302 0.0808575518956057 0.302 0.000188659220245946
17 Fhigh =~ cut9 0.203 0.0434075376943991 0.565 0.120542686414728 0.565 2.76660591365678e-06
13 Fhigh =~ cut10 0.309 0.0611134512482911 0.861 0.170236699091965 0.861 4.27057406505603e-07
14 Fhigh =~ cut11 0.340 0.057369081646909 1.009 0.170248136247861 1.009 3.11726768423268e-09
15 Fhigh =~ cut12 0.313 0.0142723522757193 1.003 0.0458291222944534 1.003 2.83786088228149e-106
25 Flow ~~ Fhigh 0.815 0.0430770822658911 0.815 0.0430770803493763 0.815 6.40010146250752e-80
27 Flow ~~ Fmid 0.952 0.0206568749937307 0.952 0.02065688154784 0.952 < 5e-300
39 Fmid ~~ Fhigh 0.926 0.0259819355667933 0.926 0.0259819419860038 0.926 2.58909938202043e-278
28 Flow ~~ MDD 0.134 0.0216660456741549 0.450 0.0730509249342471 0.450 7.12661257549846e-10
19 Fhigh ~~ MDD 0.215 0.0258500120896669 0.724 0.0871578938296937 0.724 9.9127034249167e-17
41 Fmid ~~ MDD 0.172 0.0193041619742415 0.579 0.0650874068770236 0.579 5.76402231902629e-19
1 cut1 ~~ cut1 0.003 0.00250587694295073 0.026 0.0246392294956337 0.026 0.292051345272454
5 cut2 ~~ cut2 0.002 0.00123492348146257 0.019 0.0104551757086479 0.019 0.0661308329305538
6 cut3 ~~ cut3 0.000 0.000819955692304041 0.002 0.00657561011885159 0.002 0.735980991473121
7 cut4 ~~ cut4 0.000 0.000724170050003114 -0.001 0.00577398914991679 -0.001 0.88425097276853
8 cut5 ~~ cut5 0.000 0.00093406077656999 0.004 0.00736241058107339 0.004 0.627667081893922
9 cut6 ~~ cut6 0.001 0.00095873287938332 0.006 0.00738719953750142 0.006 0.400460770264718
10 cut7 ~~ cut7 0.001 0.000954845787567866 0.010 0.00711595645016197 0.010 0.180657263365667
11 cut8 ~~ cut8 0.001 0.00120879674089395 0.008 0.00916894837158114 0.008 0.35897074443956
12 cut9 ~~ cut9 0.002 0.00154187022931137 0.018 0.0118904724419645 0.018 0.141013573545347
2 cut10 ~~ cut10 0.003 0.00240523631614619 0.020 0.0186634153853504 0.020 0.287709449249249
3 cut11 ~~ cut11 0.001 0.00303882225446567 0.005 0.0267617984509638 0.005 0.863482109736404
4 cut12 ~~ cut12 -0.001 0.00671513070663947 -0.007 0.0692383030540618 -0.007 0.919945970921058
42 MDD ~~ MDD 0.088 0.0118036802126993 1.000 0.13418670710626 1.000 9.17244429758092e-14
26 Flow ~~ Flow 1.000 1.000 1.000 NA
18 Fhigh ~~ Fhigh 1.000 1.000 1.000 NA
40 Fmid ~~ Fmid 1.000 1.000 1.000 NA

GTACCC of Negative Emotionality in relation to Continuous Neuroticism and Subjective Well-being

Produces Extended Data Figure 3 and Figure S3

load("LDSCoutput_Regenie.RData")

PREDICTORS <- -1*qnorm(c(0.85556773,0.75867410,0.64767354,
                         0.53828939,0.43558773,0.34194621,
                         0.25818629,0.19038580,0.13159727, 
                         0.08571119,0.04947786,0.02166147))


###GTACCC of Continuous Measure of Neuroticism and neuroticism cut points
summaryGLSbands(OBJECT = subSV(LDSC_OBJECT = LDSCoutput,
                INDEXVALS = indexS(MATRIX = LDSCoutput$S, R = F)[1:12,15], TYPE = "S_Stand"),
                PREDICTORS = PREDICTORS,
                INTERVALS = 40
                ,
                XLAB = "Neuroticism Severity (Z-Distribution)",
                YLAB = "Genetic Correlation w/ Continuous Neuroticism",
                QUAD = F,
                XCOORDS = c(-1,2),
                YCOORDS = c(0,1))
##           betas         pvals         SE          Z
## b0  0.973659311 2.524544e-139 0.03874959 25.1269559
## b1 -0.009179323  3.646777e-01 0.01012625 -0.9064879

Three-factor Model (Neuroticism)

The models produced below are in large part indisinguishable from the above specified models - the only difference being a continuous measure of neuroticism isthe external trait the severity-factors correlate with

three_severity_factors_model <- '
  Flow =~ NA*cut1 + cut2 + cut3 + cut4 + cut5 
  Fmid =~ NA*cut2 + cut3 + cut4 + cut5 + cut6 + cut7 + cut8 + cut9 + cut10 + cut11  
  Fhigh =~ NA*cut8 + cut9 + cut10 + cut11 + cut12

  Flow ~~ Fhigh
  Flow ~~ Fmid
  Fmid ~~ Fhigh

  MDD ~~ Flow + Fhigh + Fmid 

  Flow ~~ 1*Flow
  Fhigh ~~ 1*Fhigh
  Fmid ~~ 1*Fmid
'

three_severity_factors_fit <- usermodel(
  LDSCoutput,
  model = three_severity_factors_model
)

Three-Factor Model Fit Indices (Neuroticism)

chisq df p_chisq AIC CFI SRMR
182.112 52 0 260.112 1 0.012

Three-Factor Model Parameter Estimates (Neuroticism)

lhs op rhs Unstand_Est Unstand_SE STD_Genotype STD_Genotype_SE STD_All p_value
20 Flow =~ cut1 0.315 0.00917486072887723 0.987 0.0287695546469164 0.987 6.65101180598834e-258
21 Flow =~ cut2 0.264 0.0434484779674912 0.768 0.126421361324284 0.768 1.25688357024297e-09
22 Flow =~ cut3 0.198 0.0442819873179708 0.561 0.125400811817729 0.561 7.68521868199764e-06
23 Flow =~ cut4 0.143 0.0355950647028522 0.404 0.100509696669515 0.404 5.7153275738797e-05
24 Flow =~ cut5 0.110 0.0312395782313007 0.309 0.0877056974948271 0.309 0.00041760667366969
31 Fmid =~ cut2 0.079 0.0435027537857497 0.231 0.12657928096463 0.231 0.0677877806525361
32 Fmid =~ cut3 0.159 0.0437573048920025 0.450 0.123914981530592 0.450 0.000282256313020956
33 Fmid =~ cut4 0.215 0.0350669671866059 0.608 0.09901851438389 0.608 8.40958255966121e-10
34 Fmid =~ cut5 0.249 0.0312302291834818 0.699 0.0876794434499092 0.699 1.54647218831524e-15
35 Fmid =~ cut6 0.359 0.00731054013371482 0.997 0.0202927353432421 0.997 < 5e-300
36 Fmid =~ cut7 0.365 0.00772519097444024 0.995 0.0210891553264821 0.995 < 5e-300
37 Fmid =~ cut8 0.258 0.0299684988450407 0.710 0.0825368655146348 0.710 8.12158194464441e-18
38 Fmid =~ cut9 0.160 0.0439935663659197 0.445 0.12217010194426 0.445 0.000272511632911141
29 Fmid =~ cut10 0.050 0.0622745925128981 0.138 0.173471162853926 0.138 0.426143607208316
30 Fmid =~ cut11 -0.004 0.0583496642681591 -0.012 0.17315809606216 -0.012 0.944887916006604
16 Fhigh =~ cut8 0.110 0.0293587547169786 0.302 0.0808575518956057 0.302 0.000188659220245946
17 Fhigh =~ cut9 0.203 0.0434075376943991 0.565 0.120542686414728 0.565 2.76660591365678e-06
13 Fhigh =~ cut10 0.309 0.0611134512482911 0.861 0.170236699091965 0.861 4.27057406505603e-07
14 Fhigh =~ cut11 0.340 0.057369081646909 1.009 0.170248136247861 1.009 3.11726768423268e-09
15 Fhigh =~ cut12 0.313 0.0142723522757193 1.003 0.0458291222944534 1.003 2.83786088228149e-106
25 Flow ~~ Fhigh 0.815 0.0430770822658911 0.815 0.0430770803493763 0.815 6.40010146250752e-80
27 Flow ~~ Fmid 0.952 0.0206568749937307 0.952 0.02065688154784 0.952 < 5e-300
39 Fmid ~~ Fhigh 0.926 0.0259819355667933 0.926 0.0259819419860038 0.926 2.58909938202043e-278
28 Flow ~~ MDD 0.134 0.0216660456741549 0.450 0.0730509249342471 0.450 7.12661257549846e-10
19 Fhigh ~~ MDD 0.215 0.0258500120896669 0.724 0.0871578938296937 0.724 9.9127034249167e-17
41 Fmid ~~ MDD 0.172 0.0193041619742415 0.579 0.0650874068770236 0.579 5.76402231902629e-19
1 cut1 ~~ cut1 0.003 0.00250587694295073 0.026 0.0246392294956337 0.026 0.292051345272454
5 cut2 ~~ cut2 0.002 0.00123492348146257 0.019 0.0104551757086479 0.019 0.0661308329305538
6 cut3 ~~ cut3 0.000 0.000819955692304041 0.002 0.00657561011885159 0.002 0.735980991473121
7 cut4 ~~ cut4 0.000 0.000724170050003114 -0.001 0.00577398914991679 -0.001 0.88425097276853
8 cut5 ~~ cut5 0.000 0.00093406077656999 0.004 0.00736241058107339 0.004 0.627667081893922
9 cut6 ~~ cut6 0.001 0.00095873287938332 0.006 0.00738719953750142 0.006 0.400460770264718
10 cut7 ~~ cut7 0.001 0.000954845787567866 0.010 0.00711595645016197 0.010 0.180657263365667
11 cut8 ~~ cut8 0.001 0.00120879674089395 0.008 0.00916894837158114 0.008 0.35897074443956
12 cut9 ~~ cut9 0.002 0.00154187022931137 0.018 0.0118904724419645 0.018 0.141013573545347
2 cut10 ~~ cut10 0.003 0.00240523631614619 0.020 0.0186634153853504 0.020 0.287709449249249
3 cut11 ~~ cut11 0.001 0.00303882225446567 0.005 0.0267617984509638 0.005 0.863482109736404
4 cut12 ~~ cut12 -0.001 0.00671513070663947 -0.007 0.0692383030540618 -0.007 0.919945970921058
42 MDD ~~ MDD 0.088 0.0118036802126993 1.000 0.13418670710626 1.000 9.17244429758092e-14
26 Flow ~~ Flow 1.000 1.000 1.000 NA
18 Fhigh ~~ Fhigh 1.000 1.000 1.000 NA
40 Fmid ~~ Fmid 1.000 1.000 1.000 NA

GTACCC of Well-Being

###GTACCC of SWB and neuroticism cut points
summaryGLSbands(OBJECT = subSV(LDSC_OBJECT = LDSCoutput,
                INDEXVALS = indexS(MATRIX = LDSCoutput$S, R = F)[1:12,14], TYPE = "S_Stand"),
                PREDICTORS = PREDICTORS,
                INTERVALS = 40
                ,
                XLAB = "Neuroticism Severity (Z-Distribution)",
                YLAB = "Genetic Correlation w/ Wellbeing",
                QUAD = F,
                XCOORDS = c(-1,2),
                YCOORDS = c(-1,0))
##          betas        pvals         SE         Z
## b0 -0.65396851 3.907785e-51 0.04347675 -15.04180
## b1 -0.02114249 1.762763e-01 0.01563436  -1.35231

IRT item-difficulty estimation

Item-level Genetic Correlations with MDD Across Item Severities Plots produced by this chunk may be found in figure S4

IRT_Output <- fread("IRT_Output")
load("LDSCoutput_itemlevel.RData")
#Predictors for endorsement-rate based difficulty estimation
predictors_qnorm <- IRT_Output$qnorm_value

## --- Plot 1: qnorm_value → rgwMDD (GLS, with error bands) ---
p_qnorm <- summaryGLSbands(
  OBJECT = subSV(
    LDSC_OBJECT = LDSCoutput,
    INDEXVALS   = indexS(MATRIX = LDSCoutput$S, R = FALSE)[1:12, 13],
    TYPE        = "S_Stand"
  ),
  PREDICTORS = predictors_qnorm,
  INTERVALS  = 30,
  XLAB       = "Item Severity (Endorsement Rate-Based Index)",
  YLAB       = "Genetic Correlation w/ \n MDD",
  QUAD       = FALSE,
  BANDS      = TRUE,
  BAND_SIZE  = 1,
  XCOORDS    = c(-0.2, 1),
  YCOORDS    = c(0, 1)
)
##         betas        pvals         SE        Z
## b0 0.33425210 5.988440e-13 0.04641910 7.200745
## b1 0.09593575 1.127488e-02 0.03785849 2.534062
p_qnorm <- p_qnorm +
  geom_text_repel(
    data = IRT_Output,
    aes(x = qnorm(1-IRT_Output$prev), y = rgwMDD, label = item),
    size = 5, segment.size = 0.2, min.segment.length = 0,
    inherit.aes = FALSE
  )

print(p_qnorm)

## --- Plot 2: difficulty → rgwMDD (GLS, with error bands) ---
#Predictors for IRT based difficulty estimation
predictors_difficulty <- IRT_Output$difficulty

p_diff <- summaryGLSbands(
  OBJECT = subSV(
    LDSC_OBJECT = LDSCoutput,
    INDEXVALS   = indexS(MATRIX = LDSCoutput$S, R = FALSE)[1:12, 13],
    TYPE        = "S_Stand"
  ),
  PREDICTORS = predictors_difficulty,
  INTERVALS  = 30,
  XLAB       = "Item Severity (IRT-Based Index)",
  YLAB       = "Genetic Correlation w/ \n MDD",
  QUAD       = F,
  BANDS      = TRUE,
  BAND_SIZE  = 1,
  XCOORDS    = c(-0.2, 1.6),
  YCOORDS    = c(0, 1)
)
##        betas        pvals         SE        Z
## b0 0.3365653 7.079926e-13 0.04688922 7.177883
## b1 0.0535362 4.093507e-02 0.02618948 2.044187
p_diff <- p_diff +
  geom_text_repel(
    data = IRT_Output,
    aes(x = difficulty, y = rgwMDD, label = item),
    size = 5, segment.size = 0.2, min.segment.length = 0,
    inherit.aes = FALSE
  )


print(p_diff)

GTACCC of Negative Emotionality in Relation to 13 Psychiatric Disorders and Their Transdiagnostic Factors

Produces plots in Figure 3, and figures S5-S9

load("LDSCoutput_Regenie_withALLdisorders.RData")

###Predictors for GLS
PREDICTORS <- -1*qnorm(c(0.85556773,0.75867410,0.64767354,
                         0.53828939,0.43558773,0.34194621,
                         0.25818629,0.19038580,0.13159727, 
                         0.08571119,0.04947786,0.02166147))

##Here the user can identify the trait desired and it's index number
##needed for use in the summaryGLSbands() function, specifically within the 
##indexS() portion

colnames(LDSCoutput$S)
##  [1] "cut1"            "cut2"            "cut3"            "cut4"           
##  [5] "cut5"            "cut6"            "cut7"            "cut8"           
##  [9] "cut9"            "cut10"           "cut11"           "cut12"          
## [13] "MDD"             "Wellbeing"       "Neuroticism"     "cut1_diffactor" 
## [17] "cut11_diffactor" "ptsd"            "scz"             "anxiety"        
## [21] "bipolar"         "adhd"            "autism"          "ocd"            
## [25] "an"              "compulsive"      "internalizing"   "neurodev"       
## [29] "thought"         "addict"          "oud"             "alch"           
## [33] "ts"              "cud"
## Here place the index number of the desired trait into the
## column vector of the indexS()function as in indexS(MATRIX = LDSCoutput$S, R = F)[1:12, XXX]
## Below is an example case in which MDD is the external trait of interest. 
## Note that the axes labels may be adjusted if for the plot

summaryGLSbands(OBJECT = subSV(LDSC_OBJECT = LDSCoutput,
                INDEXVALS = indexS(MATRIX = LDSCoutput$S, R = F)[1:12, 13 ], TYPE = "S_Stand"),
                PREDICTORS = PREDICTORS,
                INTERVALS = 40
                ,
                XLAB = "Neuroticism Severity (Z-Distribution)",
                YLAB = "Genetic Correlation w/ \n MDD",
                QUAD = F,
                XCOORDS = c(-1,2),
                YCOORDS = c(0,1))
##        betas        pvals         SE         Z
## b0 0.5439339 8.405734e-24 0.05407569 10.058752
## b1 0.1071030 1.922670e-04 0.02872190  3.728966

GTACCC of Negative Emotionality in Relation to Alternative Binary Depression Phenotypes from Cai et al., (2020)

Produces plots found in Extended Data Figure 4

load("LDSCoutput_Cutpoints.MinimalPhenotypes.RData")

###Predictors for GLS
PREDICTORS <- -1*qnorm(c(0.85556773,0.75867410,0.64767354,
                         0.53828939,0.43558773,0.34194621,
                         0.25818629,0.19038580,0.13159727, 
                         0.08571119,0.04947786,0.02166147))

## Here place the index number of the desired trait into the
## column vector of the indexS()function as in indexS(MATRIX = LDSCoutput$S, R = F)[1:12, XXX]
## Below is an example case in which DepAll is the external trait of interest. 
## Note that the axes labels may be adjusted if for the plot

colnames(LDSCoutput$S)
##  [1] "cut1"        "cut2"        "cut3"        "cut4"        "cut5"       
##  [6] "cut6"        "cut7"        "cut8"        "cut9"        "cut10"      
## [11] "cut11"       "cut12"       "MDD"         "DepAll"      "GPsy"       
## [16] "LifetimeMDD" "MDDrecur"    "Psypsy"      "SelfRepDep"  "ICD10Dep"   
## [21] "GPNoDep"
summaryGLSbands(OBJECT = subSV(LDSC_OBJECT = LDSCoutput,
                INDEXVALS = indexS(MATRIX = LDSCoutput$S, R = F)[1:12,14], TYPE = "S_Stand"),
                PREDICTORS = PREDICTORS,
                INTERVALS = 40
                ,
                XLAB = "Neuroticism Severity (Z-Distribution)",
                YLAB = "Genetic Correlation w/ \n DepAll",
                QUAD = F,
                XCOORDS = c(-1,2),
                YCOORDS = c(0,1))
##        betas        pvals         SE         Z
## b0 0.5556310 1.578388e-40 0.04168713 13.328597
## b1 0.0364221 5.948308e-02 0.01932611  1.884606

GTACCC of Negative Emotionality in Relation to Alternative Depression and Neuroticism Phenotypes

Produces plots found in Extended Data Figure 5 and 6

load("20250713_LDSCoutput_resubmit_phens_w.wrayreduced.RData")

###Predictors for GLS
PREDICTORS <- -1*qnorm(c(0.85556773,0.75867410,0.64767354,
                         0.53828939,0.43558773,0.34194621,
                         0.25818629,0.19038580,0.13159727, 
                         0.08571119,0.04947786,0.02166147))

## Here place the index number of the desired trait into the
## column vector of the indexS()function as in indexS(MATRIX = LDSCoutput$S, R = F)[1:12, XXX]
## Below is an example case in which the EHR depression phenotype from Adams et al., is the external trait of interest. 
## Note that the axes labels may be adjusted if for the plot

colnames(LDSCoutput$S)
##  [1] "cut1"         "cut2"         "cut3"         "cut4"         "cut5"        
##  [6] "cut6"         "cut7"         "cut8"         "cut9"         "cut10"       
## [11] "cut11"        "cut12"        "clin"         "ehr"          "no23"        
## [16] "no23noUKB"    "Quest"        "Hyde"         "Wray_Full"    "Howard"      
## [21] "neu_gupta"    "wray_reduced"
summaryGLSbands(OBJECT = subSV(LDSC_OBJECT = LDSCoutput,
                INDEXVALS = indexS(MATRIX = LDSCoutput$S, R = F)[1:12,13], TYPE = "S_Stand"),
                PREDICTORS = PREDICTORS,
                INTERVALS = 40
                ,
                XLAB = "Neuroticism Severity (Z-Distribution)",
                YLAB = "Genetic Correlation w/ \n EHR",
                QUAD = F,
                XCOORDS = c(-1,2),
                YCOORDS = c(0,1))
##         betas        pvals         SE         Z
## b0 0.72555584 5.948035e-60 0.04442835 16.330920
## b1 0.03472392 6.704034e-02 0.01896028  1.831404