• Description:
    • Region 1
    • Region 2
    • Region 3
    • Region 4
    • Region 5

Last updated: 2024-06-27

Checks: 7 0

Knit directory: survival-data-analysis/

This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.


Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

The command set.seed(20240324) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version d80d61c. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .DS_Store
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/

Unstaged changes:
    Deleted:    analysis/logistic_gwas_asthma.Rmd
    Deleted:    analysis/susie_asthma_result.Rmd

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.


These are the previous versions of the repository in which changes were made to the R Markdown (analysis/susie_asthma_result1.Rmd) and HTML (docs/susie_asthma_result1.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File Version Author Date Message
Rmd d80d61c yunqiyang0215 2024-06-27 wflow_publish("analysis/susie_asthma_result1.Rmd")
html 85ef406 yunqiyang0215 2024-06-20 Build site.
Rmd f5dfe15 yunqiyang0215 2024-06-20 wflow_publish("analysis/susie_asthma_result1.Rmd")
html fc8e7c4 yunqiyang0215 2024-06-20 Build site.
Rmd cb1583e yunqiyang0215 2024-06-20 wflow_publish("analysis/susie_asthma_result1.Rmd")

Description:

Coxph Susie result on all asthma/ AOA/ COA in UKBiobank.

library(survival)
library(susieR)
devtools::load_all("/Users/nicholeyang/Downloads/logisticsusie")
ℹ Loading logisticsusie

Region 1

Strong signals for COA, marginal significant for AOA. rs61894547 was the most significant SNP reported by Carole’s paper, but not the most significant one in my result. However, have the largest PIP.

1. All asthma cases

region = "chr11_75500001_77400000"
res = readRDS(paste0("/Users/nicholeyang/Downloads/survivalsusie/result/asthma_self_report/result/all/fit.susie.", region, ".rds"))
gwas = readRDS(paste0("/Users/nicholeyang/downloads/survivalsusie/result/gwas_surv/all_gwas_", region, ".rds"))
fit = res[[1]]
X = res[[2]]
print(res[[3]])
      user     system    elapsed 
116070.891  65190.274   9564.771 
pip <- logisticsusie:::get_pip(fit$alpha)
effect_estimate <- colSums(fit$alpha * fit$mu)
pip.sorted = sort(pip, decreasing = TRUE)
pip.sorted[1:10]
 [1] 0.62563969 0.36681900 0.21652826 0.18795848 0.16064130 0.11232667
 [7] 0.10899708 0.07188830 0.06583735 0.05808435
class(fit) = "susie"
cs <- susie_get_cs(fit, X)
cs
$cs
$cs$L2
[1]  829 1000

$cs$L1
 [1]  943  951  952  954  961  964  965  968  979 1001


$purity
   min.abs.corr mean.abs.corr median.abs.corr
L2    0.9428971     0.9428971       0.9428971
L1    0.9003599     0.9619019       0.9510408

$cs_index
[1] 2 1

$coverage
[1] 0.9890498 0.9574096

$requested_coverage
[1] 0.95
par(mfrow = c(1,2))
snps1 = colnames(X)[cs$cs$L1]
colors <- ifelse(rownames(gwas) %in% snps1, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.8, pch = 20, main = "CS 1")

snps2 = colnames(X)[cs$cs$L2]
colors <- ifelse(rownames(gwas) %in% snps2, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.8, pch = 20, main = "CS 2")

Version Author Date
fc8e7c4 yunqiyang0215 2024-06-20
cbind(gwas[rownames(gwas) %in% snps1, ], pip[sort(cs$cs$L1)])
                   MAF missing.rate  p.value.spa p.value.norm      Stat
rs61893460_A 0.4520203            0 8.124552e-38 7.706256e-38  1613.997
rs7126418_T  0.4519984            0 1.379899e-37 1.309905e-37  1607.742
rs7110818_T  0.4512180            0 1.195983e-37 1.134530e-37  1608.185
rs7114362_T  0.4968866            0 9.059745e-36 8.818957e-36 -1570.311
rs7936070_T  0.4766971            0 1.503789e-37 1.452115e-37  1612.065
rs7936312_T  0.4766166            0 1.251399e-37 1.208159e-37  1613.890
rs7936323_A  0.4765950            0 1.070192e-37 1.033056e-37  1615.204
rs7936434_C  0.4768852            0 2.342714e-37 2.263324e-37  1607.740
rs11236791_A 0.4518935            0 8.407000e-38 7.973748e-38  1613.119
rs11236797_A 0.4510802            0 5.710907e-38 5.409022e-38  1616.420
                  Var         z           
rs61893460_A 15755.24  12.85849 0.06583735
rs7126418_T  15733.73  12.81742 0.04220489
rs7110818_T  15715.07  12.82856 0.04494202
rs7114362_T  15815.14 -12.48674 0.02049406
rs7936070_T  15838.20  12.80942 0.16064130
rs7936312_T  15838.79  12.82369 0.18795848
rs7936323_A  15834.63  12.83582 0.21652826
rs7936434_C  15838.51  12.77494 0.10899708
rs11236791_A 15744.57  12.85586 0.05808435
rs11236797_A 15735.62  12.88583 0.07188830
cbind(gwas[rownames(gwas) %in% snps2, ], pip[sort(cs$cs$L2)])
                    MAF missing.rate  p.value.spa p.value.norm     Stat
rs61894547_T 0.05155540            0 3.916273e-24 8.543620e-25 570.9379
rs55646091_A 0.05086819            0 5.110427e-25 9.672793e-26 572.9645
                  Var        z          
rs61894547_T 3083.681 10.28145 0.3668190
rs55646091_A 2983.742 10.48931 0.6256397

2. COA

region = "chr11_75500001_77400000"
res = readRDS(paste0("/Users/nicholeyang/Downloads/survivalsusie/result/asthma_self_report/result/coa/fit.susie.", region, ".rds"))
gwas = readRDS(paste0("/Users/nicholeyang/downloads/survivalsusie/result/gwas_surv/coa_gwas_", region, ".rds"))
fit = res[[1]]
X = res[[2]]

pip <- logisticsusie:::get_pip(fit$alpha)
effect_estimate <- colSums(fit$alpha * fit$mu)

class(fit) = "susie"
cs <- susie_get_cs(fit, X)
cs
$cs
$cs$L1
[1]  943  951  952  961  964  965  968  979 1001

$cs$L2
[1]  829 1000


$purity
   min.abs.corr mean.abs.corr median.abs.corr
L1    0.9491341     0.9718758       0.9511525
L2    0.9427185     0.9427185       0.9427185

$cs_index
[1] 1 2

$coverage
[1] 0.9664054 0.9999830

$requested_coverage
[1] 0.95
par(mfrow = c(1,2))
snps1 = colnames(X)[cs$cs$L1]
colors <- ifelse(rownames(gwas) %in% snps1, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.8, pch = 20, main = "CS 1")

snps2 = colnames(X)[cs$cs$L2]
colors <- ifelse(rownames(gwas) %in% snps2, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.8, pch = 20, main = "CS 2")

Version Author Date
fc8e7c4 yunqiyang0215 2024-06-20
cbind(gwas[rownames(gwas) %in% snps1, ], pip[sort(cs$cs$L1)])
                   MAF missing.rate  p.value.spa p.value.norm     Stat      Var
rs61893460_A 0.4514041            0 9.768368e-35 7.292156e-35 776.1724 3970.740
rs7126418_T  0.4513778            0 2.682089e-34 2.019865e-34 770.4486 3965.300
rs7110818_T  0.4505984            0 4.350409e-34 3.284013e-34 767.5109 3960.668
rs7936070_T  0.4760599            0 4.340150e-34 3.511120e-34 770.3164 3993.243
rs7936312_T  0.4759791            0 3.808504e-34 3.078221e-34 771.0068 3993.379
rs7936323_A  0.4759578            0 3.122270e-34 2.520459e-34 771.9348 3992.346
rs7936434_C  0.4762484            0 4.409483e-34 3.568758e-34 770.2264 3993.180
rs11236791_A 0.4512776            0 5.639224e-35 4.187165e-35 778.7052 3967.866
rs11236797_A 0.4504619            0 7.104772e-35 5.269381e-35 777.3307 3965.725
                    z           
rs61893460_A 12.31750 0.11519132
rs7126418_T  12.23505 0.05051082
rs7110818_T  12.19552 0.03266226
rs7936070_T  12.19007 0.10260708
rs7936312_T  12.20079 0.11384264
rs7936323_A  12.21706 0.13434046
rs7936434_C  12.18874 0.10180916
rs11236791_A 12.36217 0.18323304
rs11236797_A 12.34368 0.14755555
cbind(gwas[rownames(gwas) %in% snps2, ], pip[sort(cs$cs$L2)])
                    MAF missing.rate  p.value.spa p.value.norm     Stat
rs61894547_T 0.05133647            0 5.289826e-31 3.923856e-34 338.8359
rs55646091_A 0.05064274            0 1.900986e-30 1.698689e-33 329.9953
                  Var        z          
rs61894547_T 773.7703 12.18101 0.8735478
rs55646091_A 748.6081 12.06092 0.1293317

3. AOA

region = "chr11_75500001_77400000"
res = readRDS(paste0("/Users/nicholeyang/Downloads/survivalsusie/result/asthma_self_report/result/aoa/fit.susie.", region, ".rds"))
gwas = readRDS(paste0("/Users/nicholeyang/downloads/survivalsusie/result/gwas_surv/aoa_gwas_", region, ".rds"))
fit = res[[1]]
X = res[[2]]

pip <- logisticsusie:::get_pip(fit$alpha)
effect_estimate <- colSums(fit$alpha * fit$mu)

class(fit) = "susie"
cs <- susie_get_cs(fit, X)
cs
$cs
$cs$L1
 [1]  927  943  951  952  954  961  964  965  968  975  979  990  998 1001 1011


$purity
   min.abs.corr mean.abs.corr median.abs.corr
L1    0.8866167     0.9479983       0.9486556

$cs_index
[1] 1

$coverage
[1] 0.9786045

$requested_coverage
[1] 0.95
snps1 = colnames(X)[cs$cs$L1]
colors <- ifelse(rownames(gwas) %in% snps1, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.8, pch = 20, main = "CS 1")

Version Author Date
fc8e7c4 yunqiyang0215 2024-06-20
cbind(gwas[rownames(gwas) %in% snps1, ], pip[sort(cs$cs$L1)])
                   MAF missing.rate  p.value.spa p.value.norm      Stat
rs2212434_T  0.4459042            0 1.180847e-08 1.177059e-08  544.8959
rs61893460_A 0.4500073            0 1.006457e-08 1.003220e-08  548.8137
rs7126418_T  0.4499901            0 1.225724e-08 1.221881e-08  545.2186
rs7110818_T  0.4492170            0 9.716273e-09 9.684795e-09  548.6854
rs7114362_T  0.4987621            0 2.299994e-09 2.292086e-09 -573.7888
rs7936070_T  0.4747088            0 9.528090e-09 9.499610e-09  551.3555
rs7936312_T  0.4746268            0 8.941648e-09 8.914703e-09  552.3982
rs7936323_A  0.4746036            0 8.621435e-09 8.595334e-09  552.9179
rs7936434_C  0.4748971            0 1.240088e-08 1.236504e-08  547.0490
rs4494327_T  0.4991553            0 5.440181e-09 5.423273e-09 -561.0342
rs11236791_A 0.4498759            0 1.375619e-08 1.371371e-08  543.5143
rs10160518_G 0.4992504            0 4.477610e-09 4.463369e-09 -564.2245
rs2155219_T  0.4996877            0 7.540926e-09 7.518396e-09 -556.0470
rs11236797_A 0.4490612            0 1.245507e-08 1.241597e-08  545.0017
rs7930763_A  0.4986349            0 7.366281e-09 7.344218e-09 -555.4006
                  Var         z           
rs2212434_T  9128.890  5.703016 0.04631248
rs61893460_A 9173.030  5.730184 0.05289074
rs7126418_T  9160.160  5.696645 0.04455621
rs7110818_T  9149.649  5.736159 0.05423549
rs7114362_T  9220.138 -5.975626 0.19137559
rs7936070_T  9228.392  5.739429 0.05519014
rs7936312_T  9228.713  5.750183 0.05833582
rs7936323_A  9226.294  5.756348 0.06022704
rs7936434_C  9228.343  5.694615 0.04388581
rs4494327_T  9249.114 -5.833635 0.08829031
rs11236791_A 9166.331  5.676926 0.04027700
rs10160518_G 9251.549 -5.866036 0.10499977
rs2155219_T  9258.281 -5.778915 0.06624745
rs11236797_A 9161.653  5.693914 0.04374122
rs7930763_A  9224.176 -5.782858 0.06706628
rm()

Region 2

No GWAS significant signal for COA, marginal significant for AOA.

Result: for all asthma and COA, no CS found. For AOA, there is one CS. rs56389811_T was the top significant signal reported by Carole’s paper, and also the top significant one found in AOA survival gwas. PIP = 0.2

1. All asthma cases

region = "chr12_46000001_48700000"
res = readRDS(paste0("/Users/nicholeyang/Downloads/survivalsusie/result/asthma_self_report/result/all/fit.susie.", region, ".rds"))
gwas = readRDS(paste0("/Users/nicholeyang/downloads/survivalsusie/result/gwas_surv/all_gwas_", region, ".rds"))
fit = res[[1]]
X = res[[2]]
print(res[[3]])
     user    system   elapsed 
58385.924 32800.712  4853.472 
pip <- logisticsusie:::get_pip(fit$alpha)
effect_estimate <- colSums(fit$alpha * fit$mu)
pip.sorted = sort(pip, decreasing = TRUE)
pip.sorted[1:10]
 [1] 0.16103664 0.15329430 0.13547046 0.12479556 0.12383535 0.06007752
 [7] 0.05636508 0.04679767 0.04200849 0.03698004
class(fit) = "susie"
cs <- susie_get_cs(fit, X)
cs
$cs
NULL

$coverage
NULL

$requested_coverage
[1] 0.95

2. COA

region = "chr12_46000001_48700000"
res = readRDS(paste0("/Users/nicholeyang/Downloads/survivalsusie/result/asthma_self_report/result/coa/fit.susie.", region, ".rds"))
gwas = readRDS(paste0("/Users/nicholeyang/downloads/survivalsusie/result/gwas_surv/coa_gwas_", region, ".rds"))
fit = res[[1]]
X = res[[2]]

pip <- logisticsusie:::get_pip(fit$alpha)
effect_estimate <- colSums(fit$alpha * fit$mu)

class(fit) = "susie"
cs <- susie_get_cs(fit, X)
cs
$cs
NULL

$coverage
NULL

$requested_coverage
[1] 0.95

3. AOA

region = "chr12_46000001_48700000"
res = readRDS(paste0("/Users/nicholeyang/Downloads/survivalsusie/result/asthma_self_report/result/aoa/fit.susie.", region, ".rds"))
gwas = readRDS(paste0("/Users/nicholeyang/downloads/survivalsusie/result/gwas_surv/aoa_gwas_", region, ".rds"))
fit = res[[1]]
X = res[[2]]

pip <- logisticsusie:::get_pip(fit$alpha)
effect_estimate <- colSums(fit$alpha * fit$mu)

class(fit) = "susie"
cs <- susie_get_cs(fit, X)
cs
$cs
$cs$L1
[1] 760 785 787 808 812 814 818 828 829


$purity
   min.abs.corr mean.abs.corr median.abs.corr
L1    0.9091167     0.9655223       0.9615681

$cs_index
[1] 1

$coverage
[1] 0.954442

$requested_coverage
[1] 0.95
snps1 = colnames(X)[cs$cs$L1]
colors <- ifelse(rownames(gwas) %in% snps1, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.8, pch = 20, main = "CS 1")

Version Author Date
fc8e7c4 yunqiyang0215 2024-06-20
cbind(gwas[rownames(gwas) %in% snps1, ], pip[sort(cs$cs$L1)])
                   MAF missing.rate  p.value.spa p.value.norm      Stat
rs73107980_T 0.2409617            0 2.672868e-09 2.622406e-09 -487.2183
rs73107993_T 0.2486728            0 7.335992e-10 7.177201e-10 -512.9498
rs55726902_A 0.2423548            0 9.548414e-10 9.338834e-10 -504.9693
rs11168244_T 0.2389714            0 2.196870e-10 2.135056e-10 -520.8690
rs11168245_G 0.2391581            0 2.026872e-10 1.969276e-10 -522.1264
rs11168246_A 0.2092811            0 2.426769e-09 2.369254e-09 -446.0817
rs56389811_T 0.2389045            0 1.324647e-10 1.284605e-10 -527.2831
rs11168250_T 0.2389891            0 1.395104e-10 1.353413e-10 -527.3332
rs11168252_A 0.2392910            0 1.480046e-10 1.436179e-10 -525.1814
                  Var         z           
rs73107980_T 6697.022 -5.953642 0.01450955
rs73107993_T 6929.274 -6.162132 0.04253712
rs55726902_A 6807.408 -6.120319 0.03359861
rs11168244_T 6725.614 -6.351298 0.13217838
rs11168245_G 6731.772 -6.363718 0.14305372
rs11168246_A 5582.740 -5.970226 0.01800813
rs56389811_T 6726.716 -6.428983 0.21584288
rs11168250_T 6744.640 -6.421046 0.19956945
rs11168252_A 6708.584 -6.412005 0.18871140
rm(res, gwas, X, fit)

Region 3

Very strong signals for COA, very week signals for AOA.

1. All asthma cases

region = "chr17_33500001_39800000"
res = readRDS(paste0("/Users/nicholeyang/Downloads/survivalsusie/result/asthma_self_report/result/all/fit.susie.", region, ".rds"))
gwas = readRDS(paste0("/Users/nicholeyang/downloads/survivalsusie/result/gwas_surv/all_gwas_", region, ".rds"))
fit = res[[1]]
X = res[[2]]
print(res[[3]])
     user    system   elapsed 
268671.80 153861.85  22346.75 
pip <- logisticsusie:::get_pip(fit$alpha)
effect_estimate <- colSums(fit$alpha * fit$mu)
pip.sorted = sort(pip, decreasing = TRUE)
pip.sorted[1:10]
 [1] 0.82458455 0.70062370 0.30037338 0.12825964 0.11330354 0.10882188
 [7] 0.10600954 0.09703998 0.09403512 0.08914059
class(fit) = "susie"
cs <- susie_get_cs(fit, X)
cs
$cs
$cs$L1
 [1] 1467 1470 1471 1478 1479 1484 1491 1493 1501 1524

$cs$L2
[1] 3086 3350


$purity
   min.abs.corr mean.abs.corr median.abs.corr
L1    0.9748638     0.9939885       0.9991312
L2    0.8652420     0.8652420       0.8652420

$cs_index
[1] 1 2

$coverage
[1] 0.9607429 0.9994304

$requested_coverage
[1] 0.95
par(mfrow = c(1,2))
snps1 = colnames(X)[cs$cs$L1]
colors <- ifelse(rownames(gwas) %in% snps1, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.5, pch = 20, main = "CS 1")

snps2 = colnames(X)[cs$cs$L2]
colors <- ifelse(rownames(gwas) %in% snps2, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.5, pch = 20, main = "CS 2")

Version Author Date
fc8e7c4 yunqiyang0215 2024-06-20
cbind(gwas[rownames(gwas) %in% snps1, ], pip[sort(cs$cs$L1)])
                   MAF missing.rate  p.value.spa p.value.norm      Stat
rs11651596_C 0.4711579            0 9.947331e-36 9.618739e-36 -1570.315
rs12949100_A 0.4709323            0 7.665040e-36 7.409160e-36 -1572.964
rs8069176_A  0.4712393            0 1.199672e-35 1.160289e-35 -1568.498
rs4795399_C  0.4712166            0 8.096519e-36 7.827701e-36 -1573.030
rs2305480_A  0.4712096            0 9.136234e-36 8.833860e-36 -1571.808
rs11078926_A 0.4711895            0 9.466887e-36 9.153782e-36 -1571.440
rs11078927_T 0.4710100            0 1.059990e-35 1.024952e-35 -1570.176
rs12939832_A 0.4710054            0 8.304129e-36 8.027765e-36 -1572.232
rs4795400_T  0.4712255            0 7.079374e-36 6.843430e-36 -1573.472
rs9303279_C  0.4799832            0 1.659336e-35 1.610852e-35 -1565.204
                  Var         z           
rs11651596_C 15832.73 -12.47983 0.08914059
rs12949100_A 15833.46 -12.50060 0.11330354
rs8069176_A  15834.00 -12.46489 0.07470562
rs4795399_C  15845.87 -12.49623 0.10882188
rs2305480_A  15845.64 -12.48661 0.09703998
rs11078926_A 15845.40 -12.48378 0.09403512
rs11078927_T 15842.77 -12.47478 0.08500385
rs12939832_A 15834.88 -12.49422 0.10600954
rs4795400_T  15827.70 -12.50691 0.12825964
rs9303279_C  15834.02 -12.43871 0.07128375
cbind(gwas[rownames(gwas) %in% snps2, ], pip[sort(cs$cs$L2)])
                     MAF missing.rate  p.value.spa p.value.norm     Stat
rs112401631_A 0.02299226            0 3.594625e-10 2.742929e-10 227.7786
rs8067124_T   0.02203134            0 1.540705e-09 1.236998e-09 200.3148
                   Var        z          
rs112401631_A 1301.973 6.312653 0.7006237
rs8067124_T   1087.127 6.075373 0.3003734
rm(res, gwas, X, fit)

2. COA

region = "chr17_33500001_39800000"
res = readRDS(paste0("/Users/nicholeyang/Downloads/survivalsusie/result/asthma_self_report/result/coa/fit.susie.", region, ".rds"))
gwas = readRDS(paste0("/Users/nicholeyang/downloads/survivalsusie/result/gwas_surv/coa_gwas_", region, ".rds"))
fit = res[[1]]
X = res[[2]]

pip <- logisticsusie:::get_pip(fit$alpha)
effect_estimate <- colSums(fit$alpha * fit$mu)
pip.sorted = sort(pip, decreasing = TRUE)

class(fit) = "susie"
cs <- susie_get_cs(fit, X)
cs
$cs
$cs$L1
[1] 1467 1470 1471 1478 1479 1484 1491 1493


$purity
   min.abs.corr mean.abs.corr median.abs.corr
L1    0.9987126      0.999371       0.9993555

$cs_index
[1] 1

$coverage
[1] 0.9730042

$requested_coverage
[1] 0.95
snps1 = colnames(X)[cs$cs$L1]
colors <- ifelse(rownames(gwas) %in% snps1, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.5, pch = 20, main = "CS 1")

Version Author Date
fc8e7c4 yunqiyang0215 2024-06-20
cbind(gwas[rownames(gwas) %in% snps1, ], pip[sort(cs$cs$L1)])
                   MAF missing.rate  p.value.spa p.value.norm      Stat
rs11651596_C 0.4716738            0 1.020407e-85 1.761956e-86 -1244.913
rs12949100_A 0.4714520            0 2.488787e-85 4.344104e-86 -1242.038
rs8069176_A  0.4717552            0 1.071309e-85 1.855596e-86 -1244.792
rs4795399_C  0.4717367            0 7.864877e-86 1.353189e-86 -1246.273
rs2305480_A  0.4717295            0 8.198479e-86 1.411140e-86 -1246.129
rs11078926_A 0.4717091            0 1.010142e-85 1.744655e-86 -1245.442
rs11078927_T 0.4715354            0 3.534607e-85 6.223003e-86 -1241.244
rs12939832_A 0.4715304            0 3.507188e-85 6.174890e-86 -1240.954
                  Var         z           
rs11651596_C 3989.285 -19.71022 0.14754075
rs12949100_A 3989.365 -19.66450 0.06671948
rs8069176_A  3989.573 -19.70760 0.14008476
rs4795399_C  3992.592 -19.72357 0.19316881
rs2305480_A  3992.531 -19.72145 0.18545768
rs11078926_A 3992.474 -19.71072 0.15131288
rs11078927_T 3991.667 -19.64626 0.04795115
rs12939832_A 3989.644 -19.64666 0.04836222
rm(res, gwas, X, fit)

3. AOA

region = "chr17_33500001_39800000"
res = readRDS(paste0("/Users/nicholeyang/Downloads/survivalsusie/result/asthma_self_report/result/aoa/fit.susie.", region, ".rds"))
gwas = readRDS(paste0("/Users/nicholeyang/downloads/survivalsusie/result/gwas_surv/aoa_gwas_", region, ".rds"))
fit = res[[1]]
X = res[[2]]

pip <- logisticsusie:::get_pip(fit$alpha)
effect_estimate <- colSums(fit$alpha * fit$mu)
pip.sorted = sort(pip, decreasing = TRUE)

class(fit) = "susie"
cs <- susie_get_cs(fit, X)
cs
$cs
NULL

$coverage
NULL

$requested_coverage
[1] 0.95
rm(res, gwas, X, fit)

Region 4

Marginal significant signals for both COA and AOA. Combined analysis a lot more significant.

1. All asthma cases

region = "chr10_6600001_12200000"
res = readRDS(paste0("/Users/nicholeyang/Downloads/survivalsusie/result/asthma_self_report/result/all/fit.susie.", region, ".rds"))
gwas = readRDS(paste0("/Users/nicholeyang/downloads/survivalsusie/result/gwas_surv/all_gwas_", region, ".rds"))
fit = res[[1]]
X = res[[2]]
print(res[[3]])
     user    system   elapsed 
488558.23 162450.30  72205.62 
pip <- logisticsusie:::get_pip(fit$alpha)
effect_estimate <- colSums(fit$alpha * fit$mu)
pip.sorted = sort(pip, decreasing = TRUE)
pip.sorted[1:10]
 [1] 0.8606128 0.7544745 0.6370307 0.5253848 0.3574363 0.2456690 0.2340421
 [8] 0.1934556 0.1697711 0.1198051
class(fit) = "susie"
cs <- susie_get_cs(fit, X)
cs
$cs
$cs$L1
 [1] 2465 2467 2472 2473 2475 2477 2482 2483 2487 2501

$cs$L3
[1] 2435 2440 2453

$cs$L5
[1] 2365 2391

$cs$L7
[1] 2433 3018

$cs$L6
 [1] 112 113 115 116 120 121 123 124 126 127 128 130 131 132 133 134 136 137 138
[20] 139 140 142 143

$cs$L2
 [1] 1531 1533 1536 1572 1573 1587 1590 1603 1616 1621 1628 1631 1633 1636 1640
[16] 1642 1645 1649 1655 1660 1663 1664 1665 1683 1684 1695 1715 1719 1727 1780


$purity
   min.abs.corr mean.abs.corr median.abs.corr
L1    0.9962847     0.9986513       0.9989900
L3    0.9795058     0.9862823       0.9796464
L5    0.9698670     0.9698670       0.9698670
L7    0.8873465     0.8873465       0.8873465
L6    0.7206002     0.9193265       0.9841399
L2    0.6739457     0.9555920       0.9948012

$cs_index
[1] 1 3 5 7 6 2

$coverage
[1] 0.9723606 0.9998067 0.9502191 0.9583051 0.9997691 0.9518624

$requested_coverage
[1] 0.95
par(mfrow = c(3,2))
snps1 = colnames(X)[cs$cs$L1]
colors <- ifelse(rownames(gwas) %in% snps1, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.5, pch = 20, main = "CS 1")

snps2 = colnames(X)[cs$cs$L2]
colors <- ifelse(rownames(gwas) %in% snps2, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.5, pch = 20, main = "CS 2")

snps3 = colnames(X)[cs$cs$L3]
colors <- ifelse(rownames(gwas) %in% snps3, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.5, pch = 20, main = "CS 3")

snps5 = colnames(X)[cs$cs$L5]
colors <- ifelse(rownames(gwas) %in% snps5, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.5, pch = 20, main = "CS 5")


snps6 = colnames(X)[cs$cs$L6]
colors <- ifelse(rownames(gwas) %in% snps6, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.5, pch = 20, main = "CS 6")

snps7 = colnames(X)[cs$cs$L7]
colors <- ifelse(rownames(gwas) %in% snps7, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.5, pch = 20, main = "CS 7")

cbind(gwas[rownames(gwas) %in% snps1, ], pip[sort(cs$cs$L1)])
                  MAF missing.rate  p.value.spa p.value.norm      Stat      Var
rs962992_C  0.4233676            0 4.972442e-36 4.602638e-36 -1561.201 15503.66
rs962993_T  0.4233362            0 4.244034e-36 3.926629e-36 -1562.794 15504.17
rs1775553_T 0.4231244            0 5.990803e-36 5.546556e-36 -1559.401 15504.49
rs1775554_C 0.4232281            0 7.134297e-36 6.609131e-36 -1557.634 15503.81
rs1775555_C 0.4234788            0 7.750566e-36 7.183546e-36 -1557.319 15513.96
rs1663687_A 0.4231007            0 8.284016e-36 7.675699e-36 -1556.507 15510.86
rs1663680_C 0.4234260            0 8.028930e-36 7.441096e-36 -1557.934 15533.18
rs1031163_T 0.4231330            0 7.225009e-36 6.691733e-36 -1558.794 15529.36
rs1444782_A 0.4233832            0 1.026801e-35 9.522526e-36 -1555.862 15540.63
rs2197415_G 0.4226570            0 2.467360e-36 2.277027e-36 -1566.524 15471.91
                    z           
rs962992_C  -12.53839 0.09949565
rs962993_T  -12.55097 0.11098298
rs1775553_T -12.52359 0.08463862
rs1775554_C -12.50968 0.07540126
rs1775555_C -12.50306 0.07407266
rs1663687_A -12.49779 0.07108740
rs1663680_C -12.50026 0.07564261
rs1031163_T -12.50869 0.07989495
rs1444782_A -12.48063 0.06846180
rs2197415_G -12.59404 0.23404207
cbind(gwas[rownames(gwas) %in% snps2, ], pip[sort(cs$cs$L2)])
                   MAF missing.rate  p.value.spa p.value.norm      Stat
rs11255890_A 0.3814107            0 3.749630e-06 3.747194e-06  563.2263
rs11255891_C 0.3811078            0 3.538737e-06 3.536402e-06  564.8737
rs7923068_A  0.3763112            0 6.214948e-06 6.211267e-06  549.0645
rs7087891_T  0.3818541            0 6.440784e-06 6.437136e-06  551.5516
rs7079263_T  0.3818331            0 7.987898e-06 7.983613e-06  546.0822
rs7090156_A  0.3819658            0 8.062806e-06 8.058492e-06  546.2313
rs10795674_T 0.3818882            0 7.872046e-06 7.867806e-06  546.8361
rs4749820_A  0.3819934            0 7.467303e-06 7.463229e-06  548.4964
rs10795677_C 0.3820472            0 7.206742e-06 7.202777e-06  549.6670
rs3928823_A  0.3816941            0 5.804970e-06 5.801590e-06  555.2804
rs11255912_A 0.3818426            0 5.072942e-06 5.069891e-06  558.8195
rs9665552_C  0.3821648            0 6.404849e-06 6.401220e-06  553.0131
rs9665567_T  0.3820250            0 6.477949e-06 6.474288e-06  552.7773
rs11255914_G 0.3820245            0 6.537510e-06 6.533823e-06  552.5458
rs10795678_T 0.3821332            0 6.370477e-06 6.366862e-06  553.2126
rs2027105_T  0.3820411            0 6.512476e-06 6.508800e-06  552.7926
rs10905464_A 0.3820435            0 6.505651e-06 6.501978e-06  552.8488
rs1556593_A  0.3820064            0 6.505568e-06 6.501895e-06  552.9813
rs7902526_G  0.3821339            0 7.090815e-06 7.086899e-06  550.8873
rs2986300_T  0.3821073            0 7.131334e-06 7.127400e-06  550.7870
rs169693_T   0.3818606            0 7.205248e-06 7.201278e-06  549.9620
rs401305_G   0.3820712            0 6.958835e-06 6.954971e-06  551.3828
rs290356_G   0.3820713            0 7.061460e-06 7.057554e-06  550.9938
rs2483937_G  0.3816637            0 7.836184e-06 7.831953e-06  547.9434
rs2483936_A  0.3816479            0 7.570474e-06 7.566350e-06  548.8293
rs2244336_C  0.3816293            0 7.777904e-06 7.773695e-06  547.9782
rs7068268_T  0.3826083            0 8.222822e-06 8.218454e-06  545.9362
rs12769121_G 0.3825829            0 7.481953e-06 7.477879e-06  548.2823
rs2025758_C  0.4570126            0 5.477560e-14 5.455417e-14 -941.3337
rs11255938_T 0.4549757            0 8.176240e-14 8.143698e-14 -934.2983
                  Var         z            
rs11255890_A 14830.47  4.624933 0.004672291
rs11255891_C 14840.33  4.636919 0.004893552
rs7923068_A  14762.13  4.519067 0.002786895
rs7087891_T  14946.20  4.511498 0.002822863
rs7079263_T  14953.81  4.465623 0.002268905
rs7090156_A  14975.39  4.463624 0.002266996
rs10795674_T 14974.15  4.468750 0.002303589
rs4749820_A  14989.45  4.480031 0.002472998
rs10795677_C 15002.71  4.487607 0.002536395
rs3928823_A  15002.33  4.533494 0.003021933
rs11255912_A 15005.70  4.561875 0.003468309
rs9665552_C  15017.61  4.512685 0.002831367
rs9665567_T  15020.83  4.510278 0.002844658
rs11255914_G 15021.18  4.508335 0.002815488
rs10795678_T 15020.85  4.513825 0.002890039
rs2027105_T  15029.17  4.509150 0.002827371
rs10905464_A 15030.75  4.509372 0.002830226
rs1556593_A  15037.94  4.509375 0.002812388
rs7902526_G  15046.22  4.491062 0.002608556
rs2986300_T  15048.88  4.489848 0.002283002
rs169693_T   15018.52  4.487651 0.002319646
rs401305_G   15046.49  4.495062 0.002340193
rs290356_G   15046.12  4.491945 0.002306369
rs2483937_G  15028.28  4.469727 0.002215504
rs2483936_A  15027.30  4.477101 0.002294813
rs2244336_C  15019.46  4.471324 0.002235058
rs7068268_T  14987.47  4.459413 0.002217577
rs12769121_G 14980.54  4.479612 0.002437201
rs2025758_C  15667.16 -7.520531 0.525384811
rs11255938_T 15651.83 -7.467977 0.357436300
cbind(gwas[rownames(gwas) %in% snps3, ], pip[sort(cs$cs$L3)])
                    MAF missing.rate  p.value.spa p.value.norm      Stat
rs186856025_T 0.1076229            0 6.037658e-25 2.694409e-25 -809.3564
rs144536148_G 0.1076069            0 7.225399e-25 3.244956e-25 -807.9486
rs12413578_T  0.1040860            0 4.269344e-26 1.669032e-26 -819.0107
                   Var         z          
rs186856025_T 6065.624 -10.39207 0.1934556
rs144536148_G 6065.237 -10.37433 0.1697711
rs12413578_T  5909.462 -10.65407 0.6370307
cbind(gwas[rownames(gwas) %in% snps5, ], pip[sort(cs$cs$L5)])
                    MAF missing.rate  p.value.spa p.value.norm      Stat
rs72782675_T 0.01139900            0 6.101269e-11 3.132679e-11 -174.1911
rs11256010_C 0.01180942            0 1.709305e-10 9.657797e-11 -172.0632
                  Var         z          
rs72782675_T 688.1613 -6.640195 0.7544745
rs11256010_C 706.7571 -6.472213 0.2456690
cbind(gwas[rownames(gwas) %in% snps6, ], pip[sort(cs$cs$L6)])
                   MAF missing.rate  p.value.spa p.value.norm     Stat      Var
rs10905360_T 0.4582294            0 3.372720e-08 3.369323e-08 693.0571 15757.63
rs1361152_T  0.4555337            0 8.386657e-08 8.379190e-08 672.4307 15745.82
rs10905361_T 0.4555329            0 8.692650e-08 8.684954e-08 671.9820 15762.92
rs11255685_G 0.3091178            0 4.595135e-07 4.586948e-07 587.9150 13591.90
rs11255686_C 0.4561533            0 7.598315e-08 7.591467e-08 675.0231 15762.51
rs11255687_T 0.3094110            0 5.264154e-07 5.255075e-07 584.6758 13582.64
rs2388715_A  0.4032956            0 5.973238e-08 5.965609e-08 669.8155 15273.52
rs1572597_G  0.4048370            0 1.092954e-07 1.091706e-07 656.0356 15259.55
rs10905362_C 0.4035023            0 7.419131e-08 7.410009e-08 665.3488 15289.15
rs10905363_A 0.4010923            0 1.555422e-07 1.553714e-07 647.9115 15253.23
rs10905364_G 0.4035100            0 7.532316e-08 7.523077e-08 665.0269 15289.85
rs10905365_G 0.4035108            0 7.591432e-08 7.582135e-08 664.8256 15288.60
rs7082651_C  0.4035309            0 7.047604e-08 7.038864e-08 666.4066 15285.25
rs7082798_C  0.4035165            0 7.282078e-08 7.273095e-08 665.6286 15282.94
rs7082816_A  0.4035296            0 7.438961e-08 7.429818e-08 665.1285 15281.76
rs7082946_C  0.3951242            0 6.739997e-08 6.731047e-08 663.2095 15093.94
rs1338057_A  0.4035502            0 8.414118e-08 8.403993e-08 661.7569 15252.95
rs1338058_C  0.4035186            0 7.989791e-08 7.980089e-08 662.8972 15252.31
rs6602301_T  0.4033038            0 8.010373e-08 8.000636e-08 662.5799 15240.34
rs7084475_G  0.4031563            0 9.161644e-08 9.150748e-08 659.4124 15232.52
rs7088182_A  0.4032129            0 7.052657e-08 7.043891e-08 665.1473 15228.27
rs7100949_T  0.4033142            0 8.016838e-08 8.007094e-08 662.3217 15229.30
rs7100961_T  0.3110061            0 5.035850e-07 5.027192e-07 585.0123 13552.23
                    z           
rs10905360_T 5.521075 0.04873610
rs1361152_T  5.358767 0.01936750
rs10905361_T 5.352287 0.01861149
rs11255685_G 5.042831 0.01173543
rs11255686_C 5.376578 0.02099712
rs11255687_T 5.016758 0.01037761
rs2388715_A  5.419829 0.06572784
rs1572597_G  5.310757 0.03829371
rs10905362_C 5.380935 0.05444661
rs10905363_A 5.246078 0.02494528
rs10905364_G 5.378208 0.05371247
rs10905365_G 5.376800 0.05328616
rs7082651_C  5.390176 0.05716896
rs7082798_C  5.384291 0.05530267
rs7082816_A  5.380454 0.05430836
rs7082946_C  5.398206 0.06285292
rs1338057_A  5.358233 0.04768706
rs1338058_C  5.367578 0.05006357
rs6602301_T  5.367114 0.05064724
rs7084475_G  5.342828 0.04489357
rs7088182_A  5.390048 0.05739413
rs7100949_T  5.366969 0.05050511
rs7100961_T  5.025272 0.01129927
cbind(gwas[rownames(gwas) %in% snps7, ], pip[sort(cs$cs$L7)])
                    MAF missing.rate  p.value.spa p.value.norm     Stat
rs11256016_A 0.05104932            0 7.631291e-19 3.347436e-19 487.8163
rs17406680_C 0.05256578            0 2.022166e-16 1.177382e-16 466.0854
                  Var        z          
rs11256016_A 2966.354 8.956633 0.8606128
rs17406680_C 3164.512 8.285374 0.1198051
rm(res, gwas, X, fit)

2. COA

region = "chr10_6600001_12200000"
res = readRDS(paste0("/Users/nicholeyang/Downloads/survivalsusie/result/asthma_self_report/result/coa/fit.susie.", region, ".rds"))
gwas = readRDS(paste0("/Users/nicholeyang/downloads/survivalsusie/result/gwas_surv/coa_gwas_", region, ".rds"))
fit = res[[1]]
X = res[[2]]

pip <- logisticsusie:::get_pip(fit$alpha)
effect_estimate <- colSums(fit$alpha * fit$mu)
pip.sorted = sort(pip, decreasing = TRUE)

class(fit) = "susie"
cs <- susie_get_cs(fit, X)
cs
$cs
$cs$L3
[1] 2435 2440 2453

$cs$L1
 [1] 2441 2451 2456 2463 2464 2468 2492 2495 2496 2499 2503 2505 2506 2508 2510
[16] 2511


$purity
   min.abs.corr mean.abs.corr median.abs.corr
L3    0.9795657     0.9863213       0.9797085
L1    0.7377178     0.9345094       0.9879257

$cs_index
[1] 3 1

$coverage
[1] 0.9503942 0.9510038

$requested_coverage
[1] 0.95
par(mfrow = c(1,2))
snps1 = colnames(X)[cs$cs$L1]
colors <- ifelse(rownames(gwas) %in% snps1, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.5, pch = 20, main = "CS 1")
snps3 = colnames(X)[cs$cs$L3]
colors <- ifelse(rownames(gwas) %in% snps3, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.5, pch = 20, main = "CS 3")

cbind(gwas[rownames(gwas) %in% snps1, ], pip[sort(cs$cs$L1)])
                   MAF missing.rate  p.value.spa p.value.norm     Stat      Var
rs11256017_C 0.1837118            0 1.053318e-13 7.690263e-14 365.5747 2391.498
rs2589561_G  0.1872850            0 7.054922e-14 5.138640e-14 371.0499 2429.211
rs2440781_G  0.1873717            0 1.418919e-13 1.055768e-13 366.9615 2436.842
rs1775550_A  0.1873137            0 1.091430e-13 8.057847e-14 368.8975 2439.177
rs1775551_A  0.1873124            0 1.101068e-13 8.131135e-14 368.8524 2439.358
rs2797288_A  0.1873133            0 1.099283e-13 8.117563e-14 368.8745 2439.507
rs957349_A   0.1881942            0 5.070604e-14 3.664190e-14 375.2046 2455.107
rs2589559_T  0.1883260            0 9.656862e-14 7.120593e-14 371.2176 2459.239
rs6602349_C  0.1878566            0 1.483481e-13 1.106448e-13 367.8782 2453.123
rs2589563_T  0.1878786            0 1.626690e-13 1.216697e-13 367.1926 2452.287
rs725861_G   0.1881653            0 5.614297e-14 4.070136e-14 373.7336 2444.704
rs1444788_C  0.1879988            0 4.179163e-14 3.000021e-14 375.3407 2440.149
rs1444789_C  0.1870411            0 5.549340e-14 4.008932e-14 371.7285 2417.281
rs1612986_C  0.1861778            0 2.059662e-13 1.543271e-13 361.6793 2399.592
rs1342773_A  0.2601550            0 8.829267e-10 8.375701e-10 339.6538 3062.460
rs1663693_C  0.2586433            0 2.103092e-09 2.007103e-09 331.5144 3055.648
                    z           
rs11256017_C 7.475514 0.01866303
rs2589561_G  7.528348 0.04360345
rs2440781_G  7.433731 0.02546739
rs1775550_A  7.469372 0.03134852
rs1775551_A  7.468180 0.03114277
rs2797288_A  7.468400 0.03118268
rs957349_A   7.572390 0.05617388
rs2589559_T  7.485628 0.03346991
rs6602349_C  7.427529 0.02422673
rs2589563_T  7.414951 0.02246583
rs725861_G   7.558733 0.05118060
rs1444788_C  7.598319 0.06458213
rs1444789_C  7.560704 0.05158561
rs1612986_C  7.383375 0.02152232
rs1342773_A  6.137639 0.31484788
rs1663693_C  5.997231 0.14023622
cbind(gwas[rownames(gwas) %in% snps3, ], pip[sort(cs$cs$L3)])
                    MAF missing.rate  p.value.spa p.value.norm      Stat
rs186856025_T 0.1078996            0 1.348827e-10 1.012989e-10 -253.0634
rs144536148_G 0.1078832            0 1.453655e-10 1.094788e-10 -252.5957
rs12413578_T  0.1043563            0 3.721099e-11 2.614579e-11 -257.5814
                   Var         z          
rs186856025_T 1532.222 -6.465000 0.2596517
rs144536148_G 1532.129 -6.453247 0.2416930
rs12413578_T  1492.778 -6.666790 0.5000388
rm(res, gwas, X, fit)

3. AOA

region = "chr10_6600001_12200000"
res = readRDS(paste0("/Users/nicholeyang/Downloads/survivalsusie/result/asthma_self_report/result/aoa/fit.susie.", region, ".rds"))
gwas = readRDS(paste0("/Users/nicholeyang/downloads/survivalsusie/result/gwas_surv/aoa_gwas_", region, ".rds"))
fit = res[[1]]
X = res[[2]]

pip <- logisticsusie:::get_pip(fit$alpha)
effect_estimate <- colSums(fit$alpha * fit$mu)
pip.sorted = sort(pip, decreasing = TRUE)

class(fit) = "susie"
cs <- susie_get_cs(fit, X)
cs
$cs
$cs$L1
 [1] 2434 2465 2467 2472 2473 2475 2477 2482 2483 2487 2501


$purity
   min.abs.corr mean.abs.corr median.abs.corr
L1    0.9793508     0.9954759        0.998709

$cs_index
[1] 1

$coverage
[1] 0.9999378

$requested_coverage
[1] 0.95
snps1 = colnames(X)[cs$cs$L1]
colors <- ifelse(rownames(gwas) %in% snps1, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.5, pch = 20, main = "CS 1")

cbind(gwas[rownames(gwas) %in% snps1, ], pip[sort(cs$cs$L1)])
                   MAF missing.rate  p.value.spa p.value.norm      Stat
rs61840192_A 0.4326324            0 1.430427e-17 1.398486e-17 -807.2539
rs962992_C   0.4247181            0 3.369056e-17 3.292047e-17 -802.2441
rs962993_T   0.4246893            0 3.309988e-17 3.234219e-17 -802.4573
rs1775553_T  0.4244705            0 3.322487e-17 3.246321e-17 -802.4365
rs1775554_C  0.4245711            0 2.879613e-17 2.813033e-17 -804.0078
rs1775555_C  0.4248292            0 4.624008e-17 4.520674e-17 -798.9812
rs1663687_A  0.4244477            0 5.136474e-17 5.022077e-17 -797.7244
rs1663680_C  0.4247732            0 3.977349e-17 3.887450e-17 -801.1533
rs1031163_T  0.4244853            0 5.527005e-17 5.404475e-17 -797.3789
rs1444782_A  0.4247341            0 5.303978e-17 5.186418e-17 -798.1299
rs2197415_G  0.4240186            0 3.115001e-17 3.042797e-17 -802.3379
                  Var         z           
rs61840192_A 8945.147 -8.535253 0.15296259
rs962992_C   9044.196 -8.435708 0.09475255
rs962993_T   9044.562 -8.437780 0.09638994
rs1775553_T  9045.029 -8.437343 0.09528988
rs1775554_C  9044.573 -8.454077 0.10970193
rs1775555_C  9050.354 -8.398539 0.06943893
rs1663687_A  9048.524 -8.386177 0.06336783
rs1663680_C  9061.384 -8.416245 0.08166521
rs1031163_T  9059.334 -8.377541 0.05926589
rs1444782_A  9065.912 -8.382389 0.06152557
rs2197415_G  9026.608 -8.444910 0.12070713
rm(res, gwas, X, fit)

Region 5

COA: pval = 1e-40, AOA no significant signals.

1. All asthma cases

region = "chr1_150600001_155100000"
res = readRDS(paste0("/Users/nicholeyang/Downloads/survivalsusie/result/asthma_self_report/result/all/fit.susie.", region, ".rds"))
gwas = readRDS(paste0("/Users/nicholeyang/downloads/survivalsusie/result/gwas_surv/all_gwas_", region, ".rds"))
fit = res[[1]]
X = res[[2]]
print(res[[3]])
     user    system   elapsed 
184006.62  58608.10  26540.03 
pip <- logisticsusie:::get_pip(fit$alpha)
effect_estimate <- colSums(fit$alpha * fit$mu)
pip.sorted = sort(pip, decreasing = TRUE)
pip.sorted[1:10]
 [1] 0.99999963 0.99948127 0.07271724 0.02383702 0.01798914 0.01246631
 [7] 0.01168878 0.01120367 0.01083357 0.01072128
class(fit) = "susie"
cs <- susie_get_cs(fit, X)
cs
$cs
$cs$L1
[1] 1501

$cs$L2
[1] 1951


$purity
   min.abs.corr mean.abs.corr median.abs.corr
L1            1             1               1
L2            1             1               1

$cs_index
[1] 1 2

$coverage
[1] 0.9994803 0.9999996

$requested_coverage
[1] 0.95
par(mfrow = c(1,2))
snps1 = colnames(X)[cs$cs$L1]
colors <- ifelse(rownames(gwas) %in% snps1, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.5, pch = 20, main = "CS 1")

snps2 = colnames(X)[cs$cs$L2]
colors <- ifelse(rownames(gwas) %in% snps2, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.5, pch = 20, main = "CS 2")

cbind(gwas[rownames(gwas) %in% snps1, ], pip[sort(cs$cs$L1)])
                     [,1]      [,2]
MAF          4.788278e-02 0.9994813
missing.rate 0.000000e+00 0.9994813
p.value.spa  4.503327e-19 0.9994813
p.value.norm 1.822049e-19 0.9994813
Stat         4.864129e+02 0.9994813
Var          2.905775e+03 0.9994813
z            9.023481e+00 0.9994813
cbind(gwas[rownames(gwas) %in% snps2, ], pip[sort(cs$cs$L2)])
                     [,1]      [,2]
MAF          2.309783e-02 0.9999996
missing.rate 0.000000e+00 0.9999996
p.value.spa  8.920365e-16 0.9999996
p.value.norm 3.297129e-16 0.9999996
Stat         2.918389e+02 0.9999996
Var          1.278500e+03 0.9999996
z            8.161930e+00 0.9999996
rm(res, gwas, X, fit)

2. COA

region = "chr1_150600001_155100000"
res = readRDS(paste0("/Users/nicholeyang/Downloads/survivalsusie/result/asthma_self_report/result/coa/fit.susie.", region, ".rds"))
gwas = readRDS(paste0("/Users/nicholeyang/downloads/survivalsusie/result/gwas_surv/coa_gwas_", region, ".rds"))
fit = res[[1]]
X = res[[2]]

pip <- logisticsusie:::get_pip(fit$alpha)
effect_estimate <- colSums(fit$alpha * fit$mu)
pip.sorted = sort(pip, decreasing = TRUE)

class(fit) = "susie"
cs <- susie_get_cs(fit, X)
cs
$cs
$cs$L1
[1] 1951

$cs$L2
[1] 1501


$purity
   min.abs.corr mean.abs.corr median.abs.corr
L1            1             1               1
L2            1             1               1

$cs_index
[1] 1 2

$coverage
[1] 1.0000000 0.9751952

$requested_coverage
[1] 0.95
par(mfrow = c(1,2))
snps1 = colnames(X)[cs$cs$L1]
colors <- ifelse(rownames(gwas) %in% snps1, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.5, pch = 20, main = "CS 1")
snps2 = colnames(X)[cs$cs$L2]
colors <- ifelse(rownames(gwas) %in% snps2, "red", "black")
plot(-log10(gwas[, "p.value.spa"]), col = colors, xlab = "SNP", ylab = "-log10(p-value)", cex = 0.5, pch = 20, main = "CS 2")

cbind(gwas[rownames(gwas) %in% snps1, ], pip[sort(cs$cs$L1)])
                     [,1] [,2]
MAF          2.310334e-02    1
missing.rate 0.000000e+00    1
p.value.spa  3.660474e-46    1
p.value.norm 1.173288e-55    1
Stat         2.820831e+02    1
Var          3.221547e+02    1
z            1.571610e+01    1
cbind(gwas[rownames(gwas) %in% snps2, ], pip[sort(cs$cs$L2)])
                     [,1]      [,2]
MAF          4.777014e-02 0.9752416
missing.rate 0.000000e+00 0.9752416
p.value.spa  1.618538e-33 0.9752416
p.value.norm 2.895329e-37 0.9752416
Stat         3.448891e+02 0.9752416
Var          7.310485e+02 0.9752416
z            1.275576e+01 0.9752416
rm(res, gwas, X, fit)

3. AOA

region = "chr1_150600001_155100000"
res = readRDS(paste0("/Users/nicholeyang/Downloads/survivalsusie/result/asthma_self_report/result/aoa/fit.susie.", region, ".rds"))
gwas = readRDS(paste0("/Users/nicholeyang/downloads/survivalsusie/result/gwas_surv/aoa_gwas_", region, ".rds"))
fit = res[[1]]
X = res[[2]]

pip <- logisticsusie:::get_pip(fit$alpha)
effect_estimate <- colSums(fit$alpha * fit$mu)
pip.sorted = sort(pip, decreasing = TRUE)

class(fit) = "susie"
cs <- susie_get_cs(fit, X)
cs
$cs
NULL

$coverage
NULL

$requested_coverage
[1] 0.95

sessionInfo()
R version 4.1.1 (2021-08-10)
Platform: x86_64-apple-darwin20.6.0 (64-bit)
Running under: macOS Monterey 12.0.1

Matrix products: default
BLAS:   /usr/local/Cellar/openblas/0.3.18/lib/libopenblasp-r0.3.18.dylib
LAPACK: /usr/local/Cellar/r/4.1.1_1/lib/R/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] logisticsusie_0.0.0.9004 testthat_3.1.0           susieR_0.12.35          
[4] survival_3.2-11          workflowr_1.6.2         

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.8.3       lattice_0.20-44    prettyunits_1.1.1  ps_1.6.0          
 [5] rprojroot_2.0.2    digest_0.6.28      utf8_1.2.2         R6_2.5.1          
 [9] plyr_1.8.6         RcppZiggurat_0.1.6 evaluate_0.14      highr_0.9         
[13] ggplot2_3.4.3      pillar_1.9.0       rlang_1.1.1        rstudioapi_0.13   
[17] irlba_2.3.5        whisker_0.4        callr_3.7.3        jquerylib_0.1.4   
[21] Matrix_1.5-3       rmarkdown_2.11     desc_1.4.0         devtools_2.4.2    
[25] splines_4.1.1      stringr_1.4.0      munsell_0.5.0      mixsqp_0.3-43     
[29] compiler_4.1.1     httpuv_1.6.3       xfun_0.27          pkgconfig_2.0.3   
[33] pkgbuild_1.2.0     htmltools_0.5.5    tidyselect_1.2.0   tibble_3.1.5      
[37] matrixStats_0.63.0 reshape_0.8.9      fansi_0.5.0        crayon_1.4.1      
[41] dplyr_1.0.7        withr_2.5.0        later_1.3.0        grid_4.1.1        
[45] jsonlite_1.7.2     gtable_0.3.0       lifecycle_1.0.3    git2r_0.28.0      
[49] magrittr_2.0.1     scales_1.2.1       Rfast_2.0.6        cli_3.6.1         
[53] stringi_1.7.5      cachem_1.0.6       fs_1.5.0           promises_1.2.0.1  
[57] remotes_2.4.2      bslib_0.4.1        ellipsis_0.3.2     generics_0.1.2    
[61] vctrs_0.6.3        tools_4.1.1        glue_1.4.2         purrr_0.3.4       
[65] parallel_4.1.1     processx_3.8.1     pkgload_1.2.3      fastmap_1.1.0     
[69] yaml_2.2.1         colorspace_2.0-2   sessioninfo_1.1.1  memoise_2.0.1     
[73] knitr_1.36         usethis_2.1.3      sass_0.4.4