Meta-analysis是对多个GWAS分析结果进行综合评价。METAL是GWAS meta分析最常用的工具之一，说明文档见： METAL_Documentation

该软件的安装非常简单，直接下载编译好的二进制文件即可，安装过程如下：

wget http://csg.sph.umich.edu/abecasis/metal/download/Linux-metal.tar.gz
tar xzvf Linux-metal.tar.gz
cd generic-metal/
在安装目录，有一个名为metal的可执行文件，该程序用法很简单，只需要编写一个配置文件，然后执行即可，所以关键在于配置文件的编写。在软件的安装目录，有一个名为example的文件夹，提供了两个示例，其中的metal.txt就是配置文件。
 
基本应用介绍：
 
1. 输入文件的分隔符：
 
METAL希望将每组结果汇总在一个表中。此表必须存储在文本文件中，但METAL在列分隔符、列标题等细节方面非常灵活。这确实意味着在进行meta分析之前需要的基本信息是每个输入文件的描述。应该指定的第一件事就是列分隔符。默认情况下，METAL假定列由空格 (由空格和制表符的任意组合组成) 分隔。也可以进行指定：
 
SEPARATOR  WHITESPACE    # 默认
SEPARATOR  COMMA         # 用于在某些平台流行的以逗号分隔的文件
SEPARATOR  TAB           # 由单个制表符分隔的列，因此连续的制表符表示空列
2. 软件支持以下两种meta分析的算法：
 
1) pvalue
 2) standard error
 
第一种是基于p值；第二种是基于标准误，我们知道标准误指的是某个统计量的分布，在使用第二种算法时，需要提供对应的统计量，即Effect, 以逻辑回归/线性回归为例，Effect对应的就是回归系数BETA, 标准误对应的就是回归系数的SE。
 
每种算法要求的gwas分析结果的格式稍有不同，其中以下3列是必须有的：
 1) SNP对应的id或者rs号
 2) test allele
 3) other allele
 
在关联分析的结果中，会有OR值来表征关联强弱，而OR值是一个比值，分子除以分母，分子对应的allele为test allele, 分母对应的allele为other allele。
 
2.1 基于pvalue的算法，额外要求以下3列
 1) Pvalue
 2) 效应方向：表示test allele和疾病关联方向的列，有正相关和负相关两种，以OR值为例，大于1为危险因素，小于1为保护因素，为了能够区分正负，OR值需要取log。
 3) 可选的列：表示样本的大小，根据每个数据集的样本大小来进行加权
 
2.2 基于标准误的算法，额外要求以下2列
 1) effect
 2) standard error
 
前文已经给过解释，effect对应回归分析中的回归系数beta值，standard error对应回归系数的SE。
 
在配置文件中，需要指定每个study的GWAS结果中上述列对应的标题，以及文件分隔符等选项，这样才能保证软件正确的识别所需的信息，一个配置文件 (文件名：metal.config.txt ) 的示例如下：
 
# (1) 第一个study的GWAS结果
# MARKER  指定SNP ID对应的列标题
# ALLELE  指定test allele和other allele对应的列标题
# PVALUE  指定P值对应的列标题
# EFFECT  指定效应方向对应的列标题
# WEIGHT  指定样本量对应的列标题
# PROCESS 指定gwas结果对应的文件
MARKER   SNP
ALLELE   RefAllele NonRefAllele
PVALUE   P-value
EFFECT   Effect
WEIGHT   N
PROCESS  inputfile1.txt
# (2) 第2个study的gwas结果
# 当格式和前一个格式相同时，可以省略不写
PROCESS  inputfile2.txt
# (3) 指定输出结果
# OUTFILE 指定输出结果的名称（注意：此为空格分隔的前缀和后缀）
# ANALYZE 表示进行分析
OUTFILE METAANALYSIS .TBL
ANALYZE
注：可以对inputfile的所有maker使用相同的权重 (在这种情况下，可以使用DEFAULTWEIGHT命令设置固定的权重)。WEIGHTLABEL命令优先于DEFAULTWEIGHT命令，因此使用中的WEIGHT列标签不能与inputfile中的任何列匹配。（WEIGHTLABEL可简写为WEIGHT）
 
