setwd("/home/mchopra/Documents/PhD-Year1/PRS/")
library(RegParallel)
## Loading required package: doParallel
## Loading required package: foreach
## Loading required package: iterators
## Loading required package: parallel
## Loading required package: data.table
## Loading required package: stringr
## Loading required package: survival
## Loading required package: arm
## Loading required package: MASS
## Loading required package: Matrix
## Loading required package: lme4
## 
## arm (Version 1.13-1, built: 2022-8-25)
## Working directory is /home/mchopra/Documents/PhD-Year1/PRS
library(ggplot2)

GWAS of cellular proportion in GTEx samples and calculating PGS on UKBB data. Further, need to find if cellular proportion is associated with the distensibility. Can perform correlations as well.

aao_dist <- read.table("/home/mchopra/Documents/PhD-Year1/GWAS/norm_pheno_aao.txt", header = T, sep = ",")
dao_dist <-  read.table("/home/mchopra/Documents/PhD-Year1/GWAS/norm_pheno_dao.txt", header = T, sep = ",")

rownames(aao_dist) <- aao_dist[,1]
rownames(dao_dist) <- dao_dist[,1]

intersect <- intersect(rownames(aao_dist), rownames(dao_dist))

aao_dist <- aao_dist[intersect, ] 

dao_dist <- dao_dist[intersect,]

distensibility <- merge(aao_dist, dao_dist, by = "IID", all = TRUE)
#for VSMCI
##Case1: where PGS was calculated (without additional optimisation)
case1 <- read.table("cellp_pgs/VSMCI_nopheno_nocovar.all_score", header = T)

##distensibility_results:
#distensibility

#i want both of them in one table
#making rows
rownames(case1) <- case1[,1]
rownames(distensibility) <- distensibility[,1]

intersect <- intersect(rownames(case1), rownames(distensibility))

case1 <- case1[intersect, ] 

distensibility <- distensibility[intersect,]

case1_data <- merge(case1, distensibility, by = "IID", all = TRUE)
colnames(case1_data)
##  [1] "IID"                "FID"                "Pt_0.001"          
##  [4] "Pt_0.05"            "Pt_0.1"             "Pt_0.2"            
##  [7] "Pt_0.3"             "Pt_0.4"             "Pt_0.5"            
## [10] "Pt_1"               "FID.x"              "AAo_distensibility"
## [13] "FID.y"              "DAo_distensibility"
# Define the response variables
distens <- c("AAo_distensibility", "DAo_distensibility")

# Initialize a list to store the model results
macro_res_list <- list()

# Use a for loop to fit linear regression models for each response variable
for (ctypes in distens) {
  formula <- paste(ctypes, '~ as.numeric([*])')

  # Note: Replace 'as.numeric([*])' with the actual predictor variable
  # you want to use. I've used a placeholder here.

  macro_res <- RegParallel(
    data = case1_data,
    formula = formula,
    FUN = function(formula, data)
      lm(formula = formula, data = data),
    FUNtype = 'lm',
    variables = colnames(case1_data)[3:10],
    blocksize = 8,
    p.adjust = "BH",
    cores = 2,
    nestedParallel = FALSE,
    conflevel = 95
  )

  # Store the results in the list
  macro_res_list[[ctypes]] <- macro_res
}
## 
## ##############################
## #RegParallel
## ##############################
## System is:
## -- Linux
## Blocksize:
## -- 8
## Cores / Threads:
## -- 2
## Terms included in model:
## -- AAo_distensibility
## First 5 formulae:
## -- AAo_distensibility ~ as.numeric(Pt_0.001)
## -- AAo_distensibility ~ as.numeric(Pt_0.05)
## -- AAo_distensibility ~ as.numeric(Pt_0.1)
## -- AAo_distensibility ~ as.numeric(Pt_0.2)
## -- AAo_distensibility ~ as.numeric(Pt_0.3)
## Processing 8 formulae, batch 1 of 1
## -- index1: 1; index2: 8
## Done!
## 
## ##############################
## #RegParallel
## ##############################
## System is:
## -- Linux
## Blocksize:
## -- 8
## Cores / Threads:
## -- 2
## Terms included in model:
## -- DAo_distensibility
## First 5 formulae:
## -- DAo_distensibility ~ as.numeric(Pt_0.001)
## -- DAo_distensibility ~ as.numeric(Pt_0.05)
## -- DAo_distensibility ~ as.numeric(Pt_0.1)
## -- DAo_distensibility ~ as.numeric(Pt_0.2)
## -- DAo_distensibility ~ as.numeric(Pt_0.3)
## Processing 8 formulae, batch 1 of 1
## -- index1: 1; index2: 8
## Done!
PGS_results <- dplyr::bind_rows(macro_res_list, .id ="distens")

