These are the slides I used for a recent talk I gave (in Italian) during a conference on ME/CFS (video recording). These same slides can be downloaded here.

Abstract

I report here on the correlation between my ME/CFS symptoms and seven environmental parameters. I experience a regular improvement in my condition during the warmer months and I show here that the improvements directly correlate with air temperature and inversely correlate with air density and dry air density. No correlation could be found between my symptoms and relative humidity, absolute humidity, atmospheric pressure, and atmospheric particulate with a diameter (PM2.5). The correlation between symptoms and air density and dry air density might be driven by the correlation with temperature because density is strongly inversely correlated with the temperature of the air.

Introduction

My ME/CFS improves during summer, in the period of the year that goes from May/June to the end of September (when I am in Rome, Italy). I don’t know why. We know that ME/CFS tends to have an oscillating course in most patients (Chu L. et al. 2019), but the presence of a seasonal pattern in this patient population has not been investigated so far, to my knowledge. And yet, if you ask directly to patients, some of them say that they feel better in summer. Unfortunately, we don’t have scientific data on that, this is an area worth investigating with some carefully done surveys. I collected a functionality score (the higher, the more I function) for 190 days and I collected seven environmental parameters during these same days; I then calculated Spearman’s rank correlation coefficient between the score and each one of these parameters.

Data collection

I registered a functionality score three times a day (value from 3 to 7, the higher the better) for 199 days. I retrieved four environmental parameters for the same days (if available): mean temperature (T), mean atmospheric pressure (P), mean relative humidity (RH), and mean PM 2.5. While T, P, and RH were collected from Wunderground, PM 2.5 was retrieved from aqicn.

Data were collected during three separate periods, in three separate geographical regions: from March 1st to June 9th, 2021 in Arona (Tenerife, Spain); from February 1st to March 7th, 2020 in Rosario (Santa Fe, Argentina); from September 1st to October 31st, 2020 in Rome (Italy). These three periods were chosen for the following reasons: they are periods in which there is no big difference between the parameters of the air indoors and outdoors; they are periods in which my symptoms had their seasonal switch (I regularly improve during summer, from June to September, when I live in Rome). The closest weather station was used for each period, and when data were missing from one weather station, another one nearby was selected.

Calculation of derived parameters

The PM 2.5 value (particulate matter with a diameter equal to or below ) was converted from AQI (air quality index) to by solving the following equation with respect to C:

1)  

Where the intervals and are the ones collected in Table 1 and taken from (Mintz D, 2018). So for example, if we have a value of PM 2.5 of 55 AQI, we select and then we solve Eq. 1 with respect to C to get the corresponding value of PM 2.5 in ().

0.0–12.0	0–50
12.1–35.4	51–100
35.5–55.4	101–150
55.5–150.4	151–200
150.5–250.4	201–300
250.5–350.4	301–400
350.5–500.4	401–500

Relative humidity is defined as follows

where is the partial pressure of water vapor in the air and is the saturation vapor pressure with respect to water (the vapor pressure in the air in equilibrium with the surface of the water). We then have

Here and in what follows we express pressure in hPa, with , while is the temperature in °C and is the temperature in , with . Absolute humidity is simply given by

where is the mass of vapor and is the total volume occupied by the mixture (Ref, chapter 4). From the law of perfect gasses and considering that, according to Dalton’s law, in moisture, the partial pressure of each gas is the pressure it would have if it occupied the whole volume (Ref), we can also write

with . Then, by substituting and in , we have

2)   

It is also possible to show that the density of dry air (which is due to all the gasses that constitute the air, but water vapor) is given by

3)  

Where we assume . These calculations are reviewed in more detail in (Maccallini P, 2021). We can then find AD:

4)  

Statistical analysis

The correlation between each environmental parameter and the mean daily score of functionality was measured by the Spearman’s rank correlation coefficient, calculated as:

5)  

where S is the score of functionality and R(S) is its rank, Y is one of the environmental parameters and R(Y) is its rank. The statistical significance for the correlation parameters was calculated by assuming that the null distribution of Spearman’s coefficient is approximated by a symmetric law of the second type, therefore it is possible to show that

6)   

follows a t-student distribution with n – 2 degrees of freedom, where n is the number of measures available for S, Y (the number of days for which we have both S and Y) (Kendall MG, “The Advanced Theory of Statistics”, vol 2, par 31.19), (Maccallini P. 2021). The correction for multiple comparisons was performed according to the Benjamini-Hochberg test. All the calculations and plotting were performed by a custom Octave script (see at the bottom of the article). R was used for the correction for multiple comparisons, by using the following commands:

pval<-c(0.000357, 0.5356, 0.5151, 0.238, 0.00201, 0.00536, 0.1428)
p.adjust(pval, method = "hochberg") 

Parameter (daily mean)		t	p	p’
T (°C)	0.2563	3.637		0.0025
P (mbar)	-0.04522	-0.6206	0.5356	0.54
RH (%)	-0.04752	-0.6522	0.5151	0.54
AH ()	0.08595	1.183	0.238	0.54
AD ()	-0.223	-3.133		0.012
DAD ()	-0.2013	-2.817		0.027
PM2.5 ()	-0.1067	-1.472	0.1428	0.54

