关联分析全程，理解+过程，2.0版

 参考文章：

[GWAS]plink数据格式的转化以及缺失质控 - 知乎 (zhihu.com)

GWAS学习笔记 1 -- 如何对基因型数据进行严格的质量控制，以保证后续GWAS的成功进行 - 知乎 (zhihu.com)

GWAS质控要点 - 简书 (jianshu.com)

参考文献：

基于SNP芯片的丫杈猪保种群体遗传结构研究 - 中国知网 (cnki.net) 

基于简化基因组测序的永登七山羊遗传多样性分析 - 中国知网 (cnki.net) 

华西牛胴体及原始分割肉块重量性状遗传参数估计与全基因组关联分析 - 中国知网 (cnki.net) 

一，质控的原因是： 

1. 样本质量控制：

   • DNA质量：低质量的DNA样本可能导致基因分型错误或缺失数据。因此，需要检查样本的DNA浓度、纯度和完整性。
   • 样本重复性：通过检查重复样本之间的一致性来评估样本的重复性和可靠性。

2. 基因型质量控制：

   • 缺失数据：检查每个基因位点的缺失率，如果某个位点缺失的频率过高，则可能会影响结果的可靠性。
   • 杂合率：杂合率较高的位点可能会导致基因型分析中的错误。因此，需要排除杂合率异常的位点。
   • 稳定性：检查基因型数据在不同实验重复中的稳定性，以确定是否存在技术误差或其他问题。

3. 人群结构控制：

   • 人群结构：人群结构的存在可能引入假阳性结果，因为某些基因变异可能与特定人群相关。因此，需要使用适当的统计方法来检测和纠正人群结构的影响。

4. 表型数据质量控制：

   • 异常值：检查表型数据中是否存在异常值，这些异常值可能是数据记录错误或实验操作问题导致的。
   • 缺失数据：排除缺失值较多的表型数据，以避免对结果的影响。
   • 数据一致性：确保表型数据的收集和记录方式一致，避免不一致性引入误差。

 主要来说，一个意思。

1. 数据准确性：质控可以排除由于实验误差、技术问题或样本问题引起的数据错误。通过验证基因型和表型数据的准确性，可以提高结果的可信度。

2. 数据一致性：质控可以确保使用的基因型和表型数据在整个研究中保持一致。一致的数据收集和处理方式有助于减少潜在的偏见和混乱因素，从而提高结果的可靠性和解释性。

3. 排除混淆因素：质控过程可以检测和纠正可能导致假阳性结果的混淆因素，例如人群结构、样本污染或技术偏倚等。排除这些混淆因素可以确保GWAS结果更加准确和可靠。

4. 可重复性：质控确保研究数据的一致性和准确性，从而增加了该研究的可重复性。其他研究团队可以使用相同的数据集进行验证或复制研究结果，从而增强科学研究的可靠性。

二，质控分析的过程

1，先对指定样本的提取，然后再质控，

plink --bfile c --keep G2.txt --make-bed --out xian1 --chr-set 24

其中，c,vcf的二进制文件，G2，FID ID文件

plink -bfile xian1 -recode -out xian1 --chr-set 24

#bim/bed/fam转ped/map

wc -l xian1.map xian1.ped

注意查看多少个位点，每次质控分析统计多少个位点，方便观察和记录

2.geno0.02

先过滤个体缺失率高于2%的SNP，排除个体中"缺失率"高的SNP，低质量和无信息的SNP，无信息的SNP是指在研究样本中没有多态性或变异的SNP。大概就是，每个人都有，没有存在的必要

plink --bfile xian1 --geno 0.02 --make-bed --out xian1.geno.0.02 --chr-set 24

#转化为map和ped的形式，查看一下过滤后的行数，之前的为：
plink -bfile xian1.geno.0.02 -recode -out xian1.geno.0.02 --chr-set 24
wc -l xian1.geno.0.02.map xian1.geno.0.02.ped

 3.mind0.02

过滤SNP缺失率高于2%的个体，排除基因型缺失率高的个体

plink --bfile xian1.geno.0.02 --mind 0.02 --make-bed --out xian1.geno.0.02.mind.0.02 --chr-set 24
plink -bfile xian1.geno.0.02.mind.0.02 -recode -out xian1.geno.0.02.mind.0.02 --chr-set 24
wc -l xian1.geno.0.02.mind.0.02.map xian1.geno.0.02.mind.0.02.ped

4.maf0.05

大概的意思是，

可能我们研究的内容是，常见的基因变异在人群中出现的频率相对较高。

通过去除MAF较低的基因变异，可以减少对罕见变异的关注，常见的基因变异更容易在人群中观察到，这样我们可以用较小的样本数量检测到其与性状的相关性。从而提高研究的效率和可靠性。
根据最小等位基因频率（MAF）去筛选，

最小等位基因频率怎么计算？某个一个位点有AA或AT或TT，那么就可以计算A的基因频率和T的基因频率，qA + qT = 1，这里谁比较小，谁就是最小等位基因频率，qA = 0.3，qT = 0.7，那么这个位点的MAF为0.3。之所以用这个过滤标准，是因为MAF如果非常小，大约少于0.02，那么意味着大部分位点都是相同的基因型，这些位点贡献的信息非常少，增加了假阳性。 MAF为0，那就是所有位点只有一种基因型，这些位点没有贡献信息，放在计算中增加计算量，没有意义，所以要根据MAF进行过滤。

（对于一个位点“AA AT TT”，其中A的频率为基因频率，AA为基因型频率。MAF直接是对基因频率进行筛选，而哈温平衡检验，则是根据基因型推断出理想的（AA，AT，TT）的分布，然后和实际观察的进行适合性检验，然后得到P值，根据P值进行筛选。即P值越小，说明该位点越不符合哈温平衡。）

plink -bfile xian1.geno.0.02.mind.0.02 -maf 0.05 -make-bed -out xian1.geno.0.02.mind.0.02.maf.0.05 --chr-set 24
plink -bfile xian1.geno.0.02.mind.0.02.maf.0.05 -recode -out xian1.geno.0.02.mind.0.02.maf.0.05 --chr-set 24
wc -l xian1.geno.0.02.mind.0.02.maf.0.05.map xian1.geno.0.02.mind.0.02.maf.0.05.ped

 5.hwe1e5

笔记 GWAS 操作流程2-4：哈温平衡检验_哈温平衡的p值-CSDN博客

了确保研究样本的基因频率符合遗传学原理。 

什么是哈温平衡？

在满足哈迪温伯格的情况下，AB出现的概率为：
                                                P(A)×P(B)

若出现了连锁不平衡，则AB出现的概率为：
                                                P(AB)

满足条件：①种群足够大；②种群个体间随机交配；③没有突变；④没有选择；⑤没有迁移；⑥没有遗传漂变

怎么做哈温平衡？

「卡方适合性检验！」，对于卡方检验的理解，可以去看一下生物统计的课本，或者是阅读一下这些文章

遗传平衡定律(又称Hardy Weinberg定律)与卡方检验 - 知乎 (zhihu.com)

 四格表统计中该用Fisher确切概率法还是卡方检验？ (360doc.com)

哈迪温伯格原理：揭开群体遗传学之谜 - (lambdageeks.com)

 简单来说，如果不符合哈温平衡，可能无法准确评估基因座的遗传特征：哈代平衡原理是基因座遗传学研究的基础，无法准确评估基因座的遗传特征，例如等位基因频率、遗传模式等。