WEIGHTLABEL     DONTUSECOLUMN
DEFAULTWEIGHT   1000
3. meta分析
 
配置好之后，只需执行以下命令即可进行分析：
 
# Linux 
metal metal.config.txt
meta分析后会生成两个文件，分别是 METAANALYSIS1.TBL 和 METAANALYSIS1.TBL.info
 
以METAL软件自带的测试数据进行说明：
 
 用测试数据进行meta分析：
 
# 软件示例文件如下：
cd generic-metal/examples/GlucoseExample/ 
ls -l
# DGI_three_regions.txt    # inputfile1
# MAGIC_FUSION_Results.txt.gz   # inputfile2
# magic_SARDINIA.tbl  # inputfile3
# metal.txt  # 配置文件
meta metal.txt
测试数据及结果文件：
 
 1) METAANALYSIS1.TBL 是meta分析的结果文件：通过输出结果中的pvalue,来筛选显著关联的位点，通过gwas meta-analysis, 可以达到增加样本量，提高检验效能的目的，而且有助于发现新的关联位点。
 
 结果文件理解：
 
① P-value值提示marker与性状的关联强度证据：此例中，rs560887 (靠近G6PC2) 在该区域显示出最强有力的关联证据。
 
② Direction显示不同研究中效应值的方向情况：几乎所有排在前10位的marker都在总共6806个个体中进行了核查，并在研究中表现出了一致的效应方向 ( 通过加号(+)和减号(-)的模式来说明)。例外的是rs853789，该marker在DGI的研究中没有进行检测，因此相应的权重较小，在效应列方向的第一项中有一个问号(?)。当分析大量的样本时，效应的方向可以让你一目了然地检查，是否所有的研究都支持你最强有力的发现。如果你发现有一项研究对你的最强信号的影响方向一直是相反的，这可能意味着该研究对等位基因1和2的标记不一致。
 
2) METAANALYSIS1.TBL.info 是meta分析的说明文档，比如 Marker 指的是什么。
 
 
4. 其他额外的分析选项
 
4.1 选择分析方案
 
SCHEME SAMPLESIZE     # 默认方法：使用p值和效应方向，根据样本大小加权
SCHEME STDERR         # 经典方法：使用效应大小估计和标准误差，根据标准误的倒数加权
STDERR SE             # 指定标准错误列的标签
注：虽然基于标准误差的权重在生物统计学文献中更常见，但需确保效应大小估计(beta系数)和标准误差在所有研究中使用相同的单位 (即确保在每个研究中使用了完全相同的性状，并应用了相同的转换)。不同研究中测量单位的使用不一致是这两种分析策略之间差异的最常见原因。
 
4.2 Genomic Control (GC) 校正
 
GENOMICCONTROL OFF   # 默认方式：即不校正
GENOMICCONTROL ON    # 自动校正检验统计量，以解释少量的人口分层或未解释的亲缘关系
GENOMICCONTROL [value]  # 使用指定的膨胀系数校正检验统计量
METAL可以对所有输入文件应用GC校正。通过将中位数检验统计量与随机期望的统计量进行比较来估计检验统计量的膨胀，然后将GC校正应用于p值 (SAMPLESIZE加权meta分析) 或标准误差 (STDERR加权meta分析)。这应仅适用于具有整个基因组数据的文件(即不应用于结果仅可用于候选位点或为GWAS结果的后续选择的少量SNP的设置)。
 
4.3 样本重叠校正
 
样本量加权元分析中的样本重叠校正(由Sebanti Sengupta开发，Daniel Taliun执行)。首先，METAL根据每个研究的z统计数据估计了两个或两个以上研究中常见的个体数量。然后，在计算总体z统计量时，METAL通过用估计的共有个体数量校正权重来调整样本重叠。为了校正样本大小加权meta分析中的样本重叠，可使用OVERLAP ON命令(仅对SCHEME SAMPLESIZE有效)。默认情况下，METAL使用Z-statistics <1来估计研究中常见的个体数量。使用ZCUTOFF [number]命令更改该阈值。
 