Table 2. Spearman correlation between my daily functionality score and several atmospheric parameters. T: temperature; P: pressure; RH: relative humidity; AH: absolute humidity; AD: air density; DAD: dry air density. P values after correction for multiple comparisons (Benjamini-Hochberg) are collected in the last column.

Results

The Spearman’s correlation coefficients between my score of functionality and each one of the environmental parameters considered are reported in Table 2. The only parameters to show a statistically significant correlation with my symptoms are temperature (T), air density (AD), and dry air density (DAD). While temperature is positively correlated with my score (the higher the temperature, the better I feel), air density and dry air density are inversely correlated (the less the density, the better I feel). The statistical significance of these three correlations survives after correction for multiple comparisons. No correlation was observed between my symptoms and atmospheric pressure, relative humidity (RH), absolute humidity (AH), and the concentration of atmospheric particulate with diameter (PM2.5). The distribution of good days (green) and bad days (red) on the plane AH-T and AD-T is reported in Figures 1, and 2. In Figures 3 and 4, it is caught an improvement and its relationship with temperature and air density, respectively. In Figures 5 and 6 a worsening of symptoms can be seen while temperature decreases and air density increases.

Figure 1. The mean daily temperature is on the x-axis while absolute humidity is on the y-axis. Points at the same relative humidity are indicated by the tick black curves. Green dots indicate days with a functionality score above 5; red dots indicate days with a score below 5; blue dots indicate days with a score of 5.

Figure 2. The mean daily temperature is on the x-axis while air density is on the y-axis. Green dots indicate days with a functionality score above 5; red dots indicate days with a score below 5; blue dots indicate days with a score of 5.

Figure 3. Improvement with increasing temperature. Green dots indicate days with a functionality score above 5; red dots indicate days with a score below 5; blue dots indicate days with a score of 5. Vertical dashed lines indicate the end-beginning of months.

Figure 4. Improvement with decreasing air density. Green dots indicate days with a well-being score above 5; red dots indicate days with a score below 5; blue dots indicate days with a score of 5. Vertical dashed lines indicate the end-beginning of months.

Figure 5. Worsening of symptoms with decreasing temperature. Green dots indicate days with a functionality score above 5; red dots indicate days with a score below 5; blue dots indicate days with a score of 5. Vertical dashed lines indicate the end-beginning of months.

Figure 6. Worsening of symptoms with increasing air density. Green dots indicate days with a functionality score above 5; red dots indicate days with a score below 5; blue dots indicate days with a score of 5. Vertical dashed lines indicate the end-beginning of months.

Discussion

The direct correlation between my degree of functionality and air temperature is the strongest (even if the correlation is low) and the most significant. Since air density and dry air density have a strong inverse dependence on temperature (see Eq. 3, 4), it seems likely that the correlation between my symptoms and AD and DAD is driven by the correlation with temperature. A more accurate statistical analysis (multinomial regression) might elucidate this point. Yet, further experiments (not shown here) have shown that increasing air temperature is not able to induce, by itself, an improvement in my symptoms, so this correlation might very well not indicate causation. On the other hand, a more complex experimental setting has allowed inducing, during winter, the improvement of symptoms that I would normally experience during summer. So I am currently designing better experiments in order to define the combination of parameters that are responsible for the seasonal pattern of my disease. The observation that the significant correlations here reported have very low coefficients is perhaps due to the fact that my score of functionality is particularly flat, and does not depict accurately the range of variability of my symptoms.

Notes

We calculated the daily mean absolute humidity by substituting in Eq. 2 the daily mean values for temperature, relative humidity, and pressure. In general, given a function x = x(t) and a function y = y(x), it is not true that the mean over t of y is given by y(mean of x over t), as it is easy to verify. But in our case, there is no big difference between the daily mean absolute humidity calculated by using Eq. 2 with daily mean values of T, P, and RH and the daily mean calculated from data taken every half an hour, as I verified for Rome, 2019 (see Figure 7). The same is true for air density (Figure 8). For further details see (Maccallini P. 2021)

Figure 7. Absolute humidity as a function of the day of the year for the year 2019, for the city of Rome (2019). The black line is the mean daily value calculated from measures of T, P, and RH taken every half an hour. The blue one is the value calculated from the mean daily values of T, P, and RH.

Figure 8. Air density as a function of the day of the year for the year 2019, for the city of Rome (2019). Black line: the daily mean of the values calculated from the actual measures for T, P, and RH. Blue line: the value calculated by using the daily means for T, P, and RH.

Supplementary material

The environmental measures and the symptom scores used for this study are available for download:

The script

Note: this script is pretty long for the following reasons: Octave is less concise than R in IO-related instructions; it plots and calculates more output than what is shown above; I calculate the p values without using a built-in function for that (I tend to avoid built-in functions, but I am now changing my mind).

