使用GWAS数据估计跨族裔遗传相关 popcorn – Cross-ancestry genetic correlation

popcorn 背景介绍

回顾(在单一族裔中估计遗传相关的方法):通过Bivariate LD Score regression估计遗传相关性 genetic correlation

随着GWAS等遗传变异关联研究在不同人群中不断开展,相关研究数量不断增加,这为我们提供了宝贵的机会以研究遗传结构在不同人种中的差异,有助于解决这个在医疗研究与群体遗传学中十分重要的问题。

该文章基于Brielin等的研究,简单介绍一种估计跨族裔遗传相关的方法,popcorn,用于估计在不同群体中常见SNP的因果效应量的相关性。该方法纳入了研究中所有SNP的关联信息(区别于其他需要事先进行LD-pruning的方法),而且仅使用GWAS的概括数据,减少了计算负担,并避免了无法获取基因型数据的问题。

通常情况下,在单一群体中,两个表性之间的遗传相关的定义为对于两个表型的SNP效应量大小之间的相关系数。而在多群体的情况时,因为不同群体内等位基因频率有显著差异,定义并不统一。因为某个变异可能在某个特定群体中有较高的效应量但较低的频率,所以本方法不仅考虑了效应量大小的相关,同时也考虑遗传冲击(allelic impact)的相关性

该文中,跨族裔遗传效应相关的定义为SNP效应量的相关系数,跨族裔遗传冲击相关的定义为特定群体的等位方差正态化(allele-variance-normalized)的SNP效应量之间的相关系数

直觉上,遗传效应相关性(genetic-effect correlation)衡量的是同一个SNP造成表型改变的程度差异,而遗传冲击相关性(genetic-impact correlation)则是分别在各个群体中,相比于的稀有allele,给予common allele更高的权重。举个例子,如果只看效应相关性,计算得出的相关系数是不完全的,没有反应出两个群体间的效应的完整信息,因为即使效应量的遗传相关是1,但其在两个群体中,对其各自表型分布的冲击也可能由于allele frequency的不同而不同。

所以跨族裔遗传相关的两个核心问题就被引出:

1 同一个变异在两个群体中效应量大小的差异?

2 同一个变异的效应量在两个群体中的差异是否被不同群体的allele frequency所稀释或加重?

为了估计遗传相关,popcorn采用了贝叶斯方法,假设基因型是分别从各个群体内抽取的,并且效应量服从正态的先验分布(无穷小模型)。无穷小模型的假设会产生一个观测统计量(Z scores)的多元正态分布,其中该分布的协方差是遗传力与遗传相关的一个方程。这使得我们可以直接将Z score的膨胀纳入模型中,而不需要进行LD pruning。然后popcorn会最大化近似加权似然方程,来估计遗传力与遗传相关性。

方程推导可参考原文章,这里我们要估计的是以下两个值:

遗传效应相关:pge=Cor(β1,β2)

遗传冲击相关:pgi=Cor(δ1β1,δ2β2)

其中δ为allele的标准差

使用方法:

该文章作者的github仓库: https://github.com/brielin/Popcorn

要求:

popcorn只能使用python2,建议用conda新建一个环境
需要以下python包:
numpy 1.14.2
scipy 1.0.1
pandas 0.22.0
pysnptools 0.3.9
bottleneck 1.0.0
statsmodels 0.8.0
(to use --plot_likelihood) matplotlib 1.5.1

从github clone之后,安装:

cd Popcorn
python setup.py install

可以下载使用千人基因组项目提前计算好的LD分数 (用于EUR-EAS):

https://www.dropbox.com/sh/37n7drt7q4sjrzn/AAAa1HFeeRAE5M3YWG9Ac2Bta .

popcorn 方法分为两步,

第一步为计算LD分数(使用下载好的数据可以跳过这一步,作者提供了EUR-EAS的LD分数)

第二步为计算遗传相关性

第一步,计算LD分数(可选)

python -m Popcorn compute -v 2 --bfile1 Popcorn/test/EUR_ref --bfile2 Popcorn/test/EAS_ref --gen_effect Popcorn/test/EUR_EAS_test_ge.cscore

-v 表示 输出log的详细程度

—bfile表示 输入的plink bfile文件前缀

—gen_effec 表示计算遗传效应相关 (计算遗传冲击相关时不加这个flag)

最后是输出的路径和文件名

第一步log文件

Popcorn version 0.9.6
(C) 2015-2016 Brielin C Brown
University of California, Berkeley
GNU General Public License v3