PGS_results
##                distens Variable                 Term        Beta StandardError
##                 <char>   <char>               <char>       <num>         <num>
##  1: AAo_distensibility Pt_0.001 as.numeric(Pt_0.001)    39.57533      25.36870
##  2: AAo_distensibility  Pt_0.05  as.numeric(Pt_0.05)   455.94933     221.70091
##  3: AAo_distensibility   Pt_0.1   as.numeric(Pt_0.1)   414.04947     323.68739
##  4: AAo_distensibility   Pt_0.2   as.numeric(Pt_0.2)   297.46584     475.53976
##  5: AAo_distensibility   Pt_0.3   as.numeric(Pt_0.3)   173.79724     615.44473
##  6: AAo_distensibility   Pt_0.4   as.numeric(Pt_0.4)   136.69057     728.77421
##  7: AAo_distensibility   Pt_0.5   as.numeric(Pt_0.5)   157.12139     843.76856
##  8: AAo_distensibility     Pt_1     as.numeric(Pt_1)   230.04442    1315.66968
##  9: DAo_distensibility Pt_0.001 as.numeric(Pt_0.001)    16.07625      25.67877
## 10: DAo_distensibility  Pt_0.05  as.numeric(Pt_0.05)   244.08982     224.41310
## 11: DAo_distensibility   Pt_0.1   as.numeric(Pt_0.1)    39.11156     327.64227
## 12: DAo_distensibility   Pt_0.2   as.numeric(Pt_0.2)  -397.13080     481.34023
## 13: DAo_distensibility   Pt_0.3   as.numeric(Pt_0.3)  -504.98462     622.94981
## 14: DAo_distensibility   Pt_0.4   as.numeric(Pt_0.4)  -747.25118     737.65805
## 15: DAo_distensibility   Pt_0.5   as.numeric(Pt_0.5)  -605.49052     854.05894
## 16: DAo_distensibility     Pt_1     as.numeric(Pt_1) -1007.42858    1331.71419
##              t          P            OR       ORlower       ORupper  P.adjust
##          <num>      <num>         <num>         <num>         <num>     <num>
##  1:  1.5600062 0.11876514  1.539379e+17  3.921672e-05  6.042544e+38 0.4750605
##  2:  2.0565966 0.03973055 1.038191e+198  2.014286e+09            NA 0.3178444
##  3:  1.2791647 0.20084539 6.597791e+179  1.977857e-96            NA 0.5355877
##  4:  0.6255331 0.53162429 1.540896e+129  0.000000e+00            NA 0.8611985
##  5:  0.2823929 0.77764348  3.014262e+75  0.000000e+00            NA 0.8611985
##  6:  0.1875623 0.85122059  2.311858e+59  0.000000e+00            NA 0.8611985
##  7:  0.1862138 0.85227787  1.725644e+68  0.000000e+00            NA 0.8611985
##  8:  0.1748497 0.86119850  8.072764e+99  0.000000e+00            NA 0.8611985
##  9:  0.6260522 0.53128375  9.590191e+06  1.330506e-15  6.912542e+28 0.6071814
## 10:  1.0876808 0.27674162 1.015921e+106  9.684669e-86 1.065701e+297 0.6071814
## 11:  0.1193727 0.90498057  9.681282e+16 1.248319e-262 7.508276e+295 0.9049806
## 12: -0.8250522 0.40934616 3.375093e-173  0.000000e+00            NA 0.6071814
## 13: -0.8106345 0.41757968 4.874880e-220  0.000000e+00            NA 0.6071814
## 14: -1.0130049 0.31106305  0.000000e+00  0.000000e+00            NA 0.6071814
## 15: -0.7089564 0.47835507 1.093476e-263  0.000000e+00            NA 0.6071814
## 16: -0.7564901 0.44935921  0.000000e+00  0.000000e+00            NA 0.6071814
# List of variables to fit the linear model
variables <- c("Pt_0.001","Pt_0.05","Pt_0.1","Pt_0.2","Pt_0.3","Pt_0.4","Pt_0.5","Pt_1")

# Empty list to store the model summaries
model_summaries <- list()

# Loop through each variable and fit a linear model
for (var in variables) {
  # Construct the formula dynamically
  formula <- as.formula(paste("AAo_distensibility ~", var))
  
  # Fit the linear model
  fit <- lm(formula, data = case1_data)
  
  # Store the summary of the model
  model_summary <- summary(fit)
  model_summaries[[var]] <- model_summary
}