% file name = Correlation
% date of creation = 08/11/2021
%-------------------------------------------------------------------------------
% It uses file correlation.xlsx
%-------------------------------------------------------------------------------
clear all
close all
pkg load io
missing = 999999999;                                                          % missing measure
%-------------------------------------------------------------------------------
% Physical parameters
%-------------------------------------------------------------------------------
R = 8.31;                                                                              % constant for ideal gasses (J/K*mole)
mu_v = 18.01528;                                                                % molecular weight of water (g/mole)
mu_d = 28.97;                                                                     % molecular weight of dry air (g/mole)
mu_N2 = 28.013;                                                                 % molecular weight of N2 (g/mole)
mu_O2 = 31.999;                                                                 % molecular weight of O2 (g/mole)
mu_CO2 = 44.010;                                                               % molecular weight of CO2 (g/mole)
mf_N2 = 0.78084;                                                                % molar fraction of N2
mf_O2 = 0.20946                                                                 % molar fraction of O2
mf_CO2 = 0.00039445                                                         % molar fraction of CO2
%-------------------------------------------------------------------------------
% We read the name of the sheet and its number of lines directly from the .xlsx file
%-------------------------------------------------------------------------------
file_name = 'Correlation.xlsx';
[filetype, sh_names, fformat, nmranges] = xlsfinfo (file_name);
sheet = char(sh_names(1,1));                                                % name of the first sheet
n_str = substr(char(sh_names(1,2)),5)                                   % number of lines as string
n_num = str2num(n_str)                                                       % number of lines as number
%-------------------------------------------------------------------------------
% We read column F (mean daily temperature, Celsius)
%-------------------------------------------------------------------------------
range = strcat('F2:F',n_str)
[num,tex,raw,limits] = xlsread (file_name, sheet, range);
t = num;
%-------------------------------------------------------------------------------
% We read column I (mean daily relative humidity, percentage)
%-------------------------------------------------------------------------------
range = strcat('I2:I',n_str)
[num,tex,raw,limits] = xlsread (file_name, sheet, range);
RH = num;
%-------------------------------------------------------------------------------
% We read column K (max daily pressure, mbar) 
%-------------------------------------------------------------------------------
range = strcat('K2:K',n_str)
[num,tex,raw,limits] = xlsread (file_name, sheet, range);
pmax = num;
%-------------------------------------------------------------------------------
% We read column L (min daily pressure, mbar)
%-------------------------------------------------------------------------------
range = strcat('L2:L',n_str)
[num,tex,raw,limits] = xlsread (file_name, sheet, range);
pmin = num;
%-------------------------------------------------------------------------------
% We read column M, N, O (score 1, 2, 3 respectively)
%-------------------------------------------------------------------------------
range = strcat('M2:M',n_str)
[num,tex,raw,limits] = xlsread (file_name, sheet, range);
score1 = num;
range = strcat('N2:N',n_str)
[num,tex,raw,limits] = xlsread (file_name, sheet, range);
score2 = num;
range = strcat('O2:O',n_str)
[num,tex,raw,limits] = xlsread (file_name, sheet, range);
score3 = num;
%-------------------------------------------------------------------------------
% We read column P (min daily PM 2.5, AQI)
%-------------------------------------------------------------------------------
range = strcat('P2:P',n_str)
[num,tex,raw,limits] = xlsread (file_name, sheet, range);
pm25 = num;
%-------------------------------------------------------------------------------
% We calculate mean daily pressure, absolute humidity, air density
%-------------------------------------------------------------------------------
for i = 1:n_num-1
  if ( t(i)!= missing )
    p(i) = ( pmin(i) + pmax(i) )/2;
    score(i) = ( score1(i) + score2(i) + score3(i) )/3;
    AH(i) = 1000*(mu_v/1000)*(RH(i)/(R*(273.15+t(i))))*(1.0016+3.15*10^(-6)*p(i)-0.074/p(i))*6.112*e^(17.62*t(i)/(243.12+t(i)));
    ep (i) = (RH(i)/100)*(1.0016+3.15*10^(-6)*p(i)-0.074/p(i))*6.112*e^(17.62*t(i)/(243.12+t(i)));
    AD_dry (i) = 1000*100*(p(i) - ep(i))*(mu_d/1000)/(R*(273.15+t(i)));
    AD (i) = AH(i) + AD_dry (i);
    N2D(i) = 1000*100*(p(i) - ep(i))*mf_N2*(mu_N2/1000)/(R*(273.15+t(i)));
    O2D(i) = 1000*100*(p(i) - ep(i))*mf_O2*(mu_O2/1000)/(R*(273.15+t(i)));
    CO2D(i) = 1000*100*(p(i) - ep(i))*mf_CO2*(mu_CO2/1000)/(R*(273.15+t(i)));
  else
    p(i) = missing;
    score(i) = missing;
    AH(i) = missing;
    AD_dry(i) = missing;
    AD (i) = missing;
    N2D(i) = missing;
    O2D(i) = missing;
    CO2D(i) = missing;
  endif
endfor
%-------------------------------------------------------------------------------
% We put all these arrays into a single array: this will help us to write a more
% compact script
%-------------------------------------------------------------------------------
k = 1;
for i = 1:n_num-1
  if ( t(i)!= missing )
    array(1,k) = t(i);
    array(2,k) = p(i);
    array(3,k) = pmin(i);
    array(4,k) = pmax(i);
    array(5,k) = RH(i);
    array(6,k) = AH(i);
    array(7,k) = AD(i);
    array(8,k) = N2D(i);
    array(9,k) = O2D(i);
    array(10,k) = CO2D(i);
    array(11,k) = score(i);
    k = k+1;
  endif
endfor
dim = size(array);
n = k - 1;
k = 1;
for i = 1:n_num-1
  if ( pm25(i)!= missing )
    array2(1,k) = score(i);
    array2(2,k) = pm25(i);
    k = k+1;
  endif
