一行python画Manhattan plot与QQ plot

之前介绍过曼哈顿图与QQ图的画法,但自己画终究还是有点麻烦,有很多数据的时候就很头疼,于是自己写了一个非常简单的python package,一行代码画好对齐的Manhattan plot 与 QQ plot,并基于一个给定的滑动窗口自动检测top SNP并注释。

Package: gwaslab

安装方法: pip install gwaslab

使用方法:myplot = gl.mqqplot(insumstats,"CHR","POS","PVALUE")

当前版本:0.0.5

依赖的包:scipy,numpy,matplotlib,pandas,seaborn

(注意:非专业开发人员,完全用爱发电,自用的scripts改成的包,难免有bug,欢迎在github上指正,或提意见或建议, https://github.com/Cloufield/gwaslab )

使用示例:

实例数据下载: http://jenger.riken.jp/result

ID 1 Atrial Fibrillation

读入数据

输入文件是一般的gwas sumstats就可以,有chr,pos,pvalue三列就可以画:

import pandas as pd
import gwaslab as gl

#读取输入文件
insumstats = pd.read_csv("AF_1000G_rsq09_FINAL.txt.gz","\\s+")

本次使用的示例数据如下

使用mqqplot函数画图:

输入gwas sumstats的DataFrame,并指定染色体,碱基位置,p值的列名即可:

myplot = gl.mqqplot(insumstats,"CHR","POS","PVALUE")

#log
Basic settings:
  - Genome-wide significance level is set to 5e-08 ...
  - Raw input contains 5018048 variants...
Start conversion and QC:
  - P values are being converted to -log10(P)...
  - Sanity check: 0 variants with P value outside of (0,1] will be removed...
  - Sanity check: 0 na/inf/-inf variants will be removed...
  - Maximum -log10(P) values is 134.014573525917 .
Plotting 5018048 variants:
  - Found 14 significant variants with a sliding window size of 500 kb...
  - Skip annotating
  - Created Manhattan plot successfully!
  - Created QQ plot successfully!

lambda gc 自动计算并注释, qq plot与manhattan plot的y轴对应,,看着是舒服一点了,

可是突出的那个loci太高,其他的loci都快看不见了,有点不协调啊,显著的loci也没有注释,怎么办?

那么加几个option好了:

1.用cut 截断一下,cut以上的轴按cutfactor成比例缩小,这样就能在不丢失信息的情况下看清所有显著的loci

2.anno选项是用来自动注释lead snp的,默认基于一个500kb的滑动窗口,可以通过windowsize改变

gl.mqqplot(insumstats,"CHR","POS","PVALUE",cut=20,cutfactor=20,anno=True)

#log
Basic settings:
  - Genome-wide significance level is set to 5e-08 ...
  - Raw input contains 5018048 variants...
Start conversion and QC:
  - P values are being converted to -log10(P)...
  - Sanity check: 0 variants with P value outside of (0,1] will be removed...
  - Sanity check: 0 na/inf/-inf variants will be removed...
  - Maximum -log10(P) values is 134.014573525917 .
  - Minus log10(P) values above 20 will be shrunk with a shrinkage factor of 20...
Plotting 5018048 variants:
  - Found 14 significant variants with a sliding window size of 500 kb...
  - Annotating using column CHR:POS...
  - Created Manhattan plot successfully!
  - Created QQ plot successfully!

嗯嗯,看着舒服多了,

如果 只指定anno=True, 那么注释就只是chr: pos,当然也可以指定用于注释的列名,这里用SNP这一列注释:

gl.mqqplot(insumstats,"CHR","POS","PVALUE",cut=20,cutfactor=20,anno="SNP")

颜色不好看怎么办:

用colors选项,传个颜色的list进去,给图换个颜色,

qqplot的颜色,为了保持视觉统一,设定为曼哈顿图chr1的颜色,

例如传个彩虹进去:

myplot = gl.mqqplot(insumstats,"CHR","POS","PVALUE",
                    cut=20,cutfactor=20,anno=True,
                    colors=["#ff0000","#fc4444","#fc6404","#fcd444","#8cc43c","#029658","#1abc9c","#5bc0de","#6454ac","#fc8c84"])

其他可以选的 option

加title,调大小,什么的,目前还没有优化,后续会精简一下

mqqplot(insumstats,
          chrom,
          pos,
          p,
          scaled=False,
          cut=0,
          cutfactor=10,
          cut_line_color="#ebebeb",
          windowsizekb=500,
          anno=None,
          sig_level=5e-8,
          sig_line_color="grey",
          suggestive_sig_level=5e-6,
          title =None,
          mtitle=None,
          qtitle=None,
          figsize =(15,5),
          fontsize = 10,
          colors = ["#000042", "#7878BA"],
          layout="qqm",
          verbose=True,
          repel_force=0.03,
          gc=True,
          save=None,
          saveargs={"dpi":300,"facecolor":"white"}        
          )

保存的话 save=“./yourpath/myplot.png” ,或者save=True就行

分别画Manhattan plot 与 QQ plot

当然也可以分开画 mplot

mymplot = gl.mplot(insumstats,"CHR","POS","PVALUE",cut=20,cutfactor=10,anno=True)

qqplot

参考:

https://github.com/Cloufield/gwaslab

GWAS Gallery: gwaslab

使用corrplot绘制遗传相关矩阵

统计遗传学研究中经常会需要通过GWAS的sumstats来计算遗传相关,并以一下这种简单易懂的相关矩阵形式展示,本文就以两篇高质量文章为例,来介绍遗传相关矩阵的绘图方法。

回顾:GWASLab:通过Bivariate LD Score regression估计遗传相关性 genetic correlation

上图为An atlas of genetic correlations across human diseases and traits(文章链接:https://www.nature.com/articles/ng.3406)中的遗传相关矩阵的图,可以看到图中蓝色方块表示正向相关,红色方块表示负向相关,其中相关显著的表型(经过多重检验调整后)之间还会被 * 标记。这幅图中,方块的大小与颜色均表示相关系数,但这幅图还有个特别之处是FDR<0.01的方块被填满。(注意,这里的方块大小依然是相关系数,但近来的文章中,方块大小通常会被用来表示p值的大小,后续会介绍)

接下来我们尝试复现这幅图:

源数据下载:https://static-content.springer.com/esm/art%3A10.1038%2Fng.3406/MediaObjects/41588_2015_BFng3406_MOESM37_ESM.csv

源数据格式如下:

$ head 41588_2015_BFng3406_MOESM37_ESM.csv
Trait1,Trait2,rg,se,z,p
ADHD,Age at Menarche,-0.153,0.08218,-1.858,0.063
ADHD,Age at Smoking,-0.036,0.2427,-0.147,0.883
ADHD,Alzheimer's,-0.055,0.2191,-0.249,0.803
ADHD,Anorexia,0.192,0.1162,1.649,0.099
ADHD,Autism Spectrum,-0.164,0.1438,-1.144,0.253
ADHD,BMI,0.287,0.08913,3.222,0.001
ADHD,BMI 2010,0.258,0.08606,2.993,0.003
ADHD,Bipolar,0.265,0.1539,1.722,0.085
ADHD,Birth Length,-0.055,0.1679,-0.329,0.742

使用R的官方corrplot来绘制,示例代码如下:

install.packages("corrplot")
library(corrplot)
#读取数据
ldsc = read.csv(file = '41588_2015_BFng3406_MOESM37_ESM.csv',header=TRUE)

#提取需要画的表型
to_plot_col<-c('Age at Menarche','Alzheimer\'s','Anorexia','Autism Spectrum','BMI','Bipolar','Birth Length','Birth Weight','Childhood Obesity','Coronary Artery Disease','Crohn\'s Disease','Depression','Ever/Never Smoked','Fasting Glucose','HDL','Height','Infant Head Circumference','LDL','Rheumatoid Arthritis','Schizophrenia','T2D','Triglycerides','Ulcerative Colitis','Years of Education')

#构建一个rg矩阵,和一个p值矩阵,并命名row和col
gcm_rg<-matrix(1, nrow = length(to_plot_col), ncol = length(to_plot_col))
gcm_p<-matrix(1, nrow = length(to_plot_col), ncol = length(to_plot_col))
rownames(gcm_rg) <- to_plot_col
colnames(gcm_rg) <- to_plot_col
rownames(gcm_p) <-to_plot_col
colnames(gcm_p) <- to_plot_col

#将数据填进矩阵
for (row in 1:nrow(ldsc)) {
    if(ldsc[row, "Trait1"]%in%to_plot_col && ldsc[row, "Trait2"]%in%to_plot_col){
    Trait1 <- ldsc[row, "Trait1"]
    Trait2  <- ldsc[row, "Trait2"]
    gcm_rg[Trait1,Trait2]<-ldsc[row, "rg"]
    gcm_p[Trait1,Trait2]<-p<-ldsc[row, "p"]
    gcm_rg[Trait2,Trait1]<-ldsc[row, "rg"]
    gcm_p[Trait2,Trait1]<-p<-ldsc[row, "p"]
    }
    }

#使用corrplot绘图:(方法选择square)
#sig.level = 0.05/300 为显著水平(这里是多重检验调整的)
#gcm_rg为相关系数矩阵
#gcm_p为p值矩阵
#order="hclust" 为通过分层聚类对矩阵进行排序
#其他选项均为细节调整 可以参考 https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html 

corrplot(gcm_rg, p.mat=gcm_p, method="square",order="hclust",sig = "pch",sig.level = 0.05/300 ,pch="*",pch.cex=2,
         tl.col="black",tl.srt=45,addgrid.col="grey90")

对比与原图,已经接近90%的相似度了,不同之处在于这幅图还有个特别之处是FDR<0.01的方块被填满。这需要对底层代码进行修改。

这里我们就略过此步骤,因为目前主流方法是,方块大小应该表示p值大小, 以提高遗传相关矩阵的信息含量,而不是途中混合的表示方式。方块大小表示p值大小的示例,如下图所示:

(文章链接:https://www.nature.com/articles/s41588-018-0047-6

这张图中颜色表示相关系数,方块大小表示p值(fdr调整后)

这篇文章的作者提供了修改后的代码,可以通过github在R中安装修改后的corrplot:

github: https://github.com/mkanai/ldsc-corrplot-rg

install.packages("dplyr")
devtools::install_github("mkanai/corrplot")

安装后重新绘制:

#读取数据
ldsc = read.csv(file = '41588_2015_BFng3406_MOESM37_ESM.csv',header=TRUE)

#计算FDR调整后p值
ldsc$fdr<-p.adjust(ldsc$p, method ="fdr")

to_plot_col<-c('Age at Menarche','Alzheimer\'s','Anorexia','Autism Spectrum','BMI','Bipolar','Birth Length','Birth Weight','Childhood Obesity','Coronary Artery Disease','Crohn\'s Disease','Depression','Ever/Never Smoked','Fasting Glucose','HDL','Height','Infant Head Circumference','LDL','Rheumatoid Arthritis','Schizophrenia','T2D','Triglycerides','Ulcerative Colitis','Years of Education')

gcm_rg<-matrix(NA, nrow = length(to_plot_col), ncol = length(to_plot_col))
gcm_p<-matrix(NA, nrow = length(to_plot_col), ncol = length(to_plot_col))
rownames(gcm_rg) <- to_plot_col
colnames(gcm_rg) <- to_plot_col
rownames(gcm_p) <-to_plot_col
colnames(gcm_p) <- to_plot_col

#p值矩阵改为fdr矩阵
for (row in 1:nrow(ldsc)) {
    if(ldsc[row, "Trait1"]%in%to_plot_col && ldsc[row, "Trait2"]%in%to_plot_col){
    Trait1 <- ldsc[row, "Trait1"]
    Trait2  <- ldsc[row, "Trait2"]
    gcm_rg[Trait1,Trait2]<-ldsc[row, "rg"]
    gcm_p[Trait1,Trait2]<-ldsc[row, "fdr"]
    gcm_rg[Trait2,Trait1]<-ldsc[row, "rg"]
    gcm_p[Trait2,Trait1]<-ldsc[row, "fdr"]
    }
    }

#事先对矩阵进行排序(此修改的版本直接使用order="hclust"时有bug)
gcm_rg.hclust<-corrMatOrder(gcm_rg,order ="hclust")
gcm_rg_hclust<-gcm_rg[gcm_rg.hclust,gcm_rg.hclust]
gcm_p_hclust<-gcm_p[gcm_rg.hclust,gcm_rg.hclust]

options(repr.plot.width=15, repr.plot.height=15)
#使用修改后的corrplot绘图:(注意!!方法选择 psquare)
corrplot(gcm_rg_hclust ,p.mat=gcm_p_hclust, diag=FALSE,method="psquare",sig="pch",insig="label_sig",sig.level = 0.01 ,pch.cex=2,tl.col="black",tl.srt=45,addgrid.col="grey90")

修改后的图图例更为统一,对于p值大小的表现更为清楚。

对于图中的排序方法,字体大小,旋转,颜色, 色调,colorbar等等细节可以参考 https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html 进行调整。

参考:

https://www.nature.com/articles/ng.3406#MOESM37

https://www.rdocumentation.org/packages/corrplot/versions/0.92/topics/corrplot

https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html

https://www.nature.com/articles/s41

使用LocusZoom绘制Regional plot详解 – GWAS可视化

Locuszoom是一款非常好用的绘制region plot的软件,是GWAS研究中的必备软件之一。

Regional plot示例: https://my.locuszoom.org/gwas/881707/region/?chrom=3&start=45634967&end=46034967

Regional plot:横轴为染色体位置,左边纵轴为-log10(p),marker的颜色表示与leading SNP 的LD的大小,右侧纵轴则表示重组率。下方还标有gwas catalog中已知的与疾病相关的hit,以及该区域中基因的范围。该图主要用来用于精确定位因果SNP。

locuszoom官网: http://locuszoom.org/

该软件提供有最初的网页版(底层还是命令行版本),近来新发布的网页交互版(基于JavaScript),以及单独可下载的linux命令行版本(基于R/Python)。

网页版方便直接,命令行自定义程度更高。

本文以命令行版本为基础简要介绍locuszoom的使用方法,同时与网页版一一对照来方便理解与使用。

本文目录:

1. 命令行版本简介与软件与参考数据准备(**只看网页版的话可以跳过 1**)
2. 输入文件格式 
3. 指定绘制的区间范围
4. 指定LD文件,族裔群体,参考基因组版本
	**-》以上内容对应的网页版本**
5. 批量模式 ! Batch mode(命令行版本的优势所在)
	**-》批量模式对应的网页版本**
6. 自制LD文件
7 额外功能
	7.1指定多个参考SNP
	7.2 使用GWAS catalog进行注释
	7.3 精确定位的可信集
	7.4 对多个SNP的注释
	7.5 绘制额外的BED文件轨道
参考url

使用方法:

  1. 命令行版本简介 与 软件与参考数据准备 (只看网页版的话可以跳过 1

命令行版本总体上与网页版一一对应,但可自定义设置的程度更高,目前仅支持linux版本,windows和mac用户只能使用网页版或通过虚拟机。

下载地址:https://github.com/statgen/locuszoom-standalone

注意我们需要下载: 程序 + 数据库 + LD信息

(*单独下载软件无法绘制图像,或是需要使用自制的LD信息)

目前的最新版本为1.4版,发布于2017-05-01,共有23G。

除此之外,我们还需要准备以下有版本要求的python和R:

1.Python 2.7+ (注意要使用PYTHON2!)

2.R 3.0+.

locuzoom下载完成后,解压并添加进环境后即可使用:

cd <directory where you want to place locuszoom> 
tar zxf /path/to/locuszoom.tgz
ln -s bin/locuszoom /usr/local/bin/locuszoom #根据自己的环境变量制定

locuszoom的文件结构主要包括的内容如下所示:

locuszoom/
	bin/
		locuszoom (this is the locuszoom "executable")
		locuszoom.R (the R script which is used by locuszoom for creating the plots)
		dbmeister.py (script for creating custom user databases)
		lzupdate.py (script for creating an updated copy of the provided locuszoom database)
	conf/ (configuration file located here)
	data/
		database/ (SQLite file located here)
		hapmap/ (hapmap genotype files)
		1000G/ (1000G genotype files)
	src/ (source code for locuszoom)

命令行locuszoom语法

locuszoom [输入文件] [选项]

一个简单的例子:

locuszoom --metal Kathiresan_2009_HDL.txt --refgene FADS1

  1. 输入文件格式

locuszoom主要使用METAL 或者 EPACTS的文件格式,这里主要介绍METAL格式,因为其准备起来更加简单便捷。

METAL格式: 文件必须包含以tab分隔的下两列 1.markers 与 2.p-values

如下所

locuzoom在输入时的选项:
--metal 指定metal格式的输入文件名
--delim (可选)可以指定分隔符,默认为tab
--markercol  指定输入文件表示marker的列名,默认为"MarkerName"
--pvalcol 指定输入文件表示p值的列名,默认为"P-value"
--no-transform 当p已经转换为-log10(p)时,可以使用此选项跳过log转换

对应网页区块:

首先上传输入文件 , marker与P value的列名,以及指定分隔符

  1. 指定绘制的区间范围

在命令行版本中我们可以通过多种方式指定区间,包括

#2.1指定SNP与两侧的范围
--refsnp <your snp> --flank 500kb

#2.2指定SNP与区间的起止位点
--refsnp <your snp> --chr # --start <base position> --end <base position>

#2.3指定SNP与基因(绘制基因所在的范围)
--refsnp <your snp> --refgene <your gene>

#2.4 指定基因与两侧范围
--refgene <your gene> --flank 250kb

#2.5 指定基因与区间的起止位点
--refgene <your gene> --chr # --start <base position> --end <base position>

#2.6 区间的起止位点
--chr # --start <base position> --end <base position>


--flank 选项在这里指从起始和终止位点向外侧计算的距离,单位为kb, 在没有指定参考SNP时,locuzoom会自动选择最显著的SNP作为参考SNP。

网页版中指定绘制的区间范围所对应区块,填写方法与命令行版本一致:

4 指定LD文件,族裔群体,参考基因组版本

目前locuszoom内置了以下的组合,可以根据自己需要进行指定:

Genotype files available for: 
--source hapmap
  --build hg18
    --pop YRI
    --pop CEU
    --pop JPT+CHB

--source 1000G_March2012
  --build hg19
    --pop AMR
    --pop ASN
    --pop AFR
    --pop EUR

--source 1000G_June2010
  --build hg18
    --pop YRI
    --pop CEU
    --pop JPT+CHB

--source 1000G_Nov2014
  --build hg19
    --pop AMR
    --pop ASN
    --pop AFR
    --pop EUR
  --build hg38
    --pop SAS
    --pop EAS
    --pop AMR
    --pop AFR
    --pop EUR

一个简单的使用例:

--pop ASN --build hg19 --source 1000G_March2012
#指定1000G_March2012数据,以hg19为参考基因坐标,来计算ASN族裔的LD

以上内容对应的网页版:

Genome Build/LD Population中选择所使用的参考基因组版本与LD文件

有了以上信息就可以进行绘制。命令行中运行命令或是在网页版中点击Plot Data按钮。根据所绘制SNP数量多少,运行时间会有不同。等待即可。


5 批量模式 ! Batch mode(命令行版本的优势所在)

当我们有多个基因座,需要绘制多张regional plot,可以使用批量模式。

只需要通过 --hitspec 选项来指定一个包含多个区域信息的文件即可。

文件需要包含以下列:

  1. snp 2.chr 3.start 4.end 5.flank 6.run 7.m2zargs

该文件格式如下:

选项1-5与之单次模式的内容相同,run 是指是否指定该行进行绘制,m2zargs 则是对该行信息绘制时额外的参数。

对应的网页版区块与内容:

6 自制LD文件 (命令行版本)

可以通过--LD选项来指定自己的LD文件,文件格式如下:

dprime列可以为空,但Rsq列必须为有效的数据。该文件以空格分隔,且必须包含header。

7 额外功能 (未完待续,持续更新)

7.1指定多个参考SNP

7.2 使用GWAS catalog进行注释

7.3 精确定位的可信集

7.4 对多个SNP的注释

7.5 绘制额外的BED轨道

参考:

http://locuszoom.org/

https://genome.sph.umich.edu/wiki/LocusZoom_Standalone#User-supplied_LD

分位数-分位数图 QQ plot

QQ plot 是 GWAS结果可视化的非常有力的工具之一。 通过QQplot 我们可以直观地判断GWAS结果统计量是否存在inflation,是否存在过多的假阳性等,从而了解我们的分析是否存在群体分层等系统性的问题。

QQ plot 全称是 quantile-quantile plot (分位数-分位数图),是在统计学中,通过比较两个概率分布的分位数对这两个概率分布进行比较的概率图方法。分位数是指将一个随机变量的概率分布范围分为几个等份的数值点,常用的有中位数(即二分位数)、四分位数、百分位数等。

GWAS的QQplot中,作图的对象为关联分析时得到的p值,与曼哈顿图一样,通常我们将p之转化为-log10(p)。qqplot的x轴为期望分布(expected distribution),在这里实际为(0,1)的均匀分布的分位数,y轴则是观察分布(observed distribution),也就是我们实际得到的p值的分位数。

作图的目的我们可以简单理解为,实际得到的p值是否是服从均匀分布。因为如果检验的所有snp都与表型不相关,那么p值就会成(0,1)的随机分布。

下面用python现实如何作图:

使用python绘制GWAS的曼哈顿图 Mahattan plot

曼哈顿图对GWAS结果进行可视化最常用的图之一,其横轴可理解为snp在相应染色体上的位置,纵轴则为-log10(P)。因其酷似纽约曼哈顿的天际线而被命名为曼哈顿图。

本文将介绍利用python作图的方法。

要使用的package:

pandas , numpy, matplotlib, seaborn


具体代码如下:

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np

#读取GWAS结果,至少要包含 染色体号(CHR),P值(p.value),以及确保snp被正确排序
gwas_sumstats = pd.read_csv("test.txt")  

#构建一个辅助列,作为各个snp在曼哈顿图x轴上的坐标
gwas_sumstats["i"]=gwas_sumstats.index

# (optional) 如果p值没有转换成 -log10(p) 那么需要
gwas_sumstats["LOG10_P"]=-np.log10(gwas_sumstats["p.value"])

#开始做图
plot = sns.relplot(data=gwas_sumstats, x='i', y='LOG10_P', aspect=2.3, 
                   hue='CHR', palette = 'dark', s=4, legend=None) 

#处理x轴上每个chr标注的位置,这里使用每条染色体上SNP索引的中间值
chrom_df=gwas_sumstats.groupby('CHR')['i'].median()

#在上一步计算得到的位置添加x轴的标注
plot.ax.set_xticks(chrom_df)
plot.ax.set_xticklabels(chrom_df.index)

#最后添加各种标注,辅助线,调整图像范围
plot.ax.set_xlabel('CHR')
plot.ax.set_ylim(0,20)
plot.ax.axhline(y=-np.log10(5e-8), linewidth = 2,linestyle="--",color="grey")
plot.fig.suptitle('Manhattan plot')

这样我们就得到了如下的曼哈顿图: