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

GWAS Sumstats Harmonization GWAS数据的协调统一

背景介绍

在进行meta分析之前,我们首先要对gwas的sumstats进行预处理,这一步看似简单,但却是高质量meta分析中必不可少的一步,本文将简单介绍预处理过程中需要注意的事项,并简单介绍一款实用的预处理软件。

不同GWAS sumstats通常格式千差万别,单单统一格式就已经让人烦恼,而更令人万念俱灰的是,不同sumstats的参考的基因组版本并不相同,allele的正负链也不统一,SNPID有的用rsID(rsID的版本又是个问题~_~),有的又用chrpos…

针对以上种种的不统一,我们都要逐一进行处理,要注意事项主要包括:

  1. 参考基因组的版本:hg19还是hg38 或其他
  2. 表示变异等位的正负链是否统一
  3. 如果是a1,a2这种表示形式,则需要明确effect allele与non-efffect allele
  4. 对于palindromic SNP (A/T, 或是G/C如何处理?)
  5. 表示效应量的是beta还是OR

Pipeline

预处理简要的pipeline:

  1. 确定参考所用的VCF文件, 最好包含相应群体的allele frequency(AF)数据
  2. 对于基因组版本不一致的sumstats首先liftover转换成参考vcf的基因组版本
  3. 对sumstats进行整理并统一beta与OR的使用,
    1. 数量表型用beta与SE
    2. 二分类表型用OR与OR_95L,OR_95U
  4. 对不同群体的sumstats,使用相应的AF,利用harmoniser工具分别进行harmonise
    1. 设定依据AF判断正负链的阈值
  5. 对于处理后的gwas sumstats进行统一的质控 QC
    1. INFO
    2. MAF等

Genetics-Sumstats-Harmoniser工具使用教程

下载与安装

harmonise可以使用的工具

该工具由opentargets项目组成员开发

github链接:https://github.com/opentargets/genetics-sumstat-harmoniser

# Download 从github下载
git clone <https://github.com/opentargets/sumstat_harmoniser.git>
cd sumstat_harmoniser

# Install dependencies into isolated environment 创建环境
conda env create -n sumstat_harmoniser --file environment.yaml

# Activate environment 激活环境
source activate sumstat_harmoniser

# Run tests 检验是否安装成功
cd tests
./run_tests.py

环境要求 Requirements

参考数据 Reference data

可以使用的参考数据:

以global biobank meta-analysis initiative为例,其参考数据使用了gnomAD的数据:

https://gnomad.broadinstitute.org/downloads

示例:

INFO字段中包含了不同族裔中的allele frequency的信息:

#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
1	10108	rs62651026	CAACCCT	C	46514.3	RF	AF_eas=0;AF_nfe=0;AF_fin=0.0227273
1	10109	rs376007522	AACCCT	A	89837.3	RF	AF_eas=0;AF_nfe=0.0512821;AF_fin=0.0588235

使用方法

输入选项

需要两个文件:一是gwas的sumstats,二是参考的vcf文件

Input files:
  --sumstats <file>     GWAS summary statistics file
  --vcf <file>          Reference VCF file. Use # as chromosome wildcard.

其中gwas sumstats还需要指定对应的列名:

Input column names:
  --chrom_col <str>     Chromosome column 染色体
  --pos_col <str>       Position column 碱基位置
  --effAl_col <str>     Effect allele column 效应等位
  --otherAl_col <str>   Other allele column 非效应等位
  --beta_col <str>      beta column
  --or_col <str>        Odds ratio column 
  --or_col_lower <str>  Odds ratio lower CI column
  --or_col_upper <str>  Odds ratio upper CI column
  --eaf_col <str>       Effect allele frequency column
  --rsid_col <str>      rsID column in the summary stat file

输出选项

输出3个文件,Harmonise后的文件,Harmonise过程中的统计文件,以及forward/reverse/palindromic variants的数量的文件,通过以下选项指定:

Output files:
  --hm_sumstats <file>  Harmonised sumstat output file (use 'gz' extension to
                        gzip)
  --hm_statfile <file>  Statistics from harmonisation process output file.
                        Should only be used in conjunction with --hm_sumstats.
  --strand_counts <file>
                        Output file showing number of variants that are
                        forward/reverse/palindromic

其他选项

Other args:
  --only_chrom <str>    Only process this chromosome
  --in_sep <str>        Input file column separator (default: tab)
  --out_sep <str>       Output file column separator (default: tab)
  --na_rep_in <str>     How NA are represented in the input file (default: "")
  --na_rep_out <str>    How to represent NA values in output (default: "")
  --chrom_map <str> [<str> ...]
                        Map summary stat chromosome names, e.g. `--chrom_map
                        23=X 24=Y`