endfor
m = k - 1;
%clear p score AD N2D O2D CO2D pm25                      % we free some memory
%-------------------------------------------------------------------------------
% We calculate mean daily PM 2.5 density (microg/m^3) from daily PM 2.5 AQI
%-------------------------------------------------------------------------------
CLH = [0, 12, 12.1, 35.4, 35.5, 55.4, 55.5, 150.4, 150.5, 250.4, 250.5, 350.4, 350.5, 500.4];
ILH = [0, 50, 51, 100, 101, 150, 151, 200, 201, 300, 301, 400, 401, 500];
for k=1:m                                                                        % k defines the day
  for i=1:5005                                                                  % i defines the mass
    mass = (i-1)/10;
    for j=1:13                                                                    % j defines the index
      if ( and(mass>=CLH(j),mass<=CLH(j+1) ) )
        Itemp = ( ( ( ILH(j+1) - ILH(j) )/( CLH(j+1) - CLH(j) ) )*( mass - CLH(j) ) ) + ILH(j);
      endif
    endfor
    if (abs(Itemp - array2(2,k))<=1)
      array2(3,k) = mass;
    endif
  endfor
endfor
%-------------------------------------------------------------------------------
% We Compute Spearman’s rank correlation coefficient rho between each mean daily
% athmosferic parameter and the mean daily score for the same day. In order to do
% that, we need to create arrays without the days with missing measures. tvalue
% has a student's t-distribution with n-2 degrees of freedom (see par 31.19 of Kendall Vol. 2)
% We Compute also the coefficients a and b of the regression line
%-------------------------------------------------------------------------------
x = array(dim(1),1:n);
for h = 1:dim(1)-1
  y = array(h,1:n);
  rho(h) = spearman (x, y);
  tvalue(h) = rho(h)*sqrt( (n - 2)/(1 - rho(h)^2) );
  a(h) = cov(x,y)/var(x);
  b(h) = mean(y) - a(h)*mean(x);
endfor
clear x
x = array2(1,1:m);
for h = dim(1):dim(1)+1
  y = array2(h-dim(1)+2,1:m);
  rho(h) = spearman (x, y);
  tvalue(h) = rho(h)*sqrt( (n - 2)/(1 - rho(h)^2) );
  a(h) = cov(x,y)/var(x);
  b(h) = mean(y) - a(h)*mean(x);
endfor
%-------------------------------------------------------------------------------
% We Compute the p value (two-tailed) related to each rho 
%-------------------------------------------------------------------------------
ni = n - 2;                                                     % degrees of freedom of t-distribution
for h = 1:dim(1)+1
  s(h) = ni/(tvalue(h)^2 + ni);                       % betainc computes B(x,a,b)/B(a,b)
  pval(h) = 2*( 0.5*betainc(s(h), ni/2, 1/2) );  % so it is in fact what is usually referred to as
endfor                                                           % regularized incomplete beta function!
%-------------------------------------------------------------------------------
% We plot the physical parameters against the mean score
%-------------------------------------------------------------------------------
label_y = cell(dim(1)+1);
label_y(1) = 'daily mean temperature (°C)';
label_y(2) = 'daily mean pressure (mbar)';
label_y(3) = 'daily min pressure (mbar)';
label_y(4) = 'daily max pressure (mbar)';
label_y(5) = 'daily mean relative humidity (%)';
label_y(6) = 'daily mean absolute humidity (g/m^{3})';
label_y(7) = 'daily mean air density (g/m^{3})';
label_y(8) = 'daily mean N_{2} (g/m^{3})';
label_y(9) = 'daily mean O_{2} (g/m^{3})';
label_y(10) = 'daily mean CO_{2} (g/m^{3})';
label_y(11) = 'daily mean PM 2.5 (AQI)';
label_y(12) = 'daily mean PM 2.5 ({\mu}g/m^{3})';
for h = 1:dim(1)-1
  figure(h)
  for i=1:n
    plot (array(dim(1),i), array(h,i), 'o','markersize',5, 'markerfacecolor', 'r')
    hold on
  endfor
  xlabel ('daily mean symptom score', 'fontsize',15)
  ylabel (label_y(h), 'fontsize',15)
  rhostr = num2str(rho(h));
  pvalstr = num2str(pval(h));
  tvaluestr = num2str(tvalue(h));
  titlestr = strcat('Spearman correlation {\rho} = ', rhostr,', p = ', pvalstr, ', t = ', tvaluestr, ', delay = 0');
  title (titlestr,'fontsize',15)
  grid on
  val2 = max(array(h,:));
  val1 = min(array(h,:));
  axis ([3.5,6.5,val1,val2])
  plot ([3.5,6.5], [a(h)*3.5 + b(h), a(h)*6.5 + b(h)], '--b', 'linewidth', 2)
endfor
for h = dim(1):dim(1)+1
  figure(h)
  for i=1:m
    plot (array2(1,i), array2(h-dim(1)+2,i), 'o','markersize',5, 'markerfacecolor', 'r')
    hold on
  endfor
  xlabel ('daily mean symptom score', 'fontsize',15)
  ylabel (label_y(h), 'fontsize',15)
  rhostr = num2str(rho(h));
  pvalstr = num2str(pval(h));
  tvaluestr = num2str(tvalue(h));
  titlestr = strcat('Spearman correlation {\rho} = ', rhostr,', p = ', pvalstr, ', t = ', tvaluestr, ', delay = 0');
  title (titlestr,'fontsize',15)
  grid on
  val2 = max(array2(h-dim(1)+2,:));
  val1 = min(array2(h-dim(1)+2,:));
  axis ([3.5,6.5,val1,val2])
  plot ([3.5,6.5], [a(h)*3.5 + b(h), a(h)*6.5 + b(h)], '--b', 'linewidth', 2)
