Last updated: 2024-06-20
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 cb1583e. 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 | cb1583e | yunqiyang0215 | 2024-06-20 | wflow_publish("analysis/susie_asthma_result1.Rmd") |
Coxph Susie result on all asthma/ AOA/ COA in UKBiobank.
library(survival)
library(susieR)
devtools::load_all("/Users/nicholeyang/Downloads/logisticsusie")
ℹ Loading logisticsusie
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.
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")
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
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)
pip.sorted = sort(pip, decreasing = TRUE)
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")
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
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)
pip.sorted = sort(pip, decreasing = TRUE)
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")
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()
No GWAS significant signal for COA, marginal significant for AOA.
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
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)
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
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)
pip.sorted = sort(pip, decreasing = TRUE)
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")
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)
Very strong signals for COA, very week signals for AOA.
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")
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)
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")
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)
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)
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