GWAS和GS的第20篇学习笔记，建议阅读原始英文文献（有些点看英文更简单），翻译能力有限做不到信达雅，数组基础水平有限，部分点希望与大家一起交流。

1 introduction

AWM: Association weight matrix 关联权重矩阵。这是一种利用多表型GWAS结果结合网络推理算法构建基因关联网络的新方法。

原作者开发AWM方法的目的有两个，其一是希望可以有效地处理众多GWAS结果；另外一个则是因为在畜牧生产中往往可以获得一系列的相关表型，希望这个方法可以捕获单性状GWAS可能遗漏的相关联的基因。

虽然一个单一的表型可能在一个特定的研究项目中有其主要的关注点，但根据多元分析的启示，原作者寻找一组与主要关注表型有联系的表型，可以增加统计分析的准确性。

AWM方法最初由Frotes等人在肉牛性成熟期的相关研究上被实施[1]。传统的GWAS往往通过一个单表型的GWAS结果根据P值确定最重要的候选基因来解释一个复杂的性状，但是这个结果解释的程度是有限的，因为复杂性状（如性成熟期，肥胖等）受众多基因及生理途径的影响。随后，Frotes等人利用该方法探究了第二个群体[2]以及两群体[3]的相关研究。后一个例子由于添加了RIF (Regulatory Impact Factors ) 的成分而更丰富，RIF是一种用来识别导致表型差异的调节分子的指标[4,5]。

在本章中，作者提供了一个关于如何构建和使用AWM的教程，并解释了分析过程中每一步的逻辑。作者也讨论了表型数量、SNP芯片密度、对比选择等问题对AWM的影响，从而推断因果调控的交互作用。

2 Materials

使用AWM方法必须有可靠的多表型GWAS结果。

• GWAS前的标准QC；
• GWAS结果：SNPs、表型、SNPs效应值、P value;
  例如：如果正在研究糖尿病相关课题，一系列的表型将包括 空腹血糖、胰岛素释放 和 胰岛素敏感性的估计、体重、体重指数BMI、腰围、血压、胆固醇水平、糖尿病持续时间等等;每一种表型的GWAS结果都是必需的。关于研究中的人群的细节是不需要的，也不需要对所讨论的表型分析的模型。

影响AWM的因素：

• AWM可以适应一系列GWAS方法，但影响GWAS结果的因素同样影响AWM的质量，比如样本量、QC标准和模型拟合优度等。
• 物种基因组注释的程度和准确性;在这方面，unmapes SNP 无法用于基因网络的生成，除非在SNP选择标准中为它们定义了一些特定的标准而保留了一些 unmapes SNP(在mehod一节中对此有更多介绍)。
• AWM方法选择一系列表型时，所研究的这些表型的数量和它们之间的相互关系是重要的。

AWM的本质是在跨表型的SNP-SNP之间共同关联的基础上生成基因网络。衡量这种关联强度的基本度量是相关系数。由于相关系数估计的标准误差与表型数量的平方根的倒数成正比，因此只有当表型数量大量存在时，相关性的绝对值才会看起来很小，为~0.5。我们建议至少有10种表现型(见note1)。在糖尿病的例子中，病例对照研究将是主要的表现型，而其他提到的表现型(体重指数、胰岛素敏感性等)将作为生理学相关的特征进行调查，以完成该方案。

简而言之，需要的材料是以下三个输入文件(括号内的文件名与Methods部分提供的R脚本中使用的文件名相同):

(1)带有SNP P-values (AWM_P)的文件; ##P值
(2)每个SNP (AWM_A)的加性效应估计文件;##效应值
(3)一份包含SNP列表的文件，包含CHR、pos以及SNP到最近的基因的距离bp(GenomeMap)。

3 Methods

接下来，我们将结合PCIT[6]和RIF等算法[4,5]，对如何构建和使用AWM进行“一步一步”的指导，推断基因的共同关联和调控网络。图1提供了AWM/PCIT策略开发所涉及的流程图的示意图。
RIF: Regulatory Impact Factors (RIF) analysis
PCIT: Partial Correlation and Information Theory (PCIT) analysis