4.4 Strand信息
 
USESTRAND   ON
STRANDLABEL StrandColumnHeading
输入文件可以包含一个列，该列指示等位基因在哪个链上编码(给定为+/-)。如果出现此列，应该使用USESTRAND ON命令，并使用STRANDLABEL命令指定适当的列名。如果USESTRAND是关闭的，则假定该链对所有snp都是“+”的，尽管明显的链问题被METAL识别并适当处理(例如，当一项研究提供了A/G等位基因，而另一项研究提供了C/T等位基因)。
 
4.5 过滤
 
自定义过滤器可以用来选择包含在meta分析中的snp。例如，这可以用于在特定的最小等位基因频率范围内选择snp进行分析。
 
ADDFILTER N > 1000
ADDFILTER MAF > 0.01
这两个过滤条件表示：只考虑N列中的值大于1000且MAF列中的值大于0.01的条目。
 
可以使用<、>、<=、>=、=、!=和IN操作符定义过滤器。例如，要将分析限制为三个有趣的SNPs，可以使用(注意SNPs列表中没有空格)：
 
ADDFILTER MARKER_ID IN (rs1234,rs123456,rs123)
要删除之前定义的所有过滤器，使用命令：
 
REMOVEFILTERS
4.6 详细模式
 
VERBOSE ON
# To restric meta-analysis to two previously reported SNPs and summarize study specific results, uncomment the two lines that follow.
# ADDFILTER SNP IN (rs10830963,rs563694)
# VERBOSE ON
METAL允许完整输出所有输入文件中所有snp的单个汇总统计信息。这可能会创建一个非常大的文件，应该谨慎使用。通常，在使用此选项之前，应该创建自定义过滤器来将分析限制为感兴趣的snp。这个选项可以在许多研究中比较效应的方向，因为METAL考虑了所有的链翻转，并提供了相对于同一等位基因的效果方向。这也是一种双重检查METAL是否正确使用了预期数据的方法。
 
4.7 宽松模式
 
COLUMNCOUNTING STRICT   # 要求每行都达到预期的列数，即无缺失
COLUMNCOUNTING LENIENT  # 尝试解释列少于预期的行
默认情况下，METAL将跳过每个输入文件中没有预期列数的行。这通常是一个好主意，因为它可以避免在缺少列时产生不正确的结果。有时(例如，当每行末尾有可选的额外列时)，COLUMNCOUNTING LENIENT选项可能会很有用。
 
4.8 追踪等位基因频率
 
AVERAGEFREQ ON
MINMAXFREQ ON
METAL可以选择性地跟踪所有文件的效应等位基因频率，并报告平均、最小和最大效应等位基因频率。在METAL完成所有的链比对后，这些可以非常有用地检查不同队列的等位基因频率是否相似。不同研究中等位基因频率的巨大差异可能表明不同研究中参考等位基因的命名不一致。当这个特性被打开时，METAL要求所有的输入文件都有一个等位基因频率列。要使用FREQLABEL命令，指定等位基因频率信息的列名。
 
4.9 输入文件推荐
 
强烈建议对输入文件的所有SNP都提供两个等位基因的标签，即效应等位和非效应等位。只要每个输入文件的等位基因列都给出了，METAL就会适当地考虑不同输入文件使用不同参考等位基因的情况。等位基因可以用数字编码(A=1,C=2,G=3,T=4)或字母编码(A,C,G,T, a,c,g,t)，如果没有A/T或C/G SNP （互补配对的等位，ambiguous alleles），等位基因可以在任何一条链上。如果没有A/T或C/G SNP，等位基因可以在任何一条链上。对于A/T或C/G SNPs, METAL要求SNPs在不同输入文件中的一致链上，以便结果可解释。对于其他snp, METAL可以自动识别和解决链不一致。
 
P值：如果出现<0、>1或非数值的情况，会被认为是缺失值，并给出警告。
 
