Power and Trumpet plot¶
In [1]:
Copied!
import gwaslab as gl
import gwaslab as gl
Calculate power using .get_power¶
In [10]:
Copied!
gl.get_power(mode="b",
ncase=1000,
ncontrol=11000,
beta=0.3,
prevalence=0.1,
daf=0.2)
gl.get_power(mode="b",
ncase=1000,
ncontrol=11000,
beta=0.3,
prevalence=0.1,
daf=0.2)
Sun Jul 23 22:03:17 2023 Start to calculate statistical power... Sun Jul 23 22:03:17 2023 -Input settings (b mode): Sun Jul 23 22:03:17 2023 -Number of cases:1000 Sun Jul 23 22:03:17 2023 -Number of controls:11000 Sun Jul 23 22:03:17 2023 -Risk allele beta:0.300 Sun Jul 23 22:03:17 2023 -Disease prevalence:0.100 Sun Jul 23 22:03:17 2023 -Risk allele frequency: 0.200 Sun Jul 23 22:03:17 2023 -Significance level: 5.000e-08 Sun Jul 23 22:03:17 2023 -Alogorithm: Skol, Andrew D., et al. Nature genetics 38.2 (2006): 209-213.... Sun Jul 23 22:03:17 2023 -GRR is approximated using OR. For prevalence < 10%, GRR is very similar to OR.... Sun Jul 23 22:03:17 2023 Probability of disease : Sun Jul 23 22:03:17 2023 - Individuals with AA genotype: 0.149 Sun Jul 23 22:03:17 2023 - Individuals with AB genotype: 0.118 Sun Jul 23 22:03:17 2023 - Individuals with BB genotype: 0.088 Sun Jul 23 22:03:17 2023 Expected risk allele frequency: Sun Jul 23 22:03:17 2023 - In cases: 0.249 Sun Jul 23 22:03:17 2023 - In controls: 0.195 Sun Jul 23 22:03:17 2023 Expected power: 0.495 Sun Jul 23 22:03:17 2023 Finished calculating statistical power.
Out[10]:
0.4949100367998067
In [3]:
Copied!
gl.get_power(mode="q", eaf= 0.2, beta=0.3, n=1000)
gl.get_power(mode="q", eaf= 0.2, beta=0.3, n=1000)
Sat Feb 3 18:56:04 2024 Start to calculate statistical power... Sat Feb 3 18:56:04 2024 -Input settings (q mode): Sat Feb 3 18:56:04 2024 -Significance level: 5e-08 Sat Feb 3 18:56:04 2024 -EAF: 0.2 Sat Feb 3 18:56:04 2024 -BETA: 0.3 Sat Feb 3 18:56:04 2024 -N: 1000 Sat Feb 3 18:56:04 2024 -SNPR2: 0.028800000000000006 Sat Feb 3 18:56:04 2024 Finished calculating statistical power.
Out[3]:
0.4662311488597769
Plot power (mode="b")¶
X = MAF, colored by N_CONTROLS¶
In [4]:
Copied!
a = gl.plot_power_x(mode="b",
x="MAF",
ncases=1000,
ncontrols=[11000,15000,110000],
betas=0.3,
prevalences=0.1,
xscale="nonlog")
a = gl.plot_power_x(mode="b",
x="MAF",
ncases=1000,
ncontrols=[11000,15000,110000],
betas=0.3,
prevalences=0.1,
xscale="nonlog")
Sat Feb 3 18:56:26 2024 Start to create power plot... Sat Feb 3 18:56:26 2024 -Settings: Sat Feb 3 18:56:26 2024 -Mode: b Sat Feb 3 18:56:26 2024 -X axis: MAF Sat Feb 3 18:56:26 2024 -N_CASE: 1000 Sat Feb 3 18:56:26 2024 -N_CONTROL: [11000, 15000, 110000] Sat Feb 3 18:56:26 2024 -PREVALENCE: 0.1 Sat Feb 3 18:56:26 2024 -BETA: 0.3 Sat Feb 3 18:56:26 2024 -Significance level: 5e-08 Sat Feb 3 18:56:26 2024 Start to save figure... Sat Feb 3 18:56:26 2024 -Skip saving figure! Sat Feb 3 18:56:26 2024 Finished saving figure... Sat Feb 3 18:56:26 2024 Finished creating power plot!
X = N_CASES, colored by MAF¶
In [7]:
Copied!
a = gl.plot_power_x(mode="b",
mafs=[0.5,0.4,0.3,0.2,0.9,0.1],
betas= 0.1,
ncontrols=10000,
prevalences=0.1,
xscale="nonlog")
a = gl.plot_power_x(mode="b",
mafs=[0.5,0.4,0.3,0.2,0.9,0.1],
betas= 0.1,
ncontrols=10000,
prevalences=0.1,
xscale="nonlog")
Sat Feb 3 19:00:04 2024 Start to create power plot... Sat Feb 3 19:00:04 2024 -Settings: Sat Feb 3 19:00:04 2024 -Mode: b Sat Feb 3 19:00:04 2024 -X axis: None Sat Feb 3 19:00:04 2024 -N_CASE: 5000 Sat Feb 3 19:00:04 2024 -N_CONTROL: 10000 Sat Feb 3 19:00:04 2024 -PREVALENCE: 0.1 Sat Feb 3 19:00:04 2024 -BETA: 0.1 Sat Feb 3 19:00:04 2024 -Significance level: 5e-08 Sat Feb 3 19:00:04 2024 Start to save figure... Sat Feb 3 19:00:04 2024 -Skip saving figure! Sat Feb 3 19:00:04 2024 Finished saving figure... Sat Feb 3 19:00:04 2024 Finished creating power plot!
X = N_CASES, colored by BETA¶
In [8]:
Copied!
a = gl.plot_power_x(mode="b",
betas=[0.5,0.4,0.3,0.2],
mafs= 0.1,
xscale="nonlog")
a = gl.plot_power_x(mode="b",
betas=[0.5,0.4,0.3,0.2],
mafs= 0.1,
xscale="nonlog")
Sat Feb 3 19:00:29 2024 Start to create power plot... Sat Feb 3 19:00:29 2024 -Settings: Sat Feb 3 19:00:29 2024 -Mode: b Sat Feb 3 19:00:29 2024 -X axis: None Sat Feb 3 19:00:29 2024 -N_CASE: 5000 Sat Feb 3 19:00:29 2024 -N_CONTROL: 5000 Sat Feb 3 19:00:29 2024 -PREVALENCE: 0.1 Sat Feb 3 19:00:29 2024 -BETA: [0.5, 0.4, 0.3, 0.2] Sat Feb 3 19:00:29 2024 -Significance level: 5e-08 Sat Feb 3 19:00:29 2024 Start to save figure... Sat Feb 3 19:00:29 2024 -Skip saving figure! Sat Feb 3 19:00:29 2024 Finished saving figure... Sat Feb 3 19:00:29 2024 Finished creating power plot!
X = BETA, colored by N_CASES¶
In [20]:
Copied!
a = gl.plot_power_x(mode="b",x="BETA",ncases=[1000,5000,10000],xscale="nonlog")
a = gl.plot_power_x(mode="b",x="BETA",ncases=[1000,5000,10000],xscale="nonlog")
Sun Jul 23 22:14:45 2023 Start to create power plot... Sun Jul 23 22:14:45 2023 -Settings: Sun Jul 23 22:14:45 2023 -Mode: b Sun Jul 23 22:14:45 2023 -X axis: BETA Sun Jul 23 22:14:45 2023 -N_CASE: [1000, 5000, 10000] Sun Jul 23 22:14:45 2023 -N_CONTROL: 5000 Sun Jul 23 22:14:45 2023 -PREVALENCE: 0.1 Sun Jul 23 22:14:45 2023 -Significance level: 5e-08 Sun Jul 23 22:14:45 2023 -Skip saving figures! Sun Jul 23 22:14:45 2023 Finished creating power plot!
X = PREVALENCE, colored by N_CONTROLS¶
In [9]:
Copied!
a = gl.plot_power_x(mode="b",x="PREVALENCE",ncontrols=[100,500,1000],betas=0.1,xscale="nonlog")
a = gl.plot_power_x(mode="b",x="PREVALENCE",ncontrols=[100,500,1000],betas=0.1,xscale="nonlog")
Sat Feb 3 19:01:15 2024 Start to create power plot... Sat Feb 3 19:01:15 2024 -Settings: Sat Feb 3 19:01:15 2024 -Mode: b Sat Feb 3 19:01:15 2024 -X axis: PREVALENCE Sat Feb 3 19:01:15 2024 -N_CASE: 5000 Sat Feb 3 19:01:15 2024 -N_CONTROL: [100, 500, 1000] Sat Feb 3 19:01:15 2024 -BETA: 0.1 Sat Feb 3 19:01:15 2024 -Significance level: 5e-08 Sat Feb 3 19:01:15 2024 Start to save figure... Sat Feb 3 19:01:15 2024 -Skip saving figure! Sat Feb 3 19:01:15 2024 Finished saving figure... Sat Feb 3 19:01:15 2024 Finished creating power plot!
/home/yunye/work/gwaslab/src/gwaslab/util_in_calculate_power.py:86: RuntimeWarning: invalid value encountered in sqrt den= np.sqrt( (vcase/ncase + vcontrol/ncontrol)*0.5 )
Plot power (mode="q")¶
X = N, colored by MAF¶
In [12]:
Copied!
a = gl.plot_power_x( mafs=[0.5,0.4,0.3,0.2], betas= 0.1,xscale="nonlog")
a = gl.plot_power_x( mafs=[0.5,0.4,0.3,0.2], betas= 0.1,xscale="nonlog")
Sat Feb 3 19:03:45 2024 Start to create power plot... Sat Feb 3 19:03:45 2024 -Settings: Sat Feb 3 19:03:45 2024 -Mode: q Sat Feb 3 19:03:45 2024 -X axis: None Sat Feb 3 19:03:45 2024 -N: 10000 Sat Feb 3 19:03:45 2024 -MAF: [0.5, 0.4, 0.3, 0.2] Sat Feb 3 19:03:45 2024 -BETA: 0.1 Sat Feb 3 19:03:45 2024 -Significance level: 5e-08 Sat Feb 3 19:03:45 2024 Start to save figure... Sat Feb 3 19:03:45 2024 -Skip saving figure! Sat Feb 3 19:03:45 2024 Finished saving figure... Sat Feb 3 19:03:45 2024 Finished creating power plot!
X = N, colored by BETA¶
In [13]:
Copied!
a = gl.plot_power_x(betas=[0.5,0.4,0.3,0.2], mafs= 0.1,xscale="nonlog")
a = gl.plot_power_x(betas=[0.5,0.4,0.3,0.2], mafs= 0.1,xscale="nonlog")
Sat Feb 3 19:04:01 2024 Start to create power plot... Sat Feb 3 19:04:01 2024 -Settings: Sat Feb 3 19:04:01 2024 -Mode: q Sat Feb 3 19:04:01 2024 -X axis: None Sat Feb 3 19:04:01 2024 -N: 10000 Sat Feb 3 19:04:01 2024 -MAF: 0.1 Sat Feb 3 19:04:01 2024 -BETA: [0.5, 0.4, 0.3, 0.2] Sat Feb 3 19:04:01 2024 -Significance level: 5e-08 Sat Feb 3 19:04:01 2024 Start to save figure... Sat Feb 3 19:04:01 2024 -Skip saving figure! Sat Feb 3 19:04:01 2024 Finished saving figure... Sat Feb 3 19:04:01 2024 Finished creating power plot!
X = MAF, colored by N¶
In [14]:
Copied!
a = gl.plot_power_x(x="MAF",ns=[100,1000,10000],betas=1,xscale="nonlog")
a = gl.plot_power_x(x="MAF",ns=[100,1000,10000],betas=1,xscale="nonlog")
Sat Feb 3 19:04:17 2024 Start to create power plot... Sat Feb 3 19:04:17 2024 -Settings: Sat Feb 3 19:04:17 2024 -Mode: q Sat Feb 3 19:04:17 2024 -X axis: MAF Sat Feb 3 19:04:17 2024 -N: [100, 1000, 10000] Sat Feb 3 19:04:17 2024 -BETA: 1 Sat Feb 3 19:04:17 2024 -Significance level: 5e-08 Sat Feb 3 19:04:17 2024 Start to save figure... Sat Feb 3 19:04:17 2024 -Skip saving figure! Sat Feb 3 19:04:17 2024 Finished saving figure... Sat Feb 3 19:04:17 2024 Finished creating power plot!
X = BETA, colored by N¶
In [15]:
Copied!
a = gl.plot_power_x(x="BETA",ns=[100,1000,10000],xscale="nonlog")
a = gl.plot_power_x(x="BETA",ns=[100,1000,10000],xscale="nonlog")
Sat Feb 3 19:04:43 2024 Start to create power plot... Sat Feb 3 19:04:43 2024 -Settings: Sat Feb 3 19:04:43 2024 -Mode: q Sat Feb 3 19:04:43 2024 -X axis: BETA Sat Feb 3 19:04:43 2024 -N: [100, 1000, 10000] Sat Feb 3 19:04:43 2024 -MAF: 0.1 Sat Feb 3 19:04:43 2024 -Significance level: 5e-08 Sat Feb 3 19:04:43 2024 Start to save figure... Sat Feb 3 19:04:43 2024 -Skip saving figure! Sat Feb 3 19:04:43 2024 Finished saving figure... Sat Feb 3 19:04:43 2024 Finished creating power plot!
Plot power (Y=BETA, X=MAF)¶
In [16]:
Copied!
a = gl.plot_power()
a = gl.plot_power()
Sat Feb 3 19:05:19 2024 Start to create trumpet plot... Sat Feb 3 19:05:19 2024 -Updating eaf-beta matrix... Sat Feb 3 19:05:20 2024 -Extracting eaf-beta combinations with power = 0.2... Sat Feb 3 19:05:20 2024 -Updating eaf-beta matrix... Sat Feb 3 19:05:20 2024 -Extracting eaf-beta combinations with power = 0.4... Sat Feb 3 19:05:20 2024 -Updating eaf-beta matrix... Sat Feb 3 19:05:21 2024 -Extracting eaf-beta combinations with power = 0.6... Sat Feb 3 19:05:21 2024 -Updating eaf-beta matrix... Sat Feb 3 19:05:22 2024 -Extracting eaf-beta combinations with power = 0.8... Sat Feb 3 19:05:22 2024 Start to save figure... Sat Feb 3 19:05:22 2024 -Skip saving figure! Sat Feb 3 19:05:22 2024 Finished saving figure... Sat Feb 3 19:05:22 2024 Finished creating trumpet plot!
In [17]:
Copied!
a = gl.plot_power(mode="b",ncases=1000,ncontrols=1000, prevalences=0.1, ts=[0.1,0.2,0.3])
a = gl.plot_power(mode="b",ncases=1000,ncontrols=1000, prevalences=0.1, ts=[0.1,0.2,0.3])
Sat Feb 3 19:05:45 2024 Start to create trumpet plot... Sat Feb 3 19:05:45 2024 -Updating eaf-beta matrix... Sat Feb 3 19:05:45 2024 -GRR is approximated using OR. For prevalence < 10%, GRR is very similar to OR.... Sat Feb 3 19:05:46 2024 -Extracting eaf-beta combinations with power = 0.1... Sat Feb 3 19:05:46 2024 -Updating eaf-beta matrix... Sat Feb 3 19:05:46 2024 -GRR is approximated using OR. For prevalence < 10%, GRR is very similar to OR.... Sat Feb 3 19:05:46 2024 -Extracting eaf-beta combinations with power = 0.2... Sat Feb 3 19:05:46 2024 -Updating eaf-beta matrix... Sat Feb 3 19:05:46 2024 -GRR is approximated using OR. For prevalence < 10%, GRR is very similar to OR.... Sat Feb 3 19:05:46 2024 -Extracting eaf-beta combinations with power = 0.3... Sat Feb 3 19:05:46 2024 Start to save figure... Sat Feb 3 19:05:46 2024 -Skip saving figure! Sat Feb 3 19:05:46 2024 Finished saving figure... Sat Feb 3 19:05:46 2024 Finished creating trumpet plot!
Trumpet plot for quantitative traits¶
In [18]:
Copied!
#Sample data : Akiyama, M., Ishigaki, K., Sakaue, S., Momozawa, Y., Horikoshi, M., Hirata, M., ... & Kamatani, Y. (2019). Characterizing rare and low-frequency height-associated variants in the Japanese population. Nature communications, 10(1), 4393.
#Sample data : Akiyama, M., Ishigaki, K., Sakaue, S., Momozawa, Y., Horikoshi, M., Hirata, M., ... & Kamatani, Y. (2019). Characterizing rare and low-frequency height-associated variants in the Japanese population. Nature communications, 10(1), 4393.
In [21]:
Copied!
mysumstats = gl.Sumstats("bbj_height.tsv",fmt="auto",eaf="MAF",sep="\s+",other=["Gene"],n=160000,verbose=False)
mysumstats.basic_check(verbose=False)
mysumstats = gl.Sumstats("bbj_height.tsv",fmt="auto",eaf="MAF",sep="\s+",other=["Gene"],n=160000,verbose=False)
mysumstats.basic_check(verbose=False)
Sat Feb 3 19:06:22 2024 #WARNING! Necessary columns for .fix_allele() were not detected:EA,NEA Sat Feb 3 19:06:23 2024 #WARNING! Necessary columns for .normalize() were not detected:EA,NEA
In [22]:
Copied!
mysumstats.data
mysumstats.data
Out[22]:
| SNPID | CHR | POS | EAF | BETA | SE | P | N | STATUS | Gene | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | rs3829738 | 1 | 909309 | 0.209 | 0.023 | 0.004 | 1.800000e-08 | 160000 | 9980999 | PLEKHN1 |
| 1 | rs12748433 | 1 | 1559971 | 0.366 | -0.023 | 0.004 | 1.100000e-08 | 160000 | 9980999 | MIB2 |
| 2 | rs7551097 | 1 | 9259236 | 0.227 | -0.021 | 0.004 | 2.200000e-08 | 160000 | 9980999 | MIR34A,H6PD |
| 3 | rs118100927 | 1 | 9472981 | 0.102 | 0.037 | 0.005 | 3.500000e-12 | 160000 | 9980999 | SPSB1,LOC100506022 |
| 4 | rs3170740 | 1 | 17312743 | 0.318 | 0.041 | 0.003 | 1.500000e-33 | 160000 | 9980999 | ATP13A2 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 604 | rs140498234 | 23 | 105477754 | 0.025 | 0.068 | 0.010 | 3.960000e-11 | 160000 | 9980999 | MUM1L1,CXorf57 |
| 605 | rs3848873 | 23 | 118594422 | 0.416 | 0.024 | 0.003 | 4.540000e-17 | 160000 | 9980999 | SLC25A43,SLC25A5-AS1 |
| 606 | rs148300465 | 23 | 129325425 | 0.063 | 0.047 | 0.007 | 2.600000e-12 | 160000 | 9980999 | RAB33A,ZNF280C |
| 607 | rs7065171 | 23 | 133683590 | 0.217 | 0.027 | 0.004 | 1.880000e-13 | 160000 | 9980999 | LINC00629 |
| 608 | rs426349 | 23 | 139282091 | 0.497 | 0.023 | 0.003 | 8.120000e-12 | 160000 | 9980999 | LOC389895,SOX3 |
609 rows × 10 columns
In [25]:
Copied!
a = mysumstats.plot_trumpet(mode="q",
ts=[0.2,0.4,0.6,0.8] ,
anno="Gene",
anno_style="expand",
cmap="cool",
sig_level=5e-8,
build="19",
anno_x=0.01,
anno_y=1,
p_level=5e-8,
n_matrix=2000,
fontsize=12,
xscale="log",
repel_force=0.15,
yscale_factor=5.63,
sort="eaf",
ylim=(-6,4),
save=True)
a = mysumstats.plot_trumpet(mode="q",
ts=[0.2,0.4,0.6,0.8] ,
anno="Gene",
anno_style="expand",
cmap="cool",
sig_level=5e-8,
build="19",
anno_x=0.01,
anno_y=1,
p_level=5e-8,
n_matrix=2000,
fontsize=12,
xscale="log",
repel_force=0.15,
yscale_factor=5.63,
sort="eaf",
ylim=(-6,4),
save=True)
Sat Feb 3 19:07:23 2024 Start to create trumpet plot... Sat Feb 3 19:07:23 2024 -Settings: Sat Feb 3 19:07:23 2024 -Mode: q Sat Feb 3 19:07:23 2024 -N: N Sat Feb 3 19:07:23 2024 -BETA: BETA Sat Feb 3 19:07:23 2024 -Significance level: 5e-08 Sat Feb 3 19:07:23 2024 -Power thresholds: [0.2, 0.4, 0.6, 0.8] Sat Feb 3 19:07:23 2024 -Power line smoothness: 2000 Sat Feb 3 19:07:23 2024 -Loading column Gene for annotation... Sat Feb 3 19:07:23 2024 -Excluding variants with P values > 5e-08 Sat Feb 3 19:07:23 2024 -Plotting 609 variants... Sat Feb 3 19:07:23 2024 -Flipping 0 variants... Sat Feb 3 19:07:23 2024 -N for power calculation: 160000.0 Sat Feb 3 19:07:23 2024 -Updating eaf-beta matrix... Sat Feb 3 19:07:25 2024 -Extracting eaf-beta combinations with power = 0.2... Sat Feb 3 19:07:25 2024 -Updating eaf-beta matrix... Sat Feb 3 19:07:27 2024 -Extracting eaf-beta combinations with power = 0.4... Sat Feb 3 19:07:27 2024 -Updating eaf-beta matrix... Sat Feb 3 19:07:29 2024 -Extracting eaf-beta combinations with power = 0.6... Sat Feb 3 19:07:29 2024 -Updating eaf-beta matrix... Sat Feb 3 19:07:31 2024 -Extracting eaf-beta combinations with power = 0.8... Sat Feb 3 19:07:31 2024 Start to save figure... Sat Feb 3 19:07:32 2024 -Saved to ./gwaslab_trumpet_quant_20240203_2.png successfully! Sat Feb 3 19:07:32 2024 Finished saving figure... Sat Feb 3 19:07:32 2024 Finished creating trumpet plot!
In [27]:
Copied!
a = mysumstats.plot_trumpet(mode="q",
ts=[0.8] ,
anno="Gene",
anno_style="expand",
sig_level=5e-8,
build="19",
anno_x=0.01,
anno_y=1,
p_level=5e-8,
n_matrix=2000,
fontsize=12,
cmap="cool",
xscale="nonlog",
repel_force=0.05,
yscale_factor=5.63,
sort="beta",
ylim=(-6,4),
markercolor="#349beb")
a = mysumstats.plot_trumpet(mode="q",
ts=[0.8] ,
anno="Gene",
anno_style="expand",
sig_level=5e-8,
build="19",
anno_x=0.01,
anno_y=1,
p_level=5e-8,
n_matrix=2000,
fontsize=12,
cmap="cool",
xscale="nonlog",
repel_force=0.05,
yscale_factor=5.63,
sort="beta",
ylim=(-6,4),
markercolor="#349beb")
Sat Feb 3 19:08:12 2024 Start to create trumpet plot... Sat Feb 3 19:08:12 2024 -Settings: Sat Feb 3 19:08:12 2024 -Mode: q Sat Feb 3 19:08:12 2024 -N: N Sat Feb 3 19:08:12 2024 -BETA: BETA Sat Feb 3 19:08:12 2024 -Significance level: 5e-08 Sat Feb 3 19:08:12 2024 -Power thresholds: [0.8] Sat Feb 3 19:08:12 2024 -Power line smoothness: 2000 Sat Feb 3 19:08:12 2024 -Loading column Gene for annotation... Sat Feb 3 19:08:12 2024 -Excluding variants with P values > 5e-08 Sat Feb 3 19:08:12 2024 -Plotting 609 variants... Sat Feb 3 19:08:12 2024 -Flipping 0 variants... Sat Feb 3 19:08:12 2024 -N for power calculation: 160000.0 Sat Feb 3 19:08:12 2024 -Updating eaf-beta matrix... Sat Feb 3 19:08:14 2024 -Extracting eaf-beta combinations with power = 0.8... Sat Feb 3 19:08:14 2024 Start to save figure... Sat Feb 3 19:08:14 2024 -Skip saving figure! Sat Feb 3 19:08:14 2024 Finished saving figure... Sat Feb 3 19:08:14 2024 Finished creating trumpet plot!
Trumpet plot for binary traits¶
In [28]:
Copied!
mysumstats =gl.Sumstats("t2d_bbj.txt.gz",fmt="auto",verbose=False)
mysumstats.basic_check(verbose=False)
mysumstats =gl.Sumstats("t2d_bbj.txt.gz",fmt="auto",verbose=False)
mysumstats.basic_check(verbose=False)
In [29]:
Copied!
a = mysumstats.plot_trumpet(mode="b",
ncase=36614,
ncontrol=155150,
ts=[0.2,0.4,0.6,0.8] ,
anno="GENENAME",
anno_style="expand",
cmap="cool",
sig_level=5e-8,
build="19",
anno_x=0.01,
anno_y=0.2,
p_level=5e-8,
n_matrix=2000,
fontsize=12,
xscale="log",
repel_force=0.15,
sort="eaf",
ylim=(-3,3),
save=True)
a = mysumstats.plot_trumpet(mode="b",
ncase=36614,
ncontrol=155150,
ts=[0.2,0.4,0.6,0.8] ,
anno="GENENAME",
anno_style="expand",
cmap="cool",
sig_level=5e-8,
build="19",
anno_x=0.01,
anno_y=0.2,
p_level=5e-8,
n_matrix=2000,
fontsize=12,
xscale="log",
repel_force=0.15,
sort="eaf",
ylim=(-3,3),
save=True)
Sat Feb 3 19:12:40 2024 Start to create trumpet plot... Sat Feb 3 19:12:40 2024 -Prevalence is not given. Estimating based on scase and scontrol :0.19093260465989445... Sat Feb 3 19:12:40 2024 -Settings: Sat Feb 3 19:12:40 2024 -Mode: b Sat Feb 3 19:12:40 2024 -N_CASE: 36614 Sat Feb 3 19:12:40 2024 -N_CONTROL: 155150 Sat Feb 3 19:12:40 2024 -PREVALENCE: 0.19093260465989445 Sat Feb 3 19:12:40 2024 -BETA: BETA Sat Feb 3 19:12:40 2024 -Significance level: 5e-08 Sat Feb 3 19:12:40 2024 -Power thresholds: [0.2, 0.4, 0.6, 0.8] Sat Feb 3 19:12:40 2024 -Power line smoothness: 2000 Sat Feb 3 19:12:40 2024 -Excluding variants with P values > 5e-08 Sat Feb 3 19:12:40 2024 -Plotting 9461 variants... Sat Feb 3 19:12:40 2024 -Flipping 4518 variants... Sat Feb 3 19:12:40 2024 -Updating eaf-beta matrix... Sat Feb 3 19:12:40 2024 -GRR is approximated using OR. For prevalence < 10%, GRR is very similar to OR.... Sat Feb 3 19:12:41 2024 -Extracting eaf-beta combinations with power = 0.2... Sat Feb 3 19:12:41 2024 -Updating eaf-beta matrix... Sat Feb 3 19:12:41 2024 -GRR is approximated using OR. For prevalence < 10%, GRR is very similar to OR.... Sat Feb 3 19:12:42 2024 -Extracting eaf-beta combinations with power = 0.4... Sat Feb 3 19:12:42 2024 -Updating eaf-beta matrix... Sat Feb 3 19:12:42 2024 -GRR is approximated using OR. For prevalence < 10%, GRR is very similar to OR.... Sat Feb 3 19:12:43 2024 -Extracting eaf-beta combinations with power = 0.6... Sat Feb 3 19:12:43 2024 -Updating eaf-beta matrix... Sat Feb 3 19:12:43 2024 -GRR is approximated using OR. For prevalence < 10%, GRR is very similar to OR.... Sat Feb 3 19:12:44 2024 -Extracting eaf-beta combinations with power = 0.8... Sat Feb 3 19:12:44 2024 Start to annotate variants with nearest gene name(s)... Sat Feb 3 19:12:44 2024 -Assigning Gene name using ensembl_hg19_gtf for protein coding genes Sat Feb 3 19:12:45 2024 Finished annotating variants with nearest gene name(s) successfully! Sat Feb 3 19:12:45 2024 Start to save figure... Sat Feb 3 19:12:45 2024 -Saved to ./gwaslab_trumpet_binary_20240203_1.png successfully! Sat Feb 3 19:12:45 2024 Finished saving figure... Sat Feb 3 19:12:45 2024 Finished creating trumpet plot!
In [ ]:
Copied!