# Access the summaries of each model
for (var in variables) {
  print(var)
  print(model_summaries[[var]])
}
## [1] "Pt_0.001"
## 
## Call:
## lm(formula = formula, data = case1_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9992 -0.6616  0.0042  0.6695  4.3157 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.039440   0.006904  -5.712 1.12e-08 ***
## Pt_0.001    39.575326  25.368697   1.560    0.119    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9841 on 47037 degrees of freedom
## Multiple R-squared:  5.174e-05,  Adjusted R-squared:  3.048e-05 
## F-statistic: 2.434 on 1 and 47037 DF,  p-value: 0.1188
## 
## [1] "Pt_0.05"
## 
## Call:
## lm(formula = formula, data = case1_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0013 -0.6617  0.0041  0.6717  4.3344 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 8.023e-03  1.966e-02   0.408   0.6832  
## Pt_0.05     4.559e+02  2.217e+02   2.057   0.0397 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9841 on 47037 degrees of freedom
## Multiple R-squared:  8.991e-05,  Adjusted R-squared:  6.865e-05 
## F-statistic:  4.23 on 1 and 47037 DF,  p-value: 0.03973
## 
## [1] "Pt_0.1"
## 
## Call:
## lm(formula = formula, data = case1_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0065 -0.6620  0.0037  0.6710  4.3239 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -0.009703   0.017500  -0.554    0.579
## Pt_0.1      414.049467 323.687387   1.279    0.201
## 
## Residual standard error: 0.9841 on 47037 degrees of freedom
## Multiple R-squared:  3.479e-05,  Adjusted R-squared:  1.353e-05 
## F-statistic: 1.636 on 1 and 47037 DF,  p-value: 0.2008
## 
## [1] "Pt_0.2"
## 
## Call:
## lm(formula = formula, data = case1_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0080 -0.6621  0.0042  0.6706  4.3240 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -0.01408    0.02794  -0.504    0.614
## Pt_0.2      297.46584  475.53976   0.626    0.532
## 
## Residual standard error: 0.9841 on 47037 degrees of freedom
## Multiple R-squared:  8.319e-06,  Adjusted R-squared:  -1.294e-05 
## F-statistic: 0.3913 on 1 and 47037 DF,  p-value: 0.5316
## 
## [1] "Pt_0.3"
## 
## Call:
## lm(formula = formula, data = case1_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0079 -0.6622  0.0045  0.6705  4.3233 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -0.02370    0.02737  -0.866    0.387
## Pt_0.3      173.79724  615.44473   0.282    0.778
## 
## Residual standard error: 0.9841 on 47037 degrees of freedom
## Multiple R-squared:  1.695e-06,  Adjusted R-squared:  -1.956e-05 
## F-statistic: 0.07975 on 1 and 47037 DF,  p-value: 0.7776
## 
## [1] "Pt_0.4"
## 
## Call:
## lm(formula = formula, data = case1_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0087 -0.6620  0.0043  0.6703  4.3228 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -0.02612    0.02811  -0.929    0.353
## Pt_0.4      136.69057  728.77421   0.188    0.851
## 
## Residual standard error: 0.9841 on 47037 degrees of freedom
## Multiple R-squared:  7.479e-07,  Adjusted R-squared:  -2.051e-05 
## F-statistic: 0.03518 on 1 and 47037 DF,  p-value: 0.8512
## 
## [1] "Pt_0.5"
## 
## Call:
## lm(formula = formula, data = case1_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0089 -0.6620  0.0043  0.6703  4.3227 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -0.02608    0.02853  -0.914    0.361
## Pt_0.5      157.12139  843.76856   0.186    0.852
## 
## Residual standard error: 0.9841 on 47037 degrees of freedom
## Multiple R-squared:  7.372e-07,  Adjusted R-squared:  -2.052e-05 
## F-statistic: 0.03468 on 1 and 47037 DF,  p-value: 0.8523
## 
## [1] "Pt_1"
## 
## Call:
## lm(formula = formula, data = case1_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0089 -0.6620  0.0043  0.6703  4.3227 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)   -0.02650    0.02794  -0.948    0.343
## Pt_1         230.04442 1315.66968   0.175    0.861
## 
## Residual standard error: 0.9841 on 47037 degrees of freedom
## Multiple R-squared:  6.5e-07,    Adjusted R-squared:  -2.061e-05 
## F-statistic: 0.03057 on 1 and 47037 DF,  p-value: 0.8612

the below shows that AAo_distensibility has relationship with VSMCI

fit <- lm(AAo_distensibility ~ Pt_0.05, data = case1_data)
summary(fit)
## 
## Call:
## lm(formula = AAo_distensibility ~ Pt_0.05, data = case1_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0013 -0.6617  0.0041  0.6717  4.3344 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 8.023e-03  1.966e-02   0.408   0.6832  
## Pt_0.05     4.559e+02  2.217e+02   2.057   0.0397 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9841 on 47037 degrees of freedom
## Multiple R-squared:  8.991e-05,  Adjusted R-squared:  6.865e-05 
## F-statistic:  4.23 on 1 and 47037 DF,  p-value: 0.03973
ggplot(case1_data, aes(x = AAo_distensibility, y = Pt_0.05)) +
  geom_point(color = "lightblue") +
  geom_smooth(method = "lm", color = "blue")
## `geom_smooth()` using formula = 'y ~ x'

OBSERVATION:: PRS of cellular proportion (VSMCI) shows some association WITH AAo_distensibility

CORRELATION:

library(corrplot)
## corrplot 0.92 loaded
## 
## Attaching package: 'corrplot'
## The following object is masked from 'package:arm':
## 
##     corrplot
x <- case1_data %>%  dplyr::select(AAo_distensibility, DAo_distensibility, Pt_0.001, Pt_0.05, Pt_0.1, Pt_0.2, Pt_0.3, Pt_0.4, Pt_0.5, Pt_1)

# Calculate correlation matrix
x <- cor(x)
x
##                    AAo_distensibility DAo_distensibility    Pt_0.001
## AAo_distensibility       1.0000000000       0.7081264374 0.007192755
## DAo_distensibility       0.7081264374       1.0000000000 0.002886615
## Pt_0.001                 0.0071927554       0.0028866154 1.000000000
## Pt_0.05                  0.0094822145       0.0050150606 0.294050659
## Pt_0.1                   0.0058979229       0.0005504087 0.257896076
## Pt_0.2                   0.0028842218      -0.0038041572 0.234180869
## Pt_0.3                   0.0013020679      -0.0037376812 0.223085204
## Pt_0.4                   0.0008648197      -0.0046707539 0.215927973
## Pt_0.5                   0.0008586022      -0.0032688679 0.213198984
## Pt_1                     0.0008062039      -0.0034880346 0.210561907
##                        Pt_0.05       Pt_0.1       Pt_0.2       Pt_0.3
## AAo_distensibility 0.009482214 0.0058979229  0.002884222  0.001302068
## DAo_distensibility 0.005015061 0.0005504087 -0.003804157 -0.003737681
## Pt_0.001           0.294050659 0.2578960755  0.234180869  0.223085204
## Pt_0.05            1.000000000 0.8812401402  0.801823167  0.764512573
## Pt_0.1             0.881240140 1.0000000000  0.911640737  0.870918970
## Pt_0.2             0.801823167 0.9116407375  1.000000000  0.957840848
## Pt_0.3             0.764512573 0.8709189704  0.957840848  1.000000000
## Pt_0.4             0.745580956 0.8510620822  0.938394868  0.980802448
## Pt_0.5             0.735226314 0.8396758219  0.926360802  0.969175951
## Pt_1               0.724614189 0.8287193460  0.914734182  0.957287731
##                           Pt_0.4        Pt_0.5          Pt_1
## AAo_distensibility  0.0008648197  0.0008586022  0.0008062039
## DAo_distensibility -0.0046707539 -0.0032688679 -0.0034880346
## Pt_0.001            0.2159279729  0.2131989836  0.2105619068
## Pt_0.05             0.7455809563  0.7352263138  0.7246141889
## Pt_0.1              0.8510620822  0.8396758219  0.8287193460
## Pt_0.2              0.9383948679  0.9263608017  0.9147341820
## Pt_0.3              0.9808024476  0.9691759512  0.9572877312
## Pt_0.4              1.0000000000  0.9887484918  0.9772486191
## Pt_0.5              0.9887484918  1.0000000000  0.9885711126
## Pt_1                0.9772486191  0.9885711126  1.0000000000
# Create correlation plot
corrplot(x, method = 'number', addCoef.col = "black", tl.col = "black")

corrplot(x, method="color",  
         type="upper", order="hclust", 
         addCoef.col = "black", # Add coefficient of correlation
         tl.col="black", tl.srt=45, #Text label color and rotation
         # hide correlation coefficient on the principal diagonal
         diag=FALSE 
         )

…. NO CORRELATION FROM THE ABOVE ANALYSIS. ….

Now, to the other way around where we are going to do the same but the pgs is calculated for distensibility in GTEx data. Therefore, the base data is UKB and target data is GTEx.

#PGS of distensibility 
pgs_aao <- read.table("dis_pgs/aao.all_score", header = T)