EFFECT列可以有正值和负值(例如，来自回归的beta值)，或者只是相对于参考等位基因的效应方向，列为“+”和“-”。例如，参考等位基因A(或效应等位基因A)的“+”效应(或任何正数)代表了等位基因A的拷贝数量的增加与性状值的增加相关的情况。对于离散性状，报告比值比通常都是正的。在这种情况下，为了计算效应的方向，应该查看比值比的对数。如果指定了EFFECT log(ODDS_RATIO_COLUMN)， METAL可以计算比值比。
 
要执行基于OR的meta分析，请在脚本的开头选择SCHEME STDERR。然后，对于每个文件，提供优势比的自然对数作为EFFECT列或另一个适当的统计数据(例如逻辑回归分析中相应的回归系数)。
 
5. 其他配置文件示例
 
#THIS SCRIPT EXECUTES AN ANALYSIS OF EIGHT STUDIES
#THE RESULTS FOR EACH STUDY ARE STORED IN FILES Inputfile1.txt THROUGH Inputfile8.txt
# 1. LOAD THE FIRST EIGHT INPUT FILES
# UNCOMMENT THE NEXT LINE TO ENABLE GenomicControl CORRECTION
# GENOMICCONTROL ON
# === DESCRIBE AND PROCESS THE FIRST INPUT FILE ===
MARKER SNP
ALLELE REF_ALLELE OTHER_ALLELE
EFFECT BETA
PVALUE PVALUE 
WEIGHT N
PROCESS inputfile1.txt
# === THE SECOND INPUT FILE HAS THE SAME FORMAT AND CAN BE PROCESSED IMMEDIATELY ===
PROCESS inputfile2.txt
# === DESCRIBE AND PROCESS THE THIRD INPUT FILE ===
MARKER SNP
ALLELE A_REF OTHER_ALLELE
EFFECT BETA
PVALUE pvalue 
WEIGHT N
PROCESS inputfile3.txt
# === DESCRIBE AND PROCESS THE FOURTH INPUT FILE ===
MARKER MARKERNAME
ALLELE EFFECTALLELE NON_EFFECT_ALLELE
EFFECT EFFECT1
PVALUE PVALUE
WEIGHT NONMISS
PROCESS inputfile4.txt 
# === CARRY OUT AN INTERIM ANALYSIS OF THE FIRST FOUR FILES ===
OUTFILE METAANALYSIS_inputfile1to4_ .tbl
ANALYZE 
# 2. LOAD THE NEXT FOUR INPUT FILES
# === DESCRIBE AND PROCESS THE FIFTH INPUT FILE ===
MARKER rsid
ALLELE EFFECT_ALLELE OTHER_ALLELE
EFFECT BETA
PVALUE Add_p
WEIGHT total_N
SEPARATOR COMMAS
PROCESS inputfile5.txt
# === THE SIXTH INPUT FILE HAS THE SAME FORMAT AND CAN BE PROCESSED IMMEDIATELY ===
PROCESS inputfile6.txt
# === DESCRIBE AND PROCESS THE SEVENTH INPUT FILE ===
ALLELE ALLELE OTHER_ALLELE
MARKER SNP
EFFECT BETA
PVALUE PVALUE
WEIGHT N
SEPARATOR WHITESPACE
PROCESS inputfile7.txt
# === DESCRIBE AND PROCESS THE EIGHTH INPUT FILE ===
ALLELE BETA_ALLELE OTHER_ALLELE
MARKER SNP
EFFECT BETA
PVALUE P_VAL
WEIGHT N
PROCESS inputfile8.txt 
#for the final meta-analysis of all 8 samples only output results if the
#combined weight is greater than 10000 people
OUTFILE METAANALYSIS_inputfile1-8_ .tbl
MINWEIGHT 10000
ANALYZE 
https://cloud.tencent.com/developer/article/1554734
 https://www.cnblogs.com/chenwenyan/p/10912521.html
 https://genome.sph.umich.edu/wiki/METAL_Documentation
 https://genome.sph.umich.edu/wiki/METAL_Quick_Start