Fig1 关联权矩阵流程图:(a) 20个SNPs(行排列)在6个表型(列排列)的关联由一个矩阵表示，每个格子是SNP计算的效应值。与表型显著关联的SNP效应值由阴影格子表示，箭头所指的是关键表型，在第四列中突出显示。(b)经过一些列筛选，选择与关键表型和/或大量表型相关的10个snp来生成原始矩阵的子集AWM。(c)基于该矩阵，使用RIF和PCIT算法进行联合分析并使用绘图软件，可以生成具有10个基因（每个SNP代表一个基因）20个显著连接的基因关联网络，10个基因中有4个被标注为转录因子(这里用一个限制性的“T”表示)。(d)进一步的分析得到最终的结果，原始的网络可以被进一步简化为两个主转录因子通过9个连接连接到它们自己和各自的靶标

AWM是一个实数矩阵，行表示基因，列表示表型(体重指数、胰岛素敏感性、糖尿病等)。AWM中的单个{i, j}元素对应第i个基因SNP与第j个表型的关联。AWM的建立需要以下六个步骤:

Step 1: Primary SNP Selection

选择与关键表型相关的SNP(在糖尿病的例子中，这将是存在或没有临床糖尿病，病例对照)。在这一阶段，SNP选择的统计显著性强弱可以保持在不太严格的水平上(如:raw P-values <0.05)，因为AWM-PCIT方法的力量来自于整合一组看似独立的动作，每个动作都有关联的显著性(见note2)。这一步对应的R脚本如下:

step 1

操作
在这里我以6期的体重为表型，最后一期为 KeyPhenotype，

setwd("F:/MyDriver2.0/HUAZHONG/AWM")#设置路径
AWM_P <- read.table("AWM_Pvalue.txt", header=T, row.names=1)
AWM_A <- read.table("AWM_AdditiveEffects.txt", header=T, row.names=1)
Step1SNPs <- subset(AWM_P, BW180 <= 0.05) #根据关键表型P value筛选SNP
Step1SNP_IDs <- as.list(rownames(Step1SNPs)) #记录 Step1筛选的ID

Step 2: Exploring the dependency among phenotypes

对于上述选择的SNP，记录与之相关的非关键表型的平均数量（A_P）。为了一致性，建议在计算给定SNP所关联的表型数量时，使用与之前相同的显著性(例如，p值<0.05)。A_P是存在于可用表型之间的依赖量的指标。显然，当表型相互高度相关时，与关键表型相关的SNP也可能与许多其他表型相关。这一步对应的R脚本如下:

Step 2

Step1SNPs <- Step1SNPs[1:5]
Step2Ap <- t(apply(Step1SNPs, 1, function(Step1SNPs){table(Step1SNPs<=0.05)}))#对于上述选择的SNP，记录其与非关键表型的显著情况
Step2Ap <- t(as.data.frame(Step2Ap))
Step2Ap <- as.data.frame(Step2Ap)
Ap <- mean(Step2Ap[,2])#提取每个SNP与非关键表型显著的数目
Step2SNP_IDs <- as.list(rownames(Step2Ap))

Step 3: Secondary SNP Selection

选择至少与表型相关的SNP。在这一步中，我们通过进一步捕获那些SNPs来适当地利用表型之间存在的依赖性，这些SNPs虽然与关键表型无关，但它们与其余表型的关联大于平均方式。此步骤对应的R脚本如下:

step 3

Step3_select <- function(AWM_P){ AWM_P<=0.05 }#定义函数，提取P value阈值下的所有SNP
Step3Ap <- t(apply(AWM_P, 1, Step3_select))#二分类标签
SignificantPhenotypes <- apply(Step3Ap, 1, sum)#计数TRUE
Step3Ap <-  subset(SignificantPhenotypes, SignificantPhenotypes >= Ap)#筛选大于Ap数目的SNP，即这些SNP均至少与Ap个表型显著相关；
Step3Ap <- as.data.frame(Step3Ap)
Step3SNP_IDs <- as.list(rownames(Step3Ap))
Step3SNPs <- merge(Step1SNP_IDs, Step3SNP_IDs)#两次筛选去重合并
Step3SNPs <- t(as.data.frame(Step3SNPs))
rownames(Step3SNPs) <- Step3SNPs[,1]

Step 4: Exploiting the genome map