#cell proportion file
cell_prop <- read.csv("/home/mchopra/Documents/PhD-Year1/deconvolution/Deconvolution_results/results/CIBERSORTx_Results.csv")

#Round1 
rownames(pgs_aao) <- pgs_aao[,1]
rownames(cell_prop) <- cell_prop[,1]

intersect <- intersect(rownames(cell_prop), rownames(pgs_aao))


cell_prop <- cell_prop[intersect, ] 

pgs_aao <- pgs_aao[intersect,]

#changing the column1 name
names(pgs_aao)[1] <- "Mixture"

case2_aao <- merge(cell_prop, pgs_aao, by = "Mixture", all = TRUE)
# List of variables to fit the linear model
variables <- c("Pt_0.001","Pt_0.05","Pt_0.1","Pt_0.2","Pt_0.3","Pt_0.4","Pt_0.5","Pt_1")

# Empty list to store the model summaries
model_summaries <- list()

# Loop through each variable and fit a linear model
for (var in variables) {
  # Construct the formula dynamically
  formula <- as.formula(paste("VSMC_I ~", var))
  
  # Fit the linear model
  fit <- lm(formula, data = case2_aao)
  
  # Store the summary of the model
  model_summary <- summary(fit)
  model_summaries[[var]] <- model_summary
}

# Access the summaries of each model
for (var in variables) {
  print(var)
  print(model_summaries[[var]])
}
## [1] "Pt_0.001"
## 
## Call:
## lm(formula = formula, data = case2_aao)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27423 -0.06804  0.00040  0.07529  0.35457 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.25973    0.00767  33.862   <2e-16 ***
## Pt_0.001    -23.00487   29.02425  -0.793    0.428    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.113 on 385 degrees of freedom
## Multiple R-squared:  0.001629,   Adjusted R-squared:  -0.0009641 
## F-statistic: 0.6282 on 1 and 385 DF,  p-value: 0.4285
## 
## [1] "Pt_0.05"
## 
## Call:
## lm(formula = formula, data = case2_aao)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27594 -0.06560  0.00131  0.07568  0.35608 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.31097    0.03044  10.216   <2e-16 ***
## Pt_0.05     -329.17730  208.45690  -1.579    0.115    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1127 on 385 degrees of freedom
## Multiple R-squared:  0.006435,   Adjusted R-squared:  0.003855 
## F-statistic: 2.494 on 1 and 385 DF,  p-value: 0.1151
## 
## [1] "Pt_0.1"
## 
## Call:
## lm(formula = formula, data = case2_aao)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.26901 -0.06960 -0.00138  0.07530  0.35532 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.2889     0.0236  12.244   <2e-16 ***
## Pt_0.1      -335.6677   305.6644  -1.098    0.273    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1129 on 385 degrees of freedom
## Multiple R-squared:  0.003123,   Adjusted R-squared:  0.0005333 
## F-statistic: 1.206 on 1 and 385 DF,  p-value: 0.2728
## 
## [1] "Pt_0.2"
## 
## Call:
## lm(formula = formula, data = case2_aao)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27049 -0.06828 -0.00275  0.07488  0.35697 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.28046    0.01471  19.066   <2e-16 ***
## Pt_0.2      -494.00773  400.82453  -1.232    0.219    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1129 on 385 degrees of freedom
## Multiple R-squared:  0.00393,    Adjusted R-squared:  0.001343 
## F-statistic: 1.519 on 1 and 385 DF,  p-value: 0.2185
## 
## [1] "Pt_0.3"
## 
## Call:
## lm(formula = formula, data = case2_aao)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27475 -0.06914 -0.00143  0.07415  0.35713 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.27838    0.01217  22.884   <2e-16 ***
## Pt_0.3      -649.61454  476.74683  -1.363    0.174    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1128 on 385 degrees of freedom
## Multiple R-squared:  0.004799,   Adjusted R-squared:  0.002214 
## F-statistic: 1.857 on 1 and 385 DF,  p-value: 0.1738
## 
## [1] "Pt_0.4"
## 
## Call:
## lm(formula = formula, data = case2_aao)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27501 -0.06763 -0.00154  0.07356  0.35804 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.27493    0.01011   27.18   <2e-16 ***
## Pt_0.4      -737.08929  550.20681   -1.34    0.181    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1128 on 385 degrees of freedom
## Multiple R-squared:  0.00464,    Adjusted R-squared:  0.002055 
## F-statistic: 1.795 on 1 and 385 DF,  p-value: 0.1811
## 
## [1] "Pt_0.5"
## 
## Call:
## lm(formula = formula, data = case2_aao)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27476 -0.06786 -0.00198  0.07303  0.35830 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.724e-01  8.792e-03  30.981   <2e-16 ***
## Pt_0.5      -7.912e+02  6.117e+02  -1.293    0.197    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1128 on 385 degrees of freedom
## Multiple R-squared:  0.004326,   Adjusted R-squared:  0.00174 
## F-statistic: 1.673 on 1 and 385 DF,  p-value: 0.1967
## 
## [1] "Pt_1"
## 
## Call:
## lm(formula = formula, data = case2_aao)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27633 -0.06791 -0.00101  0.07379  0.35902 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.734e-01  8.566e-03  31.916   <2e-16 ***
## Pt_1        -1.228e+03  8.113e+02  -1.513    0.131    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1127 on 385 degrees of freedom
## Multiple R-squared:  0.005914,   Adjusted R-squared:  0.003332 
## F-statistic: 2.291 on 1 and 385 DF,  p-value: 0.131
# Define the response variables
celltypes <- c("VSMC_I", "VSMC_II")