Harmonisation mode 模式选择:

可以选择以下四种模式:[infer|forward|reverse|drop]

对于 palindromic variants 的处理方式为:

(a) infer strand from effect-allele freq, 根据AF推测

(b) assume forward strand, 假设全部都为正链

(c) assume reverse strand, 假设全部为负链

(d) drop palindromic variants 舍弃所有palindromic variants

AF字段名与阈值

如果选择通过AF推测正负链(–palin_mode infer),那么还需要制定vcf文件中af的字段名与推测的阈值 –af_vcf_field <str> VCF info field containing alt allele freq (default:AF_NFE)

指定VCF文件中INFO里包含AF信息的字段

–infer_maf_threshold <float>Max MAF that will be used to infer palindrome strand (default: 0.42)

推测palindrome strand的最大MAF的阈值

使用例

如果输入的sumstats为以下格式:

CHROMOSOME      POSITION        MARKERNAME      NEA     EA      EAF     INFO    N       BETA    SE      P       OR      OR_95L  OR_95U

那么harmonise时,可以使用如下的代码:

sumstat_harmoniser \
 --sumstats ${sumstats} \
 --vcf gnomad.vcf.gz \    #以gnomad为参考
 --eaf_col EAF \
 --af_vcf_field AF_eas \
 --chrom_col CHROMOSOME \
 --pos_col POSITION \
 --effAl_col EA \
 --otherAl_col NEA \
 --beta_col BETA \
 --or_col OR \
 --or_col_lower OR_95L \
 --or_col_upper OR_95U \
 --palin_mode infer \
 --hm_sumstats ${out}.hm \
 --hm_statfile ${out}.stats

harmonise后的结果为原文件前加上以下harmonise后的列

hm_varid hm_rsid hm_chrom hm_pos hm_other_allele hm_effect_allele hm_beta hm_OR hm_OR_lowerCI hm_OR_upperCI hm_eaf hm_code

参考

https://gnomad.broadinstitute.org/downloads

MAGMA 基因及通路分析

目录

1 背景
2 算法简介
3 使用教程
 3.1 magma软件与参考数据下载
 3.2 注释
 3.3 基因关联分析
 3.4 通路分析
4 参考

背景

通过以有生物学意义的方式整合复杂疾病的信息的基因以及通路分析,是单变异GWAS的有效补充。

本文将介绍MAGMA,一个基于多变量回归模型的基因以及通路分析工具。

算法简介

Gene-based分析

magma的gene-based分析模型采用了多元线性主成分回归(multiple linear principal components regression)

α0g为截距,Xg* 为PC矩阵,W与βg为协变量变量的效应(可选),αg为遗传效应,εg为残差

计算时:

首先根据染色体位置提取某个基因SNP矩阵,计算PC

去除掉特征值过小的PC

对剩下的主成分进行回归

最后通过F检验获得p值 (H0: αg =0 )

Gene set /Pathway 基因集 通路分析

进行通路分析时,magma首先将上一步所得到的每个基因的p值通过probit方程转化为z值,

这里的zg整体上近似服从正态分布,zg反映了该基因关联的强度。

Self-contained gene-set analysis检验一个通路上的基因的整体上是否与表型相关联。使用通路里基因的z值,可以构建一个只有截距项的回归,然后检验β是否大于0:

Competitive gene-set analysis检验在一个通路里的基因,是否比这个通路里其他基因的相关性更高。

因为通路分析中的回归模型为标准线性回归,假设残差项独立正态分布,即

但实际上,由于LD,相邻的基因可能会互相关联,而打破上述假设。所以magma采用了广义最小二乘法,并假设

R为基因-基因的关联矩阵。

概括性数据分析

在没有原始基因型数据的情况下,magma也提供分析概括性数据SNP的p值的方法,即逐个分析某个基因中的SNP,然后将SNP的p值合并成该基因的检验统计量,但是该方法需要提供一个相近群体的LD参考面板。

使用教程

magma软件与参考数据下载

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

1.软件下载,下载对应系统的版本即可

2.基因注释参考,magma提供了多种版本,下载对因自己数据的版本

3.LD参考面板,来自1000 genome,下载相应人群的文件即可

magma三步骤:

  1. 注释:将自己的SNP根据染色体位置注释到基因上
  2. Gene-based 分析
  3. 基因集/通路分析

SNP注释 Annotation

第一步需要对gwas中所包含的SNP进行注释,在这里就是将SNP根据染色体上的位置对应到相应基因上。

示例代码如下:

