Miami plot¶
In [ ]:
Copied!
!wget -O bmi_male_bbj.txt.gz http://jenger.riken.jp/2analysisresult_qtl_download/
!wget -O bmi_female_bbj.txt.gz http://jenger.riken.jp/4analysisresult_qtl_download/
!wget -O bmi_male_bbj.txt.gz http://jenger.riken.jp/2analysisresult_qtl_download/
!wget -O bmi_female_bbj.txt.gz http://jenger.riken.jp/4analysisresult_qtl_download/
Load gwaslab and sample sumstats¶
In [1]:
Copied!
import sys
sys.path.insert(0,"/home/yunye/work/gwaslab/src")
import gwaslab as gl
import sys
sys.path.insert(0,"/home/yunye/work/gwaslab/src")
import gwaslab as gl
In [2]:
Copied!
gl1 = gl.Sumstats("bbj_bmi_female.txt.gz",fmt="gwaslab",snpid="SNP",ea="REF",nea="ALT",sep="\t",build="19",verbose=False)
gl2 = gl.Sumstats("bbj_bmi_male.txt.gz",fmt="gwaslab",snpid="SNP",ea="REF",nea="ALT",sep="\t",build="19",verbose=False)
gl1 = gl.Sumstats("bbj_bmi_female.txt.gz",fmt="gwaslab",snpid="SNP",ea="REF",nea="ALT",sep="\t",build="19",verbose=False)
gl2 = gl.Sumstats("bbj_bmi_male.txt.gz",fmt="gwaslab",snpid="SNP",ea="REF",nea="ALT",sep="\t",build="19",verbose=False)
check lead variants in sumstast1¶
In [3]:
Copied!
gl1.get_lead()
gl1.get_lead()
Fri Feb 2 18:01:59 2024 Start to extract lead variants...v3.4.38 Fri Feb 2 18:01:59 2024 -Current Dataframe shape : 5961600 x 9 ; Memory usage: 326.95 MB Fri Feb 2 18:01:59 2024 -Processing 5961600 variants... Fri Feb 2 18:01:59 2024 -Significance threshold : 5e-08 Fri Feb 2 18:01:59 2024 -Sliding window size: 500 kb Fri Feb 2 18:02:02 2024 -Using P for extracting lead variants... Fri Feb 2 18:02:02 2024 -Found 948 significant variants in total... Fri Feb 2 18:02:02 2024 -Identified 20 lead variants! Fri Feb 2 18:02:02 2024 Finished extracting lead variants.
Out[3]:
| SNPID | CHR | POS | EA | NEA | BETA | SE | P | STATUS | |
|---|---|---|---|---|---|---|---|---|---|
| 5629059 | rs860295 | 1 | 155767708 | A | G | 0.04005 | 0.006973 | 9.221000e-09 | 1999999 |
| 3574382 | rs532504 | 1 | 177878933 | G | A | 0.06024 | 0.006294 | 1.062000e-21 | 1999999 |
| 5726357 | rs939584 | 2 | 621558 | C | T | 0.05446 | 0.008828 | 6.897000e-10 | 1999999 |
| 4443132 | rs713586 | 2 | 25158008 | T | C | 0.03073 | 0.005345 | 9.015000e-09 | 1999999 |
| 2166247 | rs1846974 | 5 | 87969927 | G | A | 0.02988 | 0.005307 | 1.808000e-08 | 1999999 |
| 2629989 | rs261966 | 5 | 95849587 | T | C | 0.03216 | 0.005361 | 1.991000e-09 | 1999999 |
| 32091 | rs10061263 | 5 | 124289158 | G | C | 0.02870 | 0.005244 | 4.443000e-08 | 1999999 |
| 3386429 | rs4712523 | 6 | 20657564 | A | G | -0.03863 | 0.005296 | 3.016000e-13 | 1999999 |
| 3108709 | rs3798519 | 6 | 50788778 | A | C | 0.03638 | 0.005857 | 5.218000e-10 | 1999999 |
| 958518 | rs11981973 | 7 | 69445341 | A | G | 0.03594 | 0.006569 | 4.460000e-08 | 1999999 |
| 265445 | rs10811658 | 9 | 22128600 | G | A | 0.03117 | 0.005283 | 3.613000e-09 | 1999999 |
| 2436673 | rs2237897 | 11 | 2858546 | C | T | 0.03657 | 0.005451 | 1.953000e-11 | 1999999 |
| 1777942 | rs1491850 | 11 | 27749725 | T | C | -0.04251 | 0.005376 | 2.635000e-15 | 1999999 |
| 5422141 | rs7933262 | 11 | 28666538 | C | T | -0.03033 | 0.005457 | 2.727000e-08 | 1999999 |
| 689872 | rs11602339 | 11 | 47761471 | C | T | 0.03414 | 0.005790 | 3.728000e-09 | 1999999 |
| 4952305 | rs75766425 | 14 | 52511911 | G | C | 0.05405 | 0.008145 | 3.222000e-11 | 1999999 |
| 718615 | rs11642015 | 16 | 53802494 | C | T | 0.07912 | 0.006551 | 1.375000e-33 | 1999999 |
| 5590983 | rs8098510 | 18 | 40773435 | A | C | -0.03211 | 0.005647 | 1.296000e-08 | 1999999 |
| 4126646 | rs6567160 | 18 | 57829135 | T | C | 0.05420 | 0.006344 | 1.299000e-17 | 1999999 |
| 3024448 | rs35560038 | 19 | 46175046 | A | T | -0.05872 | 0.005473 | 7.421000e-27 | 1999999 |
Create Miami plot by iteratively calling plot_mqq()¶
In [5]:
Copied!
a = gl.plot_miami2(path1= gl1,
path2= gl2)
a = gl.plot_miami2(path1= gl1,
path2= gl2)
Fri Feb 2 18:02:48 2024 Start to create miami plot v3.4.38: Fri Feb 2 18:02:48 2024 -Obtaining Sumstats1 CHR, POS, P and annotation from: ['CHR', 'POS', 'P'] Fri Feb 2 18:02:48 2024 -Loading Sumstats1 from gwaslab.Sumstats Object Fri Feb 2 18:02:48 2024 -Obtaining Sumstats2 CHR, POS, P and annotation from: ['CHR', 'POS', 'P'] Fri Feb 2 18:02:48 2024 -Loading Sumstats2 from gwaslab.Sumstats Object Fri Feb 2 18:02:49 2024 -Sanity check after conversion: 0 variants with P value outside of (0,1] will be removed... Fri Feb 2 18:02:49 2024 -Sumstats P values are being converted to -log10(P)... Fri Feb 2 18:02:50 2024 -Sanity check: 0 na/inf/-inf variants will be removed... Fri Feb 2 18:02:51 2024 -Sanity check after conversion: 0 variants with P value outside of (0,1] will be removed... Fri Feb 2 18:02:51 2024 -Sumstats P values are being converted to -log10(P)... Fri Feb 2 18:02:51 2024 -Sanity check: 0 na/inf/-inf variants will be removed... Fri Feb 2 18:02:52 2024 -Merging sumstats using chr and pos... Fri Feb 2 18:02:59 2024 -Columns in merged sumstats: P_1,scaled_P_1,TCHR+POS,P_2,scaled_P_2,CHR,POS,_ADD,i Fri Feb 2 18:02:59 2024 Start to create Manhattan plot for sumstats1... Fri Feb 2 18:02:59 2024 Start to create MQQ plot...v3.4.38: Fri Feb 2 18:02:59 2024 -Genomic coordinates version: 19... Fri Feb 2 18:02:59 2024 -Genome-wide significance level to plot is set to 5e-08 ... Fri Feb 2 18:02:59 2024 -Raw input contains 5961600 variants... Fri Feb 2 18:02:59 2024 -MQQ plot layout mode is : m Fri Feb 2 18:02:59 2024 Finished loading specified columns from the sumstats. Fri Feb 2 18:02:59 2024 Start data conversion and sanity check: Fri Feb 2 18:02:59 2024 -Sanity check will be skipped. Fri Feb 2 18:02:59 2024 -Sanity check after conversion: 0 variants with P value outside of (0,1] will be removed... Fri Feb 2 18:02:59 2024 -Sumstats P values are being converted to -log10(P)... Fri Feb 2 18:03:00 2024 -Sanity check: 0 na/inf/-inf variants will be removed... Fri Feb 2 18:03:00 2024 -Converting data above cut line... Fri Feb 2 18:03:00 2024 -Maximum -log10(P) value is 32.86169730183372 . Fri Feb 2 18:03:00 2024 Finished data conversion and sanity check. Fri Feb 2 18:03:00 2024 Start to create MQQ plot with 5961600 variants... Fri Feb 2 18:03:01 2024 -Creating background plot... Fri Feb 2 18:03:11 2024 Finished creating MQQ plot successfully! Fri Feb 2 18:03:11 2024 Start to extract variants for annotation... Fri Feb 2 18:03:11 2024 -Found 20 significant variants with a sliding window size of 500 kb... Fri Feb 2 18:03:11 2024 Finished extracting variants for annotation... Fri Feb 2 18:03:11 2024 Start to process figure arts. Fri Feb 2 18:03:11 2024 -Processing X ticks... Fri Feb 2 18:03:11 2024 -Processing X labels... Fri Feb 2 18:03:11 2024 -Processing Y labels... Fri Feb 2 18:03:11 2024 -Processing Y tick lables... Fri Feb 2 18:03:11 2024 -Processing Y labels... Fri Feb 2 18:03:11 2024 -Processing lines... Fri Feb 2 18:03:11 2024 Finished processing figure arts. Fri Feb 2 18:03:11 2024 Start to annotate variants... Fri Feb 2 18:03:11 2024 -Skip annotating Fri Feb 2 18:03:11 2024 Finished annotating variants. Fri Feb 2 18:03:11 2024 Start to save figure... Fri Feb 2 18:03:11 2024 -Skip saving figure! Fri Feb 2 18:03:11 2024 Finished saving figure... Fri Feb 2 18:03:11 2024 Finished creating plot successfully! Fri Feb 2 18:03:11 2024 Finished creating Manhattan plot for sumstats1 Fri Feb 2 18:03:11 2024 Start to create Manhattan plot for sumstats2... Fri Feb 2 18:03:11 2024 Start to create MQQ plot...v3.4.38: Fri Feb 2 18:03:11 2024 -Genomic coordinates version: 19... Fri Feb 2 18:03:11 2024 -Genome-wide significance level to plot is set to 5e-08 ... Fri Feb 2 18:03:11 2024 -Raw input contains 5961600 variants... Fri Feb 2 18:03:11 2024 -MQQ plot layout mode is : m Fri Feb 2 18:03:11 2024 Finished loading specified columns from the sumstats. Fri Feb 2 18:03:11 2024 Start data conversion and sanity check: Fri Feb 2 18:03:11 2024 -Sanity check will be skipped. Fri Feb 2 18:03:11 2024 -Sanity check after conversion: 0 variants with P value outside of (0,1] will be removed... Fri Feb 2 18:03:11 2024 -Sumstats P values are being converted to -log10(P)... Fri Feb 2 18:03:12 2024 -Sanity check: 0 na/inf/-inf variants will be removed... Fri Feb 2 18:03:12 2024 -Converting data above cut line... Fri Feb 2 18:03:12 2024 -Maximum -log10(P) value is 40.59842715432355 . Fri Feb 2 18:03:12 2024 Finished data conversion and sanity check. Fri Feb 2 18:03:12 2024 Start to create MQQ plot with 5961600 variants... Fri Feb 2 18:03:13 2024 -Creating background plot... Fri Feb 2 18:03:23 2024 Finished creating MQQ plot successfully! Fri Feb 2 18:03:23 2024 Start to extract variants for annotation... Fri Feb 2 18:03:23 2024 -Found 28 significant variants with a sliding window size of 500 kb... Fri Feb 2 18:03:23 2024 Finished extracting variants for annotation... Fri Feb 2 18:03:23 2024 Start to process figure arts. Fri Feb 2 18:03:23 2024 -Processing X ticks... Fri Feb 2 18:03:23 2024 -Processing X labels... Fri Feb 2 18:03:23 2024 -Processing Y labels... Fri Feb 2 18:03:23 2024 -Processing Y tick lables... Fri Feb 2 18:03:23 2024 -Processing Y labels... Fri Feb 2 18:03:23 2024 -Processing lines... Fri Feb 2 18:03:23 2024 Finished processing figure arts. Fri Feb 2 18:03:23 2024 Start to annotate variants... Fri Feb 2 18:03:23 2024 -Skip annotating Fri Feb 2 18:03:23 2024 Finished annotating variants. Fri Feb 2 18:03:23 2024 Start to save figure... Fri Feb 2 18:03:23 2024 -Skip saving figure! Fri Feb 2 18:03:23 2024 Finished saving figure... Fri Feb 2 18:03:23 2024 Finished creating plot successfully! Fri Feb 2 18:03:23 2024 Finished creating Manhattan plot for sumstats2 Fri Feb 2 18:03:23 2024 Start to save figure... Fri Feb 2 18:03:23 2024 -Skip saving figure! Fri Feb 2 18:03:23 2024 Finished saving figure... Fri Feb 2 18:03:23 2024 Finished creating miami plot successfully
Create customized Miami plot¶
Most options in plot_mqq() are available for plot_miami2()
Simply add 1 or 2 after the option for plot_mqq() to customize the top or bottom figure in miami plot. If no prefix, the parameter will be passed to both plots.
In [7]:
Copied!
a = gl.plot_miami2(path1= gl1,
path2= gl2,
skip=2,
cut1=20,
cut2=15,
id1="SNPID",
id2="SNPID",
anno1=True,
anno2="GENENAME",
additional_line1=[1e-14],
anno_set1=["rs3798519"],
pinpoint1=[["rs3798519","rs35560038"],["rs7933262","rs8098510"]],
pinpoint_color1=["purple","black"],
highlight1=["rs8098510"],
highlight2=[["rs8098510","rs3798519"], ["rs1491850"]],
highlight_color2=["red","yellow"],
jagged=True,
verbose1=False,
verbose2=False
)
a = gl.plot_miami2(path1= gl1,
path2= gl2,
skip=2,
cut1=20,
cut2=15,
id1="SNPID",
id2="SNPID",
anno1=True,
anno2="GENENAME",
additional_line1=[1e-14],
anno_set1=["rs3798519"],
pinpoint1=[["rs3798519","rs35560038"],["rs7933262","rs8098510"]],
pinpoint_color1=["purple","black"],
highlight1=["rs8098510"],
highlight2=[["rs8098510","rs3798519"], ["rs1491850"]],
highlight_color2=["red","yellow"],
jagged=True,
verbose1=False,
verbose2=False
)
Fri Feb 2 18:06:01 2024 Start to create miami plot v3.4.38: Fri Feb 2 18:06:01 2024 -Obtaining Sumstats1 CHR, POS, P and annotation from: ['CHR', 'POS', 'P', 'SNPID'] Fri Feb 2 18:06:01 2024 -Loading Sumstats1 from gwaslab.Sumstats Object Fri Feb 2 18:06:01 2024 -Obtaining Sumstats2 CHR, POS, P and annotation from: ['CHR', 'POS', 'P', 'SNPID'] Fri Feb 2 18:06:01 2024 -Loading Sumstats2 from gwaslab.Sumstats Object Fri Feb 2 18:06:02 2024 -Sanity check after conversion: 0 variants with P value outside of (0,1] will be removed... Fri Feb 2 18:06:03 2024 -Sumstats P values are being converted to -log10(P)... Fri Feb 2 18:06:03 2024 -Sanity check: 0 na/inf/-inf variants will be removed... Fri Feb 2 18:06:04 2024 -Sanity check after conversion: 0 variants with P value outside of (0,1] will be removed... Fri Feb 2 18:06:04 2024 -Sumstats P values are being converted to -log10(P)... Fri Feb 2 18:06:05 2024 -Sanity check: 0 na/inf/-inf variants will be removed... Fri Feb 2 18:06:06 2024 -Merging sumstats using chr and pos... Fri Feb 2 18:06:15 2024 -Columns in merged sumstats: P_1,SNPID_1,scaled_P_1,TCHR+POS,P_2,SNPID_2,scaled_P_2,CHR,POS,_ADD,i Fri Feb 2 18:06:15 2024 Start to create Manhattan plot for sumstats1... Fri Feb 2 18:06:18 2024 -Processing jagged Y axis... Fri Feb 2 18:06:18 2024 Finished creating Manhattan plot for sumstats1 Fri Feb 2 18:06:18 2024 Start to create Manhattan plot for sumstats2... Fri Feb 2 18:06:23 2024 -Processing jagged Y axis... Fri Feb 2 18:06:23 2024 Finished creating Manhattan plot for sumstats2 Fri Feb 2 18:06:23 2024 Start to save figure... Fri Feb 2 18:06:23 2024 -Skip saving figure! Fri Feb 2 18:06:23 2024 Finished saving figure... Fri Feb 2 18:06:23 2024 Finished creating miami plot successfully