从两次选择步骤中获取的snp，识别出那些(i) map 到编码区，(ii)距离已知基因位置非常接近，或(iii)距离任何编码区很远。这一步是我们第一次尝试将一些“基因组功能”引入整个过程的地方。条件(i)和(iii)是简单的，仅取决于用于GWAS的SNP芯片的密度和所考虑的物种基因组组装的质量。条件(iii)也可以通过包括unmaped snp来扩展。条件(ii)的含义更为微妙，因为它取决于连锁不平衡(LD)的程度和特定基因启动子区域的可能大小。启动子区域是[7]基因表达所必需的。如果一个特定基因的表达要被调节(即激活或抑制)，那么它的调节器必须与它的启动子区域结合。启动子区域的突变影响邻近基因的表达被称为具有顺式作用效应。这与反式作用形成了对比，反式作用指的是影响基因表达的突变，该基因位于更远的基因组上，甚至在不同的染色体上。在牛上，作者在第(ii)条件中使用2500bp作为SNP到基因的阈值距离，并使用unmapped SNP或距离超过1.5Mb定义条件(iii)。所有在2500bp到1.5Mb范围内的SNP被定义为far。这一步骤中，剩余的远端或unmapped SNP将从进一步的分析中丢弃，作者继续只从近端和非常远的SNPs中选择。

step 4

GenomeMap <- read.table("GenomeMapFile.txt", header=T, row.names=1)#距离文件需要自己整理，可使用annovar软件结果整理
SNPsofInterest <- GenomeMap[rownames(Step3SNPs),]#提取上面筛选的snp
Step4SNPclose <- subset(SNPsofInterest, bp <= 10000)#距离筛选
Step4SNPfar <- subset(SNPsofInterest, bp > 15000)
Step4SNPclose_IDs <- as.list(rownames(Step4SNPclose))
Step4SNPfar_IDs <- as.list(rownames(Step4SNPfar))

step 5:Mapping the “One SNP—One Gene” relationship

对那些代表多个SNP的基因施加一个最终的SNP选择标准。这适用于步骤4 (i)和(ii)条件下的SNP。在这些情况下，与表型数量最多的SNP被选为该基因的代表。如果每个基因仍然有一个以上的SNP，那么将选择所有表型中P value平均值最低的SNP。同样,对于拥有强LD和远离基因的区域内的SNP的SNP集合(步骤4的条件iii)建议选择一个SNP，作为LD block和基于所有表型的平均关联强度选择的代表,。
。此步骤对应的R脚本如下:

step 5

SignificantPhenotypes <- as.data.frame(SignificantPhenotypes)
Step5SNPclose <- as.data.frame(SignificantPhenotypes[rownames(Step4SNPclose),])#计数：每个close_SNP与表型显著相关的数目
rownames(Step5SNPclose) <- rownames(Step4SNPclose)
Step5Genes <- data.frame(Step4SNPclose$Gene, Step5SNPclose[,1] )
row.names(Step5Genes)<-rownames(Step5SNPclose)
colnames(Step5Genes) <- c("Gene_ID", "SignifPhenos")
Step5Sorted <-  Step5Genes[order(Step5Genes$SignifPhenos,decreasing = TRUE),] #排序：依据相关表型的数目
Step5Unique <- subset( Step5Sorted, !duplicated(Step5Sorted$Gene_ID)) #取唯一值，去重复
Unique_SNP_IDs <- as.list(rownames(Step5Unique))
Step5SNPs <- merge(Unique_SNP_IDs, Step4SNPfar_IDs)
Step5SNPs <- t(as.data.frame(Step5SNPs))
rownames(Step5SNPs) <- Step5SNPs[,1]

Step 5: Mapping the “One SNP—One Gene” relationship

使用GWAS的结果，在AWM的每个{i, j}细胞中放置第i个基因与第j个表型的关联。第i个基因与第j个表现型的关联是通过按照上述步骤(i)和(ii)条件下的SNP映射步骤分配给该基因的SNP的加性效应来测量的。对于条件(iii)中的SNP，只需用SNP代码而不是官方的基因符号标记AWM行。为了允许表型之间的直接比较，SNP加性效应应该标准化，将每个个体效应除以所讨论的第j列表型中所有SNP效应的标准差。为了避免选择偏差，SNP效应的标准差应该是从整个GWAS中使用的所有SNP中获得的，而不是从被选择包含在AWM中的集合中获得的。下面是与此步骤对应的R脚本

图片.png

Step6SNP_ID <- AWM_A[rownames(Step5SNPs),]#从效应文件提取筛选出来的snp
rownames(Step5Unique) <- Step5Unique[,1]
Step5Unique <- as.list(rownames(Step5Unique))
Step6SNPs <- merge(Step5Unique, Step4SNPfar_Ids)
Step6SNPs <- t(as.data.frame(Step6SNPs))
Step6SNPs <- t(Step6SNPs)
rownames(Step6SNPs) <- Step6SNPs[,1]
write.table(Step6SNP_ID,file="PCIT.txt",row.name= rownames(Step6SNPs))