magma \
	--annotate \
	--snp-loc [SNPLOC_FILE] \
	--gene-loc [GENELOC_FILE] \
	--out [ANNOT_PREFIX]

snp-loc 文件应当包含三列,SNPID,CHR以及POS,这里也可以直接使用plink的bim文件

rs540836310	1 861725
rs139858754	1 861728

gene-loc 可以直接使用magma提供的基因注释参考文件 (注意版本)

注释完成后会得到 .genes.annot文件,内容为第一例为geneID 之后为在这个基因内的SNP

816	rs540836310	rs139858754	rs188152259	rs13302982	rs76842830	rs13303101	rs74442310	rs6680268

基因关联分析 Gene-based analysis

第二部,进行基于基因的关联分析,

示例代码如下所示:

magma \
	--bfile [REFDATA] \
	--pval [PVAL_FILE] \
	N=[N] \
	--gene-annot [ANNOT_PREFIX].genes.annot \
	--out [GENE_PREFIX]

bfile为原始数据或是参考LD面板,如果数据量不大可以直接使用自己的plink的bed格式原始数据,在原始数据无法获得的时候可以使用magma提供的1000 genome参考数据,biobank级别的数据的情况下,可以随机抽取某个族裔无亲缘关系的一定人数(例如20000人)来构建自己的参考面板。

pval为SNP的p值文件,包含两列 SNP 以及 P, 这里N为样本量,

gene-annot为上一步注释后得到的文件。

完成后有两个文件输出;1 .genes.out 2 ..genes.raw

1 .genes.out 基因关联分析的结果

GENE       CHR      START       STOP  NSNPS  NPARAM       N        ZSTAT            P
148398       1     859993     879961    125      27  157848      -1.7291       0.9581

2 ..genes.raw 在第三步通路分析时会使用到这个文件

# VERSION = 109
# COVAR = NSAMP MAC
148398 1 859993 879961 125 27 157848 133.304 -1.72907

基因集 通路分析 gene-set analysis 、 pathway

第三步,通路分析

这一步是对显著的基因是否富集于某个基因集或通路进行检验,

如果是对经典的通路进行检验,可以在MSigDB(GSEA)下载通路文件 msigdb.v7.4.entrez.gmt

MSigDB:https://www.gsea-msigdb.org/gsea/downloads.jsp#msigdb

示例代码如下:

magma \
	--gene-results [GENE_PREFIX].genes.raw \
	--set-annot [SET_FILE] \
	--out [GS_PREFIX]

gene-results 为上一步基因关联分析所得的.genes.raw文件

set-annot 为基因集或通路的定义文件,可以直接使用MSigDB下载的文件(注意基因ID要与之前相对应)

默认输出三个.gsa.out .gsa.genes.out, .gsa.sets.genes.out

.gsa.out 为最主要的输出文件,包含了各个通路的检验结果

# TOTAL_GENES = 17198
# TEST_DIRECTION = one-sided, positive (set), two-sided (covar)
# CONDITIONED_INTERNAL = gene size, gene density, inverse mac, log(gene size), log(gene density), log(inverse mac)
VARIABLE                                TYPE  NGENES         BETA     BETA_STD           SE            P FULL_NAME
chr1p12                                  SET      18    -0.096902     -0.00381      0.30843      0.71093 chr1p12


.gsa.genes.out 为在此分析中使用到的基因关联结果

GENE       CHR      START       STOP  NSNPS  NPARAM       N        ZSTAT  ZFITTED_BASE  ZRESID_BASE
148398       1     859993     879961    125      42  12345      -1.62             0      -1.62

.gsa.sets.genes.out 为多重检验调整后, 显著的通路里的基因的信息

参考

MAGMA: Generalized Gene-Set Analysis of GWAS Data. de Leeuw CA, Mooij JM, Heskes T, Posthuma D (2015) MAGMA: Generalized Gene-Set Analysis of GWAS Data. PLOS Computational Biology 11(4): e1004219. https://doi.org/10.1371/journal.pcbi.1004219

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

多基因风险分数 PRS( Polygenic risk score)系列之七:Pathway-based PRS 通路PRS

本文是多基因风险分数PRS系列的第七篇文章,主要介绍基于通路的PRS(Pathway-based PRS)的概念,以及简要的计算方法(仅供参考)。

回顾:

基于通路的PRS (Pathway-based PRS )简介

目前主流PRS计算模型的局限

当前主流的PRS计算方法都是基于疾病的经典多基因模型,也就是假设个体处在一个从低到高风险的线性分布上,其PRS作为易感性(liability)的估计值,能够总结这个个体的遗传概况。

易感性:易感性-阈值模型 Genetic liability, Threshold model