Invoking command: python /home/brielin/CoHeritability/Popcorn/__main__.py compute -v 2 --bfile1 Popcorn/test/EUR_ref --bfile2 Popcorn/test/EAS_ref Popcorn/test/EUR_EAS_test.cscore
Beginning analysis at DATE
50000 SNPs in file 1
50000 SNPs in file 2
49855 SNPs in file 1 after MAF filter
44021 SNPs in file 2 after MAF filter
43898 SNPs common in both populations
T1
T2
TC
Analysis finished at DATE

结果的格式如下所示:

1	14930	rs75454623	A	G	0.520874751491	0.58630952381	1.4252539451	1.78489131664	1.49214493091
1	15211	rs78601809	T	G	0.731610337972	0.503968253968	1.20243142615	2.3876817626	1.29403313474
1	15820	rs200482301	T	G	0.728628230616	0.605158730159	1.56556914093	1.42934198599	1.40055911225
1	55545	rs28396308	T	C	0.813121272366	0.804563492063	2.59950434496	2.1809181648	2.06610287872

各列分别代表Chromosome, Base Position, SNP ID, Allele 1, Allele 2, Frequency of A1 in POP1, Frequency of A1 in POP2, LD score in POP1, LD score in POP2, and Cross-covariance score

第二步,计算相关

python -m Popcorn fit -v 1 --use_mle --cfile Popcorn/test/EUR_EAS_test_ge.cscore --gen_effect --sfile1 Popcorn/test/EUR_test.txt --sfile2 Popcorn/test/EAS_test.txt Popcorn/test/EAS_EUR_test_ge_result.txt

–sfile1与–sfile2为输入的两个群体的GWAS文件,注意要与LD分数中POP1和POP2的顺序相一致

输入的GWAS数据至少包含一下列(列名应与下面的介绍保持一致):

1 rsid 或者 SNP (与LD分数文件中的id保持一致)

注: 可以使用 chr 以及 pos 来替代

2 a1 和 a2

3 N 每个SNP的样本大小

4 beta和SE 或者 OR和p-value 或者 Z分数 (不要求对应a1为effect allele,但要保持所有snp一致即可)

5 AF (可选)即a2的allele frequency,用于A/T , G/C的SNP的flip

输入文件例如:

rsid	a1	a2	af	N	beta	SE
rs10458597	C	A	0.59421	50000	-0.0152782557183	0.00912151085515
rs3131972	A	C	0.83102	50000	-0.000638379272685	0.0119663325893
rs3131969	A	C	0.8654	50000	-0.000737906302463	0.0131325861467
rs3131967	A	C	0.87044	50000	-0.000189441590945	0.0133506462369

执行代码后的,第二步log文件:

Popcorn version 0.9.6
(C) 2015-2016 Brielin C Brown
University of California, Berkeley
GNU General Public License v3

Invoking command: python /home/brielin/CoHeritability/Popcorn/__main__.py fit -v 1 --cfile Popcorn/test/EUR_EAS_test_ge.cscore --gen_effect --sfile1 Popcorn/test/EUR_test.txt --sfile2 Popcorn/test/EAS_test.txt Popcorn/test/EAS_EUR_test_ge_result.txt
Beginning analysis at DATE
Analyzing a pair of traits in two populations
49999 49999 SNPs detected in input.
43897 SNPs remaining after merging with score file.
43897 SNPs remaining after filtering on AF and self-compliment SNPs.
Initial estimate:
       Val (obs)  SE
h1^2   0.172133 NaN
h2^2   0.086708 NaN
pg     0.302372 NaN
Performing jackknife estimation of the standard error using 200 blocks of size 219 .
Jackknife iter: 200
      Val (obs)        SE          Z     P (Z)
h1^2   0.172133  0.012870  13.374388      0
h2^2   0.086708  0.008555  10.135759      0
pge    0.302372  0.053021  13.157553      0
Analysis finished at DATE
Total time elapsed: T

pge 则为跨族裔遗传相关的估计值 (加–gen_effect时的结果)

pgi 为跨族裔遗传冲击的估计值(不加–gen_effect时的结果)

这里的p值的原假设是估计值等于1,p<0.05则意味着估计值显著小于1

参考:

https://pubmed.ncbi.nlm.nih.gov/27321947/

使用GWAS数据进行组织细胞特异性分析 – LDSC-SEG

目录

LDSC-SEG简介
LDSC-SEG原理
LDSC-SEG教程
参考

LDSC-SEG简介

本文主要介绍LD分数回归的一项功能扩展 ,LDSC-SEG。该方法主要通过整合基因表达数据与GWAS数据,来发现疾病相关的组织或细胞种类,其底层原理是使用分层LD分数回归(stratified LD score regression,回顾:GWASLab:基于功能分类分割遗传力 – 分层LD分数回归 Stratified LD score regression)来检测表型的遗传性是否富集于某个组织或细胞高特异表达基因的周围。

LDSC-SEG原理

LDSC-SEG的主要步骤 (如图所示)

  1. 从一个标准化后的基因表达矩阵(包含不同组织与细胞)开始
  2. 对于每个gene,首先计算其在目标组织中特异表达的t统计量
  3. 根据t统计量从大到小对所有基因进行排序,并将位于t统计量最高的前10%的基因定义为该组织特异表达的基因集。
  4. 在这些特异表达的基因的转录区域两侧加上100kb的窗口,并以此为构建组织特异的基因组注释。
  5. 最后对于以上注释后的不同组别进行分层LD回归来估计不同组织或细胞对于所研究表型遗传力的贡献

联合回归模型中同时包括

1)上述目标组织或细胞特异的基因集,

2)一个包含基因组上所有基因的基因集,

3)包含52种注释的基线模型(这些注释包括增强子区域,保守性区域等,可参考 stratified LD score regression)。

特异基因集的回归系数为正数,则表示矫正其他基因集后,该特异基因集对表型遗传性的贡献为正。

基因表达数据来源(原文中处理了五组公开的数据集,如下所示)

GTExHuman53 tissues/cell typesRNA-seq
Franke labHuman/mouse/rat152 tissues/cell typesArray
CahoyMouse3 brain cell typesArray
PsychENCODEHuman2 neuronal cell typesRNA-seq
ImmGenMouse292 immune cell typesArray

实战操作

ldsc 下载安装可参考: GWASLab:连锁不平衡分数回归 LD score regression -LDSC

教程原文:

https://github.com/bulik/ldsc/wiki/Cell-type-specific-analyses

#以上述5组数据中的Cahoy数据为例:
cts_name=Cahoy 

#下载对应的LD分数文件 (EUR)

#组织、细胞特异LD分数
wget https://data.broadinstitute.org/alkesgroup/LDSCORE/LDSC_SEG_ldscores/${cts_name}_1000Gv3_ldscores.tgz
#基线模型文件
wget https://data.broadinstitute.org/alkesgroup/LDSCORE/1000G_Phase3_baseline_ldscores.tgz
#权重文件
wget https://data.broadinstitute.org/alkesgroup/LDSCORE/weights_hm3_no_hla.tgz

#解压
tar -xvzf ${cts_name}_1000Gv3_ldscores.tgz
tar -xvzf 1000G_Phase3_baseline_ldscores.tgz
tar -xvzf weights_hm3_no_hla.tgz

#下载GWAS数据
wget https://data.broadinstitute.org/alkesgroup/UKBB/body_BMIz.sumstats.gz
#下载SNP list
wget https://data.broadinstitute.org/alkesgroup/LDSCORE/w_hm3.snplist.bz2
#解压
bunzip2 w_hm3.snplist.bz2

#第一步清洗数据 
python munge_sumstats.py \
--sumstats body_BMIz.sumstats.gz \
--merge-alleles w_hm3.snplist \
--out UKBB_BMI

#第二步 进行LDSC-SEG
cts_name=Cahoy
ldsc.py \
    --h2-cts UKBB_BMI.sumstats.gz \
    --ref-ld-chr 1000G_EUR_Phase3_baseline/baseline. \
    --out BMI_${cts_name} \
    --ref-ld-chr-cts $cts_name.ldcts \
    --w-ld-chr weights_hm3_no_hla/weights.

注意:

**--ref-ld-chr-cts 这里输入一个ldcts文件,共两列,第一列为一个便于识别的label,第二列则是逗号分隔的两个路径,第一个为被检测组织细胞的ld,第二个则为control的ld,路径格式与--ref-ld-chr的格式类似。注意这里的都是相对与解压文件夹的相对路径,使用前最好修改为绝对路径。

如下图所示:

#组织细胞标签 检验组LD分数路径,对照组LD分数路径
Astrocyte       Cahoy_1000Gv3_ldscores/Cahoy.1.,Cahoy_1000Gv3_ldscores/Cahoy.control.
Oligodendrocyte Cahoy_1000Gv3_ldscores/Cahoy.2.,Cahoy_1000Gv3_ldscores/Cahoy.control.
Neuron  Cahoy_1000Gv3_ldscores/Cahoy.3.,Cahoy_1000Gv3_ldscores/Cahoy.control.

执行完成后可以查看结果:

.cell_type_results.txt

Name    Coefficient     Coefficient_std_error   Coefficient_P_value
Neuron  7.93194788527e-09       3.02894244784e-09       0.00441303625204
Oligodendrocyte 7.32019970874e-10       3.51868270994e-09       0.417599619801
Astrocyte       -5.76220451287e-09      2.60400594455e-09       0.98654507806

第一列即为dict文件里的label,第二,三列为ldcts文件里第二例第一个路径对应的ld分数的系数与se(对应cell type specific annotation),第四列为该系数的p值。

参考

https://github.com/bulik/ldsc/wiki/Cell-type-specific-analyses

https://www.ncbi.nlm.nih.gov/pmc/ar

一行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

去哪找公开的GWAS数据 – Publicly Available GWAS Sumstats

目录

  1. 公共数据库
  2. 各大Biobank与Cohort
    • 欧美
    • 东亚
    • Global Biobank
  3. 单独的研究发表的数据
  4. 各类疾病研究组织的数据

本文将简要介绍可以简单入手的GWAS sumstats的一些来源,文中提及的来源都是公开可下载的,适合新入门的同学练手,做一些小规模的meta分析,或是post-gwas分析。

1.公共数据库:

目前收录最全的GWAS数据库,但只收录了已发表GWAS的数据,而且会有几个月到半年左右的延迟,使用时应该注意:

GWAS catalog: https://www.ebi.ac.uk/gwas/

OpenGWAS: https://gwas.mrcieu.ac.uk/

2.各大Biobank的公开数据:

规模较大的Biobank或Cohort会定期公开GWAS sumstats,这类数据有些未被收录在GWAS catalog中,一些gwas检验方法的作者也会使用其方法进行Biobank级别的GWAS分析,这类数据的表型覆盖较广:

2.1欧美:

Finngen 4 : https://r4.finngen.fi/about

Finngen 5 : https://r5.finngen.fi/about

Finngen 6 : https://r6.finngen.fi/about

UKB : https://pheweb.org/UKB-Neale/

UKB saige: https://pheweb.org/UKB-SAIGE/

UKB fastgwa-glmm: https://yanglab.westlake.edu.cn/resources/ukb_fastgwa/imp_binary/

UKB fastgwa: https://yanglab.westlake.edu.cn/resources/ukb_fastgwa/imp/

UKB TOPMed: https://pheweb.org/UKB-TOPMed/

UKB gene-based: https://genebass.org/

Pan-UKB : https://pan.ukbb.broadinstitute.org/

MGI 1 : https://pheweb.org/MGI-freeze1/

MGI 2 : https://pheweb.org/MGI-freeze2/

MGI BioUV: https://pheweb.org/MGI-BioVU/

FinMetSeq: https://pheweb.sph.umich.edu/FinMetSeq/

Generation Scotland: https://datashare.ed.ac.uk/handle/10283/844

2.2 东亚:

Biobank Japan:

JENGER: http://jenger.riken.jp/result

Pheweb: https://pheweb.jp/

ToMMo – JMorp:https://jmorp.megabank.tohoku.ac.jp/202109/gwas/

KoreanChip: https://www.koreanchip.org/downloads

KoGES Pheweb: https://koges.leelabsg.org/

2.3 全球范围Biobank的Meta分析:

近期Global Biobank项目也公开了一批全球范围biobank的常见复杂疾病meta分析

Global Biobank :http://results.globalbiobankmeta.org/

3.单独的研究发表的数据

直接在google上搜索研究论文,从Data availability里的url顺藤摸瓜找到数据,此类数据散落在各个学校,研究机构等等自家的网站上,没有好的办法,只有自己搜索。GWAS catalog 经常会出现收录不及时而漏掉很多最新数据。

例如:

Program in Complex Trait Genomics, IMB, The University of Queensland.

https://cnsgenomics.com/content/data

https://ctg.cncr.nl/software/summary_statistics

4.各类疾病研究组织的数据

这类数据通常为meta分析后的sumstats,例如:

DIAGRAM:http://www.diagram-consortium.org/downloads.html

Megastroke: https://www.megastroke.org/index.html

GIANT (Genetic Investigation of ANthropometric Traits):https://portals.broadinstitute.org/collaboration/giant/index.php/Main_Page

GLGC (Global Lipids Genetics Consortium):  http://csg.sph.umich.edu/willer/public/glgc-lipids2021/

PGC (Psychiatric Genomics Consortium): https://www.med.unc.edu/pgc/download-results/

等等

一个人的总结难免会有疏漏,欢迎大家在评论里补充!

使用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