AWM的表状结构使其成为大多数聚类分析软件程序的输入。为此，我们推荐PermutMatrix8一个用户友好的程序来执行数据集的图形化分析。它适用于多种类型的聚类分析和各种统计序列化方法。

聚类分析的结果中，列(表现型)聚类可能程度较小，除非有强大的先验知识存在；行(基因)聚集在一起的方式允许对AWM中包含的数据的质量进行初始验证，从而为后续AWM分析的结果更可信(见注释3)。

后续分析最重要的方面是生成具有预测调控能力的基因网络。然而，在建立基因网络之前，作者建议进行其他分析。对这些列进行操作，作者建议计算从列的成对相关中获得的相关性（普通相关性），并将这些值与传统定量遗传方法(如REML方法)获得的遗传相关性估计值进行比较。由于基于AWM的遗传相关估计可能受到限制，不能摆脱环境协方差的影响，这些估计应该显示出与基于多元REML方法的等效估计的巨大相似性。

对行进行操作，可以利用所有成对基因(和/或SNP)之间的相关性来重建基因网络。从本质上说，具有极端相关估计(接近-1.0或+1.0)的基因对将通过网络中的直接链接相互连接。这一点减少了需要识别哪些相关估计是显著的，进而在网络重建中建立一个连接。

由于人们对基因网络领域的巨大兴趣，加上计算数学、多元统计和计算机科学的丰富，已经提出了许多方法来生成基因网络。总的来说，这些方法可以分为两类:无监督和监督。无监督方法中，网络直接从数据中推断出来，而监督方法需要训练集的已知调节相互作用的额外信息。特别是无监督方法在算法上是多样化的，并且比监督方法的计算速度更快。在无监督方法中，我们推荐PCIT算法[6,9]，该算法使用信息论(IT)框架中的部分相关(PC)来确定基因对之间的哪些相关性是显著的。尽管如此，作者承认其他生成基因网络的方法也适用于AWM，这仍有待检验。

PCIT的一个吸引人的特点是与硬性阈值方法不同。硬性阈值方法将每个相关估计与单个阈值进行比较，在PCIT中，两个基因之间的每个相关是相对于这两个基因中的每一个和数据集中任何其他基因之间的相关进行比较。这一特征意味着，如果没有其他基因与任何一个绝对值大于0.4的基因存在相关性，那么两个基因(如|r| < 0.4)之间看似很小的相关性就可以认为是显著的。在另一个极端，PCIT认为非常接近+1或-1的估计相关系数是不显著的。