endfor
%-------------------------------------------------------------------------------
% We plot the chart temp-absolute humidity, with iso-relative humidity curves at
% an athmospheric pressure of 1,013.25 mbar
%-------------------------------------------------------------------------------
% molecular weight of water
mu_v = 18.01528;              % g/mol
% array of relative humidity (%)
RH2(1) = 10;
for i=2:20
  RH2(i) = RH2(i-1)+5;
endfor
% array of temperature (celsius)
tmin= 10;
tmax= 30;
for j=tmin:tmax
  t2(j) = j-1;
endfor
% array of pressure (hPa=mbar)
p_atm = 1013.25;
% array of absolute humidity (grams/m^3)----------------------------------------
for k=1:7
  for i=1:20
    for j=tmin:tmax
      AH2(k,i,j) = 1000*(mu_v/1000)*(RH2(i)/(R*(273.15+t2(j))))*(1.0016+3.15*10^(-6)*p_atm-0.074/p_atm)*6.112*e^(17.62*t2(j)/(243.12+t2(j)));
    endfor
  endfor
endfor
f = h;
f = f + 1;
figure(f)
for i=1:2:20
  plot (t2(tmin:tmax), AH2(k,i,tmin:tmax),'-k','LineWidth',2)
  label = num2str(RH2(i));
  label = strjoin ({label,'%'}, ' ');
  text (tmax+0.5, AH2(k,i,tmax)+0.5,label,'fontsize',16)
  hold on
endfor
axis ([tmin,tmax],"tic[xyz]", "label[xyz]")
xticks (t2(tmin:tmax))
yticks ([1:1:50])
xlabel ('temperature (°C)')
ylabel ('absolute humidity (g/{m^3})')
grid on
pressure_str = num2str(p_atm);
title_str = strcat('air pressure = ',pressure_str," mbar");
title (title_str)
%-------------------------------------------------------------------------------
% On the chart we have just plotted, we plot a red dot for a score of 4, a blue dot
% for a score of 5, a green dot for a score of 6. We consider score1, score2, score3
% instead of the mean scores. Scores are plotted assuming as coordinates temperature
% and absolute humidity of the day they refer to.
%-------------------------------------------------------------------------------
for i=1:n_num-1
  if ( t(i)!= missing )
    if (score1(i) == 4)
      plot (t(i),AH(i),'o','markersize',4,'markerfacecolor','r','markeredgecolor','r')
    endif
    if (score1(i) == 5)
      plot (t(i),AH(i),'o','markersize',4,'markerfacecolor','b','markeredgecolor','b')
    endif
    if (score1(i) == 6)
      plot (t(i),AH(i),'o','markersize',4,'markerfacecolor','g','markeredgecolor','g')
    endif
        if (score2(i) == 4)
      plot (t(i),AH(i),'o','markersize',4,'markerfacecolor','r','markeredgecolor','r')
    endif
    if (score2(i) == 5)
      plot (t(i),AH(i),'o','markersize',4,'markerfacecolor','b','markeredgecolor','b')
    endif
    if (score2(i) == 6)
      plot (t(i),AH(i),'o','markersize',4,'markerfacecolor','g','markeredgecolor','g')
    endif
        if (score3(i) == 4)
      plot (t(i),AH(i),'o','markersize',4,'markerfacecolor','r','markeredgecolor','r')
    endif
    if (score3(i) == 5)
      plot (t(i),AH(i),'o','markersize',4,'markerfacecolor','b','markeredgecolor','b')
    endif
    if (score3(i) == 6)
      plot (t(i),AH(i),'o','markersize',4,'markerfacecolor','g','markeredgecolor','g')
    endif
  endif
endfor
%-------------------------------------------------------------------------------
% On the plane temp-AD, we plot a red dot for a score of 4, a blue dot
% for a score of 5, a green dot for a score of 6. We consider score1, score2, score3
% instead of the mean scores. Scores are plotted assuming as coordinates temperature
% and absolute humidity of the day they refer to.
%-------------------------------------------------------------------------------
f = f + 1;
figure(f)
for i=1:n_num-1
  if ( t(i)!= missing )
    if (score1(i) == 4)
      plot (t(i),AD(i),'o','markersize',4,'markerfacecolor','r','markeredgecolor','r')
      hold on
    endif
    if (score1(i) == 5)
      plot (t(i),AD(i),'o','markersize',4,'markerfacecolor','b','markeredgecolor','b')
      hold on
    endif
    if (score1(i) == 6)
      plot (t(i),AD(i),'o','markersize',4,'markerfacecolor','g','markeredgecolor','g')
      hold on
    endif
    if (score2(i) == 4)
      plot (t(i),AD(i),'o','markersize',4,'markerfacecolor','r','markeredgecolor','r')
      hold on
    endif
    if (score2(i) == 5)
      plot (t(i),AD(i),'o','markersize',4,'markerfacecolor','b','markeredgecolor','b')
      hold on
    endif
    if (score2(i) == 6)
      plot (t(i),AD(i),'o','markersize',4,'markerfacecolor','g','markeredgecolor','g')
      hold on
    endif
    if (score3(i) == 4)
      plot (t(i),AD(i),'o','markersize',4,'markerfacecolor','r','markeredgecolor','r')
      hold on
    endif
    if (score3(i) == 5)
      plot (t(i),AD(i),'o','markersize',4,'markerfacecolor','b','markeredgecolor','b')
      hold on
    endif
    if (score3(i) == 6)
      plot (t(i),AD(i),'o','markersize',4,'markerfacecolor','g','markeredgecolor','g')
      hold on
    endif
  endif