plink -bfile xian1.geno.0.02.mind.0.02.maf.0.05 --hardy --chr-set 24
plink -bfile xian1.geno.0.02.mind.0.02.maf.0.05 -hwe 1e-5 -make-bed -out xian1.geno.0.02.mind.0.02.maf.0.05.hwe1e5 --chr-set 24
plink -bfile xian1.geno.0.02.mind.0.02.maf.0.05.hwe1e5 -recode -out xian1.geno.0.02.mind.0.02.maf.0.05.hwe1e5 --chr-set 24
wc -l xian1.geno.0.02.mind.0.02.maf.0.05.hwe1e5.map xian1.geno.0.02.mind.0.02.maf.0.05.hwe1e5.ped

 6.het

一般自然群体，基因型个体的杂合度过高或者过低，都不正常，我们需要根据杂合度进行过滤。偏差可能表明样品受到污染，近亲繁殖。过滤掉杂合，去除杂合子可以帮助我们更准确地识别纯合子状态下的等位基因，从而更好地分析基因型和表型之间的关系。

6.1#计算杂合度

plink -bfile xian1.geno.0.02.mind.0.02.maf.0.05.hwe1e5 -het -out R_check --chr-set 24

6.2杂合R计算

setwd("E:/GWAS/zhi")
het <- read.table("R_check.het",header=T)
het$HET_RATE = (het$N.NM. - het$O.HOM.)/het$N.NM.
het_fail = subset(het,(het$HET_RATE < mean(het$HET_RATE) - 3*sd(het$HET_RATE))|(het$HET_RATE > mean(het$HET_RATE) + 3*sd(het$HET_RATE)))
het_fail$HET_DST=(het_fail$HET_RATE - mean(het$HET_RATE))/sd(het$HET_RATE) 
write.table(het_fail, "fail-het-qc.txt",row.names=FALSE) 

6.3

sed 's/"//g' fail-het-qc.txt |awk '{print $1,$2}' > het_fail-ind.txt

6.4用remove去掉个体：

plink -bfile xian1.geno.0.02.mind.0.02.maf.0.05.hwe1e5 --remove het_fail-ind.txt -make-bed -out xian1.geno.0.02.mind.0.02.maf.0.05.hwe1e5.het --chr-set 24
plink -bfile xian1.geno.0.02.mind.0.02.maf.0.05.hwe1e5.het -recode -out xian1.geno.0.02.mind.0.02.maf.0.05.hwe1e5.het --chr-set 24
wc -l xian1.geno.0.02.mind.0.02.maf.0.05.hwe1e5.het.map xian1.geno.0.02.mind.0.02.maf.0.05.hwe1e5.het.ped

三，提取PCA

plink --bfile xian1.geno.0.02.mind.0.02.maf.0.05.hwe1e5.het --pca 20 --allow-no-sex --out pca --chr-set 24

awk '{print$1,$2,$3,$4}' pca.eigenvec > pca.txt

四，GCTA分析

在我开头给的文献里面应用，

遗传力估计 使用GCTA软件通过GBLUP模型计算所得 方差组分与遗传力结果如表4所示,分割肉块部分 中四肢及臀部区域分割肉块重属于中高遗传力性 状,遗传力范围在0.41~0.57之间;躯干部分分割 肉块重以及出栏重、胴体重、屠宰率、净肉率属于中 等遗传力性状,遗传力范围在0.19~0.39之间。

PCA和GCTA的原理需要看一下，我的另外一篇文章，这里不做多的介绍GCTA：全基因组复杂性状分析工具-CSDN博客

cta64 --bfile xian1.geno.0.02.mind.0.02.maf.0.05.hwe1e5.het --make-grm --out 1111 --thread-num 9
gcta64 --grm 1111 --make-bK-sparse 0.05 --out 1111grm
gcta64 --bfile xian1.geno.0.02.mind.0.02.maf.0.05.hwe1e5.het --grm-sparse 1111grm --pheno p.txt --out 1111gwas --qcovar pca.txt --thread-num 9 --fastGWA-mlm
awk '{print $1,$2,$3,$10}' 1111gwas.fastGWA > 1111gwas.fastGWA.txt