但这种方法也存在显著的缺陷,那就是这种概括性的估计值不能完全捕捉个体遗传特征,存在明显的信息丢失,一个最显著的例子就是无法判别不同生物学过程与通路之间遗传风险大小的变化,但这类信息通常在PRS的应用中(例如病人分层,或是治疗反应预测)起到重要的作用。

为了解决上述问题,近年来基于通路的PRS方法逐渐获得关注,本文将主要介绍Paul O’Reilly课题组开发的Pathway-based PRS方法PRSet。PRSet方法可以解释基因组中的亚结构(genomic sub-structure),该模型在经典多基因模型中进行扩展,使其能更好反映疾病的异质性。

基于通路的PRS : Pathway-based PRS

经典PRS模型计算整个基因组上risk allele与其效应估计乘积的加和,而pathway-based PRS只分别纳入k个相关通路上的SNP。也就是说,使用pathway-based PRS方法时,一个个体会有k个不同通路的通路PRS,每个通路PRS都对应基因组上一个具体的通路。(如图a)

在这里,通路可以由多种方式定义,可以是一些经典的通路数据库(如:KEGG,REACTOME等),也可以是通过基因共表达,蛋白质互作,实验扰动的功能产出等分析得到的模组。

理想状态下,通路能够反映不同生物学功能的编码。这些功能是可分割的,就像流行病学研究中,不同环境风险因素(抽烟,饮食,等等)也是分割的一样。基于这种观点,我们可以吧GWAS的结果也视为不同通路编码的生物学功能的信号的综合。(如图b)

PRSet的通路PRS方法简要原理

PRSet采用传统的C+T方法,分别对k条通路分别计算PRS,mk为每条通路上SNP数量,对于每个个体i,计算后会得到对应k个通路的k个PRS:

Clumping

与传统C+T对基因组上所有SNP进行clumping的方法不同,PRSet独立地对每条通路上的SNP进行clumping,确保能最大程度的保留信号。PRSet采用如下所示的bit-flag方法,对不同通路独立clumping。简单来说就是只有indexSNP包含在该gene set的时候才会去除掉需要被clump的SNP,否则会保留。

PRSet软件链接:https://www.prsice.info/prset_detail/

注:该方法的优势是简单快速,而且已经在软件中实现,作为探索可以尝试,但具体效果还有待检验。

(20211211:目前该文章还没有正式发表,文章里的检验方法不能令人信服,所以仅供参考)

参考:https://www.prsice.info/prset_detail/

https://europepmc.org/article/ppr/ppr362752

Linux :输入/输出重定向 >, 1>, 2>, &>, >> , <<

linux重定向可以将命令的输出或输入重新定向到其他位置或文件,以实现对输出输入的控制。默认情况下命令的输出通常为终端,如果想将输出转移到文件或其他位置,这时候就需要重定向。

文件描述符:一个命令通常都会打开三个文件,默认使用文件描述符0,1,2来指代这个三个文件

stdin   0	标准输入流 (键盘)
stdout	1	标准输出流 (终端)
stderr	2	标准错误输出流 (终端)

输出重定向

输出重定向常见形式:

command > file  #将标准输出1重定向到 file 里
command 1> file #将标准输出1重定向到 file 里,与上面的写法功能一样
command 2> file #将标准错误输出1重定向到 file 里
command &> file #将标准输出1 与 标准错误输出2 一起重定向到 file 里

最简单的用法,将命令的输出重定向到文件里:

$ ls > out.txt #将ls的结果写进out.txt里
$ cat out.txt 
file1.txt
out.txt

当然输出重定向还有其他用法,我们来看下面看这个例子:

当前目录下只有file1.txt这个文件

但我们使用cat命令查看file1.txt和file2.txt,返回终端的结果其实包含两部分

cat file1.txt的标准输出,

cat file2.txt的标准错误输出(因为文件不存在)

$ ls
file1.txt

$ cat file1.txt file2.txt 
1	apple      #标准输出
2	banana     #标准输出
3	pear       #标准输出
4	orange     #标准输出
cat: file2.txt: No such file or directory   #标准错误输出

1> 和 2> 分别重定向标准输出,与标准错误输出

# 单独使用> 或1>将命令的结果重定向到 out.txt文件中
# 此时,只有标准输出被重定向,标准错误输出被直接输出到终端
$ cat file1.txt file2.txt > out.txt
cat: file2.txt: No such file or directory #标准错误输出
$ cat out.txt 
1	apple
2	banana
3	pear
4	orange