群体分层是GWAS研究中一个比较常见的假阳性来源.
也就是说，如果数据存在群体分层，却不加以控制，那么很容易得到一堆假阳性位点。
当群体出现分层时，常规手段就是将分层的群体独立分析，最后再做meta分析。
1.如何判断群体是否分层
先用plink计算PCA，具体方法详见链接：GWAS群体分层 (Population stratification)：利用plink对基因型进...
				
Red 是一门新的编程语言，它受到了 REBOL 很大的启发，但由于它有本地代码编译器，Red 的应用领域更加广泛——下到系统编程上到高级脚本，同时提供了对现代的多核 CPU 并发编程的支持。
主要特点为：
函数式、声明式、符号式编程
支持基于原型的对象
同像性 （Red 是它自身的元语言）
类型声明可选，有丰富的数据类型（50+）
支持静态编译或 JIT 编译成本地代码
强大的并发和并行编程支持（actors、并行集合）
以内建的 Red/System DSL 提供底层系统编程的能力
支持高级脚本和 REPL 控制台
高度可嵌入
内存占用少，带有垃圾回收
磁盘占用少（1MB）
				
GWAS工具GAPIT最新版
本篇笔记主要内容是GWAS分析软件GAPIT最新版的安装和使用教程，包括常见的报错以及解决方案，主要出错位置在LDheatmap、stringi、nloptr、lme4等，测试安装的环境是东方天意的ECS云服务器（Linux centos7），R版本为4.2.2，虚拟环境使用conda。
什么是GWAS分析？
GWAS，即基因组关联分析（Genome-wide Association Study），是一种广泛应用于生物医学研究中的遗传学方法。其主要目的是在全基因组水平上，寻
 # get source
git clone --recurse-submodules git@github.com:MRCIEU/gwas2vcfweb.git
cd gwas2vcfweb
# configure host
# # store cromwell outputs here
mkdir -p /data/gwas2vcfweb
# # upload & download from here
mkdir -p /data/gwas2vcfweb/data
# # DB files
# ## Reference FASTA
mkdir -p /data/reference_genomes
cd /data/reference_genomes
wget ftp://gsapubftp-anonymous@ft
study1 = pd.read_csv("study1.csv")
study2 = pd.read_csv("study2.csv")
study3 = pd.read_csv("study3.csv")
# 将每个研究的p值进行变换，转化为z值
study1["z"] = np.sqrt(2) * sm.stats.proportion.proportions_ztest(study1["n_cases"], study1["n_total"], value=study1["OR"])[0]
study2["z"] = np.sqrt(2) * sm.stats.proportion.proportions_ztest(study2["n_cases"], study2["n_total"], value=study2["OR"])[0]
study3["z"] = np.sqrt(2) * sm.stats.proportion.proportions_ztest(study3["n_cases"], study3["n_total"], value=study3["OR"])[0]
# 将每个研究的z值和样本量进行合并
meta_data = pd.concat([study1[["SNP", "z", "n_cases", "n_total"]], study2[["SNP", "z", "n_cases", "n_total"]], study3[["SNP", "z", "n_cases", "n_total"]]])
# 计算每个SNP的z值和加权样本量
meta_data["wz"] = meta_data["z"] * np.sqrt(meta_data["n_cases"] + meta_data["n_controls"])
meta_data["w"] = np.sqrt(meta_data["n_cases"] + meta_data["n_controls"])
# 计算meta分析的z值和p值
meta_z = meta_data["wz"].sum() / meta_data["w"].sum()
meta_p = 2 * (1 - sm.stats.norm.cdf(abs(meta_z)))
# 输出meta分析结果
print("Meta-analysis result:")
print("z value: ", meta_z)
print("p value: ", meta_p)
这个脚本首先读取每个研究的GWAS结果文件，将每个研究的p值转化为z值，然后将每个研究的z值和样本量进行合并，计算每个SNP的加权z值和加权样本量，最后计算meta分析的z值和p值，输出结果。需要注意的是，不同研究的样本量、OR值等参数可能存在差异，需要根据具体情况进行调整。