从PCIT的输出文件可以用作Cytoscape的输入文件(http://www.cytoscape.org) [10]，以可视化和分析产生的网络。许多Cytoscape插件可用来描述网络的各种特征。其中,我们建议MCODE [11] (http://baderlab.org/Software/MCODE)以识别高度相互关联的基因簇，以期与该簇中基因的已知功能和机制行为一致，并对AWM-PCIT衍生的关键表型具有生物学意义。

在具有机制意义的基因连接中，最突出的是转录因子(TF)与其靶基因之间的调节相互作用。AWM-PCIT基因网络中预测的这种性质连接可以通过数据库（by mining databases of experimentally validated TF–target interactions）来验证，如TRANSFAC (http://www.biobase-international.com/gene-regulation)、Genomatix (http://www.genomatix.de)。此外，从这些数据库中获得的信息可以用于正式地检查派生网络的置信度，当这个网络与控制网络比较时，控制网络可以在改变AWM的行和列之后建立。

4 Discussion

要考虑的最后一个主题是如何实现AWM-PCIT方法的结果。由于它与QTL检测或育种价值预测等标准方法没有直接关系，因此如何应用这些发现并不明显。对于这个问题，作者的回答是，锚定在已知生物学相关性基因上的SNP更有可能脱离生物学知识，因此在群体中具有相关性。事实上，其他作者已经表明，将基因功能和调控信息与GWAS的统计关联整合起来，可能会增强对影响复杂表型的基因组机制的认识[12,13]。第二个可能不太明显的答案是，通过跨越整个AWM-PCIT基因网络大部分拓扑结构的最小(希望很小)基因相互作用集可以揭示一组SNP，可以用来开发它们明显的上位性相互作用。这一价值主张和其他价值主张仍有待探讨。

5 Notes

1

在实现AWM时值得注意的一个问题是，与大量表现型“最低限度”相关的SNP将被包含在AWM中，而不管哪个表现型是重点或关键的表现型。我们将这种现象视为两种不一定是排他的情况的结果:要么受审查的一系列表型彼此高度相关，因此不能充分跨越关键表型背后的整个生物学，要么这些“总是相关的”SNPs对应的基因在广泛的生物过程中普遍重要，如主调节因子。虽然我们没有遇到这种第二种可能性的例子，但我们预计在处理组织分化和发育时，它会出现，细胞周期的主调节因子容易影响许多所研究的表型。无论哪种情况，我们建议继续正常处理整个AWM-PCIT，同时要谨慎关注这些SNP和它们所代表的基因。

2

很可能，AWM最受期待的批评之一就是处理多重测试的问题。由于在任何给定的GWAS中都有成千上万的SNPs被测试，谨慎地认为只有那些具有非常严格的统计学显著性(即P value非常小)的SNPs才应该被考虑。这一论点在控制FDR方面有其正式的理由。然而，它似乎也与选择包含在AWM中的SNP的看似宽松的标准相冲突。如前所述，SNPs选择的统计显著性强度可以保持在一个较低的严格水平(例如，raw P value<0.05)，因为AWM-PCIT方法的力量来自于一组独立行动的整合，其中在控制FDR的数量上每个行动都有关联的重要性。显然，与关键表型关联最小的SNP(即p值小于但接近0.05)，但不与AWM中包含的任何其他表现型关联，其在表现型间的关联将是a flat profile。在AWM-PCIT基因网络的生成过程中，一个平坦的基因图谱不可能与任何其他SNP的基因图谱产生协同作用。此外，在分析网络时，作者建议用户在不同的分辨率层次上依次可视化。

3

如果在最初的聚类分析中，已知的与基因相关的性状，甚至是与生理非常接近的性状，并没有相互接近，那么数据的质量就应该受到质疑。例如，在糖尿病的例子中，一个人应该关注的结果，如果身体质量指数和体重没有聚集在一起。在这种情况下，我们的建议是后退一步，如果可能的话，回顾GWAS的结果，包括探索GWAS模型中拟合的固定效应和协变量。

References

1. Fortes MRS et al (2010) Association weight matrix for the genetic dissection of puberty in beef cattle. Proc Natl Acad Sci USA 107:13642–13647
2. Fortes MRS et al (2010) A new method for exploring genome-wide associations applied to cattle puberty. 9th World Congress on genetics applied to livestock production. German Society for Animal Science, Leipzig, Germany, pp 4–166. ISBN 978-3-00-031608-1.
3. Fortes MRS et al (2011) A single nucleotide polymorphism-derived regulatory gene network underlying puberty in 2 tropical breeds of beef cattle. J Anim Sci 89:1669–1683
4. Hudson NJ, Reverter A, Dalrymple BP (2009) A differential wiring analysis of expression data correctly identifies the gene containing the causal mutation. PLoS Comput Biol 5: e1000382
5. Reverter A et al (2010) Regulatory impact factors: unraveling the transcriptional regulation of complex traits from expression data. Bioinformatics 26:896–904
6. Reverter A, Chan EKF (2008) Combining partial correlation and an information theory approach to the reversed engineering of gene co-expression networks. Bioinformatics 24:2491–2497
7. Cookson W et al (2009) Mapping complex disease traits with global gene expression. Nat Rev Genet 10:184–194
8. Caraux G, Pinloche S (2005) Permutmatrix: a graphical environment to arrange gene expression profiles in optimal linear order. Bioinformatics 21:1280–1281
9. Watson-Haigh NS, Kadarmideen HN, Reverter A (2010) PCIT: an R package for weighted gene co-expression networks based on partial correlation and information theory approaches. Bioinformatics 26:411–413
10. Shannon P et al (2003) Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res 13:2498–2504
11. Bader GD, Hogue CW (2003) An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinformatics 4:2
12. Snelling WM et al (2012) How SNP chips will advance our knowledge of factors controlling puberty and aid in selecting replacement beef females. J Anim Sci 90(4):1152–1165. doi:10.2527/jas.2011-4581
13. Weller JI, Ron M (2011) Invited review: quantitative trait nucleotide determination in the era of genomic selection. J Dairy Sci 94:1082–1090