# 单独使用2> 时,只有标准错误输出被重定向,标准输出被直接输出到终端
$ cat file1.txt file2.txt 2>out.txt
1	apple
2	banana
3	pear
4	orange
$ cat out.txt 
cat: file2.txt: No such file or directory

&> 为一起重定向标准输出 与 标准错误输出的简便写法

# &> 是指将 标准输出 与 标准错误输出 一起重定向 
$ cat file1.txt file2.txt &> out.txt
$ cat out.txt 
1	apple
2	banana
3	pear
4	orange
cat: file2.txt: No such file or directory

#以上代码等价于 1> out.txt 2>&1
#这里 &1 指代标准输出1
#将标准错误输出2 先定向到标准输出1,然后将标准输出1定向到out.txt文件中
$ cat file1.txt file2.txt 1> out.txt 2>&1
$ cat out.txt 
1	apple
2	banana
3	pear
4	orange
cat: file2.txt: No such file or directory

#等价于
$ cat file1.txt file2.txt 2> out.txt 1>&2
$ cat out.txt 
1	apple
2	banana
3	pear
4	orange
cat: file2.txt: No such file or directory


不覆盖file的输出重定向 >>

command >> file 
#将标准输出1重定向到 file 里, 不替代原本file文件内容,只是在末尾添加

$ cat file1.txt >out.txt #将cat file1的结果输出到out.txt
$ cat out.txt 
1	apple
2	banana
3	pear
4	orange

$ cat file1.txt >>out.txt #将cat file1的结果添加到out.txt末尾,不覆盖源文件
$ cat out.txt  
1	apple
2	banana
3	pear
4	orange
1	apple
2	banana
3	pear
4	orange

>>的其他用法与 >类似, 也可以有 &>> 1>> 2>>等用法:

#将标准输出定向到out.txt,然后将标准错误输出添加到out.txt末尾
$ cat file1.txt file2.txt 1> out.txt 2>>out.txt
(base) [heyunye@gc017 redirection]$ cat out.txt 
1	apple
2	banana
3	pear
4	orange
cat: file2.txt: No such file or directory


输入重定向

与输出重定向类似,将文件内容重定向到标准输入0

command <file  #重定向标准输入
command << xxx  #here document

看下面这个例子,使用head -v 来查看文件,理解其中区别

# 使用head -v命令来查看file1.txt (-v 选项可以输出文件名)
$ head -v file1.txt 
==> file1.txt <==           
1	apple
2	banana
3	pear
4	orange

# 使用重定向,可以看到主体部分相同,但文件名为standard input,因为我们将文件重定向到了标准输入
$ head -v <file1.txt 
==> standard input <==
1	apple
2	banana
3	pear
4	orange

Here document

可以以交互形式实现快捷的多行输入。

<<xxx 这种形式被称为Here document,xxx为任意字符串,作为标签,为Here documen的起始,输入时直接在终端里输入多行内容,完成后再次输入xxx,标记输入完成。

#here document
command <<标签  
>内容
>内容
>...
>标签

$ head -v <<abc #<<abc 是指here document的起始
> 123 #内容
> 123 #内容
> 123 #内容
> abc #abc 表示结束
==> standard input <==  #Here document被定向为标准输入
123
123
123

Linux: 进程替代 >(command) <(command)

进程替换(process substitution):

将临时文件作为某个进程的输入或输出的一种重定向。

<(command) : 将括号()里的命令执行后,输出作为一个临时文件,这个临时文件通常存在 /dev/fd/目录下

>(command) : 将之前命令的输出视为一个临时文件,作为括号()里命令的输入

先看下面这个例子:

$ echo <(ls)  
/dev/fd/63

# 先执行ls命令,结果输出到临时文件中(/dev/fd/63)
# 这里<(ls)就指代这个文件的路径(/dev/fd/63)
# 这个文件的内容是ls命令的结果
# 所以 echo <(ls) 约等于
# echo /dev/fd/63

$ ls > >(cat)
#将ls命令的结果输出到/dev/fd/目录下某个临时文件中
#将临时文件作为括号内命令cat的输入 执行
#结果上和ls单独执行一样

也就是说进程替换实际上就是将进程与/dev/fd/目录下的某个临时文件相关联,然后这个临时文件可作为后续命令的输入输出。

实际操作中,使用进程替换可以大幅简化代码,

诸如下面这个例子:

#分别提取三个文件的第一行,然后paste
$ paste -d' '  <(cut -f 1 file1.txt ) <(cut -f 1 file2.txt ) <(cut -f 1 file3.txt )
1 1 pear 1 pear x
2 2 apple 2 apple y
3 3 banana 3 banana z
4 5 kiwi 5 kiwi q

进程替换可以嵌套