endfor
axis ([tmin,tmax])
xlabel (label_y(1), 'fontsize',15)
ylabel (label_y(7), 'fontsize',15)
grid on
grid minor
%-------------------------------------------------------------------------------
% We plot the environmental data for the period March 1st to June 9th 2021
% (Canary Island)
%-------------------------------------------------------------------------------
% temperature
%-------------------------------------------------------------------------------
f = f + 1;
figure(f)
for i=1:102-1
  if ( t(i)!= missing )
    if (score(i) > 5)
      plot (i,t(i),'o','markersize',4,'markerfacecolor','g','markeredgecolor','g')
      hold on
    elseif (score(i) == 5 )
      plot (i,t(i),'o','markersize',4,'markerfacecolor','b','markeredgecolor','b')
      hold on
    else
      plot (i,t(i),'o','markersize',4,'markerfacecolor','r','markeredgecolor','r')
      hold on
    endif
  endif
endfor
plot([31,31],[min(array(1,:)),max(array(1,:))],'--k','linewidth',1)
text(31-20,max(array(1,:))-2,'March', 'fontsize',15)
hold on
plot([62,62],[min(array(1,:)),max(array(1,:))],'--k','linewidth',1)
text(62-20,max(array(1,:))-2,'April', 'fontsize',15)
hold on
plot([93,93],[min(array(1,:)),max(array(1,:))],'--k','linewidth',1)
text(93-20,max(array(1,:))-2,'May', 'fontsize',15)
grid on
grid minor
xlabel ('day', 'fontsize',15)
ylabel (label_y(1), 'fontsize',15)
title ('Canary Islands, Mar-Jun 2021; green: mean score > 5; blue: mean score = 5; red: mean score < 5','fontsize',15)
axis([1,101,min(array(1,:)),max(array(1,:))])
%-------------------------------------------------------------------------------
% absolute humidity
%-------------------------------------------------------------------------------
f = f + 1;
figure(f)
for i=1:102-1
  if ( t(i)!= missing )
    if (score(i) > 5)
      plot (i,AH(i),'o','markersize',4,'markerfacecolor','g','markeredgecolor','g')
      hold on
    elseif (score(i) == 5 )
      plot (i,AH(i),'o','markersize',4,'markerfacecolor','b','markeredgecolor','b')
      hold on
    else
      plot (i,AH(i),'o','markersize',4,'markerfacecolor','r','markeredgecolor','r')
      hold on
    endif
  endif
endfor
plot([31,31],[min(array(6,:)),max(array(6,:))],'--k','linewidth',1)
text(31-20,max(array(6,:))-2,'March', 'fontsize',15)
hold on
plot([62,62],[min(array(6,:)),max(array(6,:))],'--k','linewidth',1)
text(62-20,max(array(6,:))-2,'April', 'fontsize',15)
hold on
plot([93,93],[min(array(6,:)),max(array(6,:))],'--k','linewidth',1)
text(93-20,max(array(6,:))-2,'May', 'fontsize',15)
grid on
grid minor
xlabel ('day', 'fontsize',15)
ylabel (label_y(6), 'fontsize',15)
title ('Canary Islands, Mar-Jun 2021; green: mean score > 5; blue: mean score = 5; red: mean score < 5','fontsize',15)
axis([1,101,min(array(6,:)),max(array(6,:))])
%-------------------------------------------------------------------------------
% air density
%-------------------------------------------------------------------------------
f = f + 1;
figure(f)
for i=1:102-1
  if ( t(i)!= missing )
    if (score(i) > 5)
      plot (i,AD(i),'o','markersize',4,'markerfacecolor','g','markeredgecolor','g')
      hold on
    elseif (score(i) == 5 )
      plot (i,AD(i),'o','markersize',4,'markerfacecolor','b','markeredgecolor','b')
      hold on
    else
      plot (i,AD(i),'o','markersize',4,'markerfacecolor','r','markeredgecolor','r')
      hold on
    endif
  endif
endfor
plot([31,31],[min(array(7,:)),max(array(7,:))],'--k','linewidth',1)
text(31-20,max(array(7,:))-2,'March', 'fontsize',15)
hold on
plot([62,62],[min(array(7,:)),max(array(7,:))],'--k','linewidth',1)
text(62-20,max(array(7,:))-2,'April', 'fontsize',15)
hold on
plot([93,93],[min(array(7,:)),max(array(7,:))],'--k','linewidth',1)
text(93-20,max(array(7,:))-2,'May', 'fontsize',15)
grid on
grid minor
xlabel ('day', 'fontsize',15)
ylabel (label_y(7), 'fontsize',15)
title ('Canary Islands, Mar-Jun 2021; green: mean score > 5; blue: mean score = 5; red: mean score < 5','fontsize',15)
axis([1,101,min(array(7,:)),max(array(7,:))])
%-------------------------------------------------------------------------------
% We plot the environmental data for the period February to March 2020
% (Rosario, santa Fe)
%-------------------------------------------------------------------------------
% temperature
%-------------------------------------------------------------------------------
f = f + 1;
figure(f)
for i=103-1:138-1
  if ( t(i)!= missing )
    if (score(i) > 5)
      plot (i,t(i),'o','markersize',4,'markerfacecolor','g','markeredgecolor','g')
      hold on
    elseif (score(i) == 5 )
      plot (i,t(i),'o','markersize',4,'markerfacecolor','b','markeredgecolor','b')
      hold on
    else
      plot (i,t(i),'o','markersize',4,'markerfacecolor','r','markeredgecolor','r')
      hold on
    endif
  endif