# Initialize a list to store the model results
macro_res_list <- list()

# Use a for loop to fit linear regression models for each response variable
for (ctypes in celltypes) {
  formula <- paste(ctypes, '~ as.numeric([*])')

  # Note: Replace 'as.numeric([*])' with the actual predictor variable
  # you want to use. I've used a placeholder here.

  macro_res <- RegParallel(
    data = case2_aao,
    formula = formula,
    FUN = function(formula, data)
      lm(formula = formula, data = data),
    FUNtype = 'lm',
    variables = colnames(case2_aao)[24:31],
    blocksize = 8,
    p.adjust = "BH",
    cores = 2,
    nestedParallel = FALSE,
    conflevel = 95
  )

  # Store the results in the list
  macro_res_list[[ctypes]] <- macro_res
}
## 
## ##############################
## #RegParallel
## ##############################
## System is:
## -- Linux
## Blocksize:
## -- 8
## Cores / Threads:
## -- 2
## Terms included in model:
## -- VSMC_I
## First 5 formulae:
## -- VSMC_I ~ as.numeric(Pt_0.001)
## -- VSMC_I ~ as.numeric(Pt_0.05)
## -- VSMC_I ~ as.numeric(Pt_0.1)
## -- VSMC_I ~ as.numeric(Pt_0.2)
## -- VSMC_I ~ as.numeric(Pt_0.3)
## Processing 8 formulae, batch 1 of 1
## -- index1: 1; index2: 8
## Done!
## 
## ##############################
## #RegParallel
## ##############################
## System is:
## -- Linux
## Blocksize:
## -- 8
## Cores / Threads:
## -- 2
## Terms included in model:
## -- VSMC_II
## First 5 formulae:
## -- VSMC_II ~ as.numeric(Pt_0.001)
## -- VSMC_II ~ as.numeric(Pt_0.05)
## -- VSMC_II ~ as.numeric(Pt_0.1)
## -- VSMC_II ~ as.numeric(Pt_0.2)
## -- VSMC_II ~ as.numeric(Pt_0.3)
## Processing 8 formulae, batch 1 of 1
## -- index1: 1; index2: 8
## Done!
PGS_aao_results <- dplyr::bind_rows(macro_res_list, .id ="celltypes")
PGS_aao_results
##     celltypes Variable                 Term         Beta StandardError
##        <char>   <char>               <char>        <num>         <num>
##  1:    VSMC_I Pt_0.001 as.numeric(Pt_0.001)   -23.004867      29.02425
##  2:    VSMC_I  Pt_0.05  as.numeric(Pt_0.05)  -329.177298     208.45690
##  3:    VSMC_I   Pt_0.1   as.numeric(Pt_0.1)  -335.667745     305.66439
##  4:    VSMC_I   Pt_0.2   as.numeric(Pt_0.2)  -494.007725     400.82453
##  5:    VSMC_I   Pt_0.3   as.numeric(Pt_0.3)  -649.614535     476.74683
##  6:    VSMC_I   Pt_0.4   as.numeric(Pt_0.4)  -737.089294     550.20681
##  7:    VSMC_I   Pt_0.5   as.numeric(Pt_0.5)  -791.208475     611.74571
##  8:    VSMC_I     Pt_1     as.numeric(Pt_1) -1227.894285     811.31542
##  9:   VSMC_II Pt_0.001 as.numeric(Pt_0.001)    -5.092318      21.21451
## 10:   VSMC_II  Pt_0.05  as.numeric(Pt_0.05)   -16.609081     152.74319
## 11:   VSMC_II   Pt_0.1   as.numeric(Pt_0.1)  -108.040075     223.53348
## 12:   VSMC_II   Pt_0.2   as.numeric(Pt_0.2)   -74.219052     293.30778
## 13:   VSMC_II   Pt_0.3   as.numeric(Pt_0.3)   -88.330506     349.01711
## 14:   VSMC_II   Pt_0.4   as.numeric(Pt_0.4)    76.733855     402.77798
## 15:   VSMC_II   Pt_0.5   as.numeric(Pt_0.5)   133.929869     447.72593
## 16:   VSMC_II     Pt_1     as.numeric(Pt_1)   186.362007     594.25481
##              t         P            OR       ORlower       ORupper  P.adjust
##          <num>     <num>         <num>         <num>         <num>     <num>
##  1: -0.7926084 0.4284941  1.021206e-10  2.011987e-35  5.183242e+14 0.4284941
##  2: -1.5791144 0.1151310 1.096771e-143 3.992050e-321  3.012674e+34 0.2913632
##  3: -1.0981579 0.2728219 1.664756e-146  0.000000e+00 2.531319e+114 0.3117964
##  4: -1.2324788 0.2185224 2.852141e-215  0.000000e+00            NA 0.2913632
##  5: -1.3625986 0.1738056 7.516089e-283  0.000000e+00            NA 0.2913632
##  6: -1.3396586 0.1811466 7.692602e-321  0.000000e+00            NA 0.2913632
##  7: -1.2933617 0.1966615  0.000000e+00  0.000000e+00            NA 0.2913632
##  8: -1.5134610 0.1309830  0.000000e+00  0.000000e+00            NA 0.2913632
##  9: -0.2400394 0.8104274  6.143763e-03  5.377916e-21  7.018671e+15 0.9134665
## 10: -0.1087386 0.9134665  6.120232e-08 5.909207e-138 6.338793e+122 0.9134665
## 11: -0.4833284 0.6291373  1.198924e-47 6.408385e-238 2.243029e+143 0.9134665
## 12: -0.2530415 0.8003712  5.848915e-33 1.267767e-282 2.698430e+217 0.9134665
## 13: -0.2530836 0.8003387  4.350593e-39  0.000000e+00 5.277738e+258 0.9134665
## 14:  0.1905115 0.8490087  2.113925e+33  0.000000e+00            NA 0.9134665
## 15:  0.2991336 0.7649995  1.462188e+58  0.000000e+00            NA 0.9134665
## 16:  0.3136062 0.7539900  8.629615e+80  0.000000e+00            NA 0.9134665