endfor
plot([131,131],[min(array(1,:)),max(array(1,:))],'--k','linewidth',1)
text(131-20,max(array(1,:))-2,'Feb', 'fontsize',15)
hold on
text(131+10,max(array(1,:))-2,'March', 'fontsize',15)
grid on
grid minor
xlabel ('day', 'fontsize',15)
ylabel (label_y(1), 'fontsize',15)
title ('Rosario, Feb-Mar 2020; green: mean score > 5; blue: mean score = 5; red: mean score < 5','fontsize',15)
axis([103-1,138-1,min(array(1,:)),max(array(1,:))])
%-------------------------------------------------------------------------------
% absolute humidity
%-------------------------------------------------------------------------------
f = f + 1;
figure(f)
for i=103-1:138-1
  if ( t(i)!= missing )
    if (score(i) > 5)
      plot (i,AH(i),'o','markersize',4,'markerfacecolor','g','markeredgecolor','g')
      hold on
    elseif (score(i) == 5 )
      plot (i,AH(i),'o','markersize',4,'markerfacecolor','b','markeredgecolor','b')
      hold on
    else
      plot (i,AH(i),'o','markersize',4,'markerfacecolor','r','markeredgecolor','r')
      hold on
    endif
  endif
endfor
plot([131,131],[min(array(6,:)),max(array(6,:))],'--k','linewidth',1)
text(131-20,max(array(6,:))-2,'Feb', 'fontsize',15)
hold on
text(131+10,max(array(6,:))-2,'March', 'fontsize',15)
grid on
grid minor
xlabel ('day', 'fontsize',15)
ylabel (label_y(6), 'fontsize',15)
title ('Rosario, Feb-Mar 2020; green: mean score > 5; blue: mean score = 5; red: mean score < 5','fontsize',15)
axis([103-1,138-1,min(array(6,:)),max(array(6,:))])
%-------------------------------------------------------------------------------
% air density
%-------------------------------------------------------------------------------
f = f + 1;
figure(f)
for i=103-1:138-1
  if ( t(i)!= missing )
    if (score(i) > 5)
      plot (i,AD(i),'o','markersize',4,'markerfacecolor','g','markeredgecolor','g')
      hold on
    elseif (score(i) == 5 )
      plot (i,AD(i),'o','markersize',4,'markerfacecolor','b','markeredgecolor','b')
      hold on
    else
      plot (i,AD(i),'o','markersize',4,'markerfacecolor','r','markeredgecolor','r')
      hold on
    endif
  endif
endfor
plot([131,131],[min(array(7,:)),max(array(7,:))],'--k','linewidth',1)
text(131-20,max(array(7,:))-2,'Feb', 'fontsize',15)
hold on
text(131+10,max(array(7,:))-2,'March', 'fontsize',15)
grid on
grid minor
xlabel ('day', 'fontsize',15)
ylabel (label_y(7), 'fontsize',15)
title ('Rosario, Feb-Mar 2020; green: mean score > 5; blue: mean score = 5; red: mean score < 5','fontsize',15)
axis([103-1,138-1,min(array(7,:)),max(array(7,:))])
%-------------------------------------------------------------------------------
% We plot the environmental data for the period Sept to Oct 2020
% (Rome)
%-------------------------------------------------------------------------------
% temperature
%-------------------------------------------------------------------------------
f = f + 1;
figure(f)
for i=139-1:199-1
  if ( t(i)!= missing )
    if (score(i) > 5)
      plot (i,t(i),'o','markersize',4,'markerfacecolor','g','markeredgecolor','g')
      hold on
    elseif (score(i) == 5 )
      plot (i,t(i),'o','markersize',4,'markerfacecolor','b','markeredgecolor','b')
      hold on
    else
      plot (i,t(i),'o','markersize',4,'markerfacecolor','r','markeredgecolor','r')
      hold on
    endif
  endif
endfor
plot([168,168],[min(array(1,:)),max(array(1,:))],'--k','linewidth',1)
text(138+10,max(array(1,:))-2,'Sep', 'fontsize',15)
hold on
text(168+5,max(array(1,:))-2,'Oct', 'fontsize',15)
grid on
grid minor
xlabel ('day', 'fontsize',15)
ylabel (label_y(1), 'fontsize',15)
title ('Rome, Sep-Oct 2020; green: mean score > 5; blue: mean score = 5; red: mean score < 5','fontsize',15)
axis([139-1,199-1,min(array(1,:)),max(array(1,:))])
%-------------------------------------------------------------------------------
% absolute humidity
%-------------------------------------------------------------------------------
f = f + 1;
figure(f)
for i=139-1:199-1
  if ( t(i)!= missing )
    if (score(i) > 5)
      plot (i,AH(i),'o','markersize',4,'markerfacecolor','g','markeredgecolor','g')
      hold on
    elseif (score(i) == 5 )
      plot (i,AH(i),'o','markersize',4,'markerfacecolor','b','markeredgecolor','b')
      hold on
    else
      plot (i,AH(i),'o','markersize',4,'markerfacecolor','r','markeredgecolor','r')
      hold on
    endif
  endif
endfor
plot([168,168],[min(array(6,:)),max(array(6,:))],'--k','linewidth',1)
text(138+10,max(array(6,:))-2,'Sep', 'fontsize',15)
hold on
text(168+5,max(array(6,:))-2,'Oct', 'fontsize',15)
grid on
grid minor
xlabel ('day', 'fontsize',15)
ylabel (label_y(6), 'fontsize',15)
title ('Rome, Sep-Oct 2020; green: mean score > 5; blue: mean score = 5; red: mean score < 5','fontsize',15)
axis([139-1,199-1,min(array(6,:)),max(array(6,:))])
%-------------------------------------------------------------------------------
% air density
%-------------------------------------------------------------------------------
f = f + 1;
figure(f)
for i=139-1:199-1
  if ( t(i)!= missing )
    if (score(i) > 5)
      plot (i,AD(i),'o','markersize',4,'markerfacecolor','g','markeredgecolor','g')
      hold on
    elseif (score(i) == 5 )
      plot (i,AD(i),'o','markersize',4,'markerfacecolor','b','markeredgecolor','b')
      hold on
    else
      plot (i,AD(i),'o','markersize',4,'markerfacecolor','r','markeredgecolor','r')
      hold on
    endif
  endif
endfor
plot([168,168],[min(array(7,:)),max(array(7,:))],'--k','linewidth',1)
text(138+10,max(array(7,:))-2,'Sep', 'fontsize',15)
hold on
text(168+5,max(array(7,:))-2,'Oct', 'fontsize',15)
grid on
grid minor
xlabel ('day', 'fontsize',15)
ylabel (label_y(7), 'fontsize',15)
title ('Rome, Sep-Oct 2020; green: mean score > 5; blue: mean score = 5; red: mean score < 5','fontsize',15)
axis([139-1,199-1,min(array(7,:)),max(array(7,:))])

Abstract

I report on the association between variants with a minor allele frequency included between 1% and 5% and the phenotypes: self-reported Chronic Fatigue Syndrome and chronic asthenia not otherwise specified (ICD10: R53). The analysis was performed on a metadata file generated from the UK Biobank database. Patients with self-reported CFS (both sexes) showed a significant increase (p = ) in the frequency of the alternative allele of variant rs148723539, 10:131677233:G:A on h19, an SNP of the intronic region of gene EBF3. Female patients with R53 showed a significant association (p = ) with variant rs182581502, 7:56492485:G:T on h19, which belongs to the intronic region of pseudogene RP13-492C18.2 (gene_id: ENSG00000237268.2). This analysis completes a previous one (here) performed only on common variants (minor allele frequency above 5%).

Introduction

In a previous blog post (here), I reported on the association between self-reported Chronic Fatigue Syndrome and a region of chromosome 13, in particular variant rs11147812. This result was obtained by performing a set of analyses on a metadata file generated by Neale’s Lab from the UK Biobank dataset. The analysis included about 14 million variants, 1659 patients (451 males), and 300k controls. I focused on common variants (minor allele frequency above 0.05) and I used a cut-off for statistical significance of . I performed a similar analysis on chronic asthenia not otherwise specified (ICD10: R53) and I found an association between male patients and variant rs62092652.

Here, I report on the results of the same set of analyses for variants with a minor allele frequency below 0.05 and above 0.01. I used a more stringent cut-off for statistical significance of (read this blog post for an introduction to the p value in GWAS).

Methods

I considered again the phenotypes: self-reported CFS and R53. I filtered out variants too far from the Hardy-Weinberg equilibrium (in the total UK Biobank sample) by using a cut-off of when testing the null hypothesis: the alleles are in HWE. I filtered out variants with expected minor allele count below 25, in cases.

Results

When considering self-reported CFS, the only statistically significant association was between variant rs148723539 (10:131677233:G:A on h19) and the group of patients that included both sexes (Table 1). In particular, the association is between the alternative allele and patients (in other words, patients have a higher frequency of the alternative allele A than controls). This is an intronic variant of gene EBF3. No association was found after stratification by sex.

When considering chronic asthenia not otherwise specified (ICD10: R53), the analysis for the whole sample (both sexes) gave no associations. When considering females only, variant rs182581502 (7:56492485:G:T on h19) reached statistical significance with the alternative allele more frequent in patients than in controls (Table 2). This is an intronic variant of pseudogene RP13-492C18.2 (gene_id: ENSG00000237268.2). No associations were found for males only.

Variant rs182581502 is associated with changes in gene expression of gene CHCHD2 in several tissues (see this page). The alternative allele (the one more frequent in female R53 patients) is associated with an increase in the expression of CHCHD2 in the hippocampus and basal ganglia, for instance. No effects on gene expression have been found for rs148723539.

A concise summary of the results of my analysis of the UK Biobank database is presented in the following table.