##Now, the same above but with PGS of DAo distensibility

pgs_dao <- read.table("dis_pgs/dao.all_score", header = T)
cell_prop <- read.csv("/home/mchopra/Documents/PhD-Year1/deconvolution/Deconvolution_results/results/CIBERSORTx_Results.csv")
rownames(pgs_dao) <- pgs_dao[,1]
rownames(cell_prop) <- cell_prop[,1]

intersect <- intersect(rownames(cell_prop), rownames(pgs_dao))


cell_prop <- cell_prop[intersect, ] 

pgs_dao <- pgs_dao[intersect,]

#changing the column1 name
names(pgs_dao)[1] <- "Mixture"

case2_dao <- merge(cell_prop, pgs_dao, by = "Mixture", all = TRUE)

#colnames(case2_dao)
# Define the response variables
celltypes <- c("VSMC_I", "VSMC_II")

# Initialize a list to store the model results
macro_res_list <- list()

# Use a for loop to fit linear regression models for each response variable
for (ctypes in celltypes) {
  formula <- paste(ctypes, '~ as.numeric([*])')

  # Note: Replace 'as.numeric([*])' with the actual predictor variable
  # you want to use. I've used a placeholder here.

  macro_res <- RegParallel(
    data = case2_dao,
    formula = formula,
    FUN = function(formula, data)
      lm(formula = formula, data = data),
    FUNtype = 'lm',
    variables = colnames(case2_aao)[24:31],
    blocksize = 8,
    p.adjust = "BH",
    cores = 2,
    nestedParallel = FALSE,
    conflevel = 95
  )

  # Store the results in the list
  macro_res_list[[ctypes]] <- macro_res
}
## 
## ##############################
## #RegParallel
## ##############################
## System is:
## -- Linux
## Blocksize:
## -- 8
## Cores / Threads:
## -- 2
## Terms included in model:
## -- VSMC_I
## First 5 formulae:
## -- VSMC_I ~ as.numeric(Pt_0.001)
## -- VSMC_I ~ as.numeric(Pt_0.05)
## -- VSMC_I ~ as.numeric(Pt_0.1)
## -- VSMC_I ~ as.numeric(Pt_0.2)
## -- VSMC_I ~ as.numeric(Pt_0.3)
## Processing 8 formulae, batch 1 of 1
## -- index1: 1; index2: 8
## Done!
## 
## ##############################
## #RegParallel
## ##############################
## System is:
## -- Linux
## Blocksize:
## -- 8
## Cores / Threads:
## -- 2
## Terms included in model:
## -- VSMC_II
## First 5 formulae:
## -- VSMC_II ~ as.numeric(Pt_0.001)
## -- VSMC_II ~ as.numeric(Pt_0.05)
## -- VSMC_II ~ as.numeric(Pt_0.1)
## -- VSMC_II ~ as.numeric(Pt_0.2)
## -- VSMC_II ~ as.numeric(Pt_0.3)
## Processing 8 formulae, batch 1 of 1
## -- index1: 1; index2: 8
## Done!
PGS_dao_results <- dplyr::bind_rows(macro_res_list, .id ="celltypes")
PGS_dao_results
##     celltypes Variable                 Term        Beta StandardError
##        <char>   <char>               <char>       <num>         <num>
##  1:    VSMC_I Pt_0.001 as.numeric(Pt_0.001)   47.185620      30.74294
##  2:    VSMC_I  Pt_0.05  as.numeric(Pt_0.05)   61.641115     228.71367
##  3:    VSMC_I   Pt_0.1   as.numeric(Pt_0.1)  122.438727     320.83086
##  4:    VSMC_I   Pt_0.2   as.numeric(Pt_0.2)  197.473880     456.60313
##  5:    VSMC_I   Pt_0.3   as.numeric(Pt_0.3) -136.908191     552.22005
##  6:    VSMC_I   Pt_0.4   as.numeric(Pt_0.4) -309.183116     631.54397
##  7:    VSMC_I   Pt_0.5   as.numeric(Pt_0.5) -359.923836     723.33861
##  8:    VSMC_I     Pt_1     as.numeric(Pt_1) -619.771979     987.69885
##  9:   VSMC_II Pt_0.001 as.numeric(Pt_0.001)  -31.213829      22.46645
## 10:   VSMC_II  Pt_0.05  as.numeric(Pt_0.05)  -53.288117     167.04214
## 11:   VSMC_II   Pt_0.1   as.numeric(Pt_0.1)  -18.493419     234.37166
## 12:   VSMC_II   Pt_0.2   as.numeric(Pt_0.2) -151.469693     333.48659
## 13:   VSMC_II   Pt_0.3   as.numeric(Pt_0.3) -179.762443     403.25998
## 14:   VSMC_II   Pt_0.4   as.numeric(Pt_0.4)   -2.856055     461.41218
## 15:   VSMC_II   Pt_0.5   as.numeric(Pt_0.5)   73.441313     528.47044
## 16:   VSMC_II     Pt_1     as.numeric(Pt_1)  204.669520     721.69135
##                t         P            OR       ORlower       ORupper  P.adjust
##            <num>     <num>         <num>         <num>         <num>     <num>
##  1:  1.534843905 0.1256435  3.107811e+20  2.108692e-06  4.580322e+46 0.8043260
##  2:  0.269512160 0.7876800  5.893809e+26 1.227349e-168 2.830246e+221 0.8043260
##  3:  0.381630146 0.7029463  1.494389e+53 1.209880e-220            NA 0.8043260
##  4:  0.432484723 0.6656313  5.778516e+85  0.000000e+00            NA 0.8043260
##  5: -0.247923254 0.8043260  3.479589e-60  0.000000e+00            NA 0.8043260
##  6: -0.489567048 0.6247191 5.290280e-135  0.000000e+00            NA 0.8043260
##  7: -0.497586927 0.6190593 4.864789e-157  0.000000e+00            NA 0.8043260
##  8: -0.627490840 0.5307095 6.861983e-270  0.000000e+00            NA 0.8043260
##  9: -1.389353050 0.1655284  2.779754e-14  2.091850e-33  3.693874e+05 0.9950645
## 10: -0.319010027 0.7498919  7.198881e-24 4.685331e-166 1.106088e+119 0.9950645
## 11: -0.078906378 0.9371481  9.298446e-09 2.957286e-208 2.923664e+191 0.9950645
## 12: -0.454200249 0.6499407  1.650245e-66  0.000000e+00 1.207395e+218 0.9950645
## 13: -0.445773076 0.6560117  8.514570e-79  0.000000e+00            NA 0.9950645
## 14: -0.006189812 0.9950645  5.749514e-02  0.000000e+00            NA 0.9950645
## 15:  0.138969576 0.8895469  7.855194e+31  0.000000e+00            NA 0.9950645
## 16:  0.283597027 0.7768717  7.706254e+88  0.000000e+00            NA 0.9950645