GCTA：全基因组复杂性状分析工具，1.5版本

BioXplorer

已于 2023-12-11 23:46:32 修改

阅读量1.5k

点赞数 3

文章标签：笔记

于 2023-10-11 22:46:06 首次发布

本文链接：https://blog.csdn.net/2302_79242191/article/details/133779454

版权

这是我个人的流程笔记，大概类似这种，怕忘记，最好的解决办法就写笔记。

一.PCA

1.0PCA的基本理解

PCA通过计算特征向量和特征值来进行降维。它首先对数据进行标准化处理，然后计算协方差矩阵或相关系数矩阵。接下来，通过对协方差矩阵进行特征值分解或奇异值分解，找到数据中最主要的特征向量（主成分）

PCA（主成分分析）的理解与应用 - 知乎 (zhihu.com)

【机器学习】降维——PCA（非常详细） - 知乎 (zhihu.com)

我真的才懂什么是 PCA (主成分分析) - 知乎 (zhihu.com)

PCA方法校正群体结构，GWAS该用多少个主成分？ - 知乎 (zhihu.com)

这些文章大家都可以看看，最推荐的是群体进化-GWAS分析 - 简书 (jianshu.com)

这里转载一部分内容，

PCA

分析原理
PCA(Principal Component Analysis)，即主成分分析方法，是一种使用最广泛的数据降维算法。PCA的主要思想是将n维特征映射到k维上，这k维是全新的正交特征也被称为主成分，是在原有n维特征的基础上重新构造出来的k维特征。PCA的工作就是从原始的空间中顺序地找一组相互正交的坐标轴，新的坐标轴的选择与数据本身是密切相关的。其中，第一个新坐标轴选择是原始数据中方差最大的方向，第二个新坐标轴选取是与第一个坐标轴正交的平面中使得方差最大的，第三个轴是与第1,2个轴正交的平面中方差最大的。依次类推，可以得到n个这样的坐标轴。通过这种方式获得的新的坐标轴，我们发现，大部分方差都包含在前面k个坐标轴中，后面的坐标轴所含的方差几乎为0。于是，我们可以忽略余下的坐标轴，只保留前面k个含有绝大部分方差的坐标轴。事实上，这相当于只保留包含绝大部分方差的维度特征，而忽略包含方差几乎为0的特征维度，实现对数据特征的降维处理。
简洁点来讲现在有这样的数据，100个样品，2M标记，即是2000000X100的矩阵，那么就通过数学降维的方法简化到100X3甚至100X2乘（即PC1，PC2）
分析软件
GCTA
tassel
EIGENSTRAT
结果展示
PCA结果矩阵（特征向量）
GWAS_1  0.0295707   0.0174155   -0.0245656
GWAS_10 0.0212291   -0.0552983  -0.0280335
GWAS_100    -0.0645872  0.00456635  0.00588907
GWAS_101    -0.0779853  -0.0317529  0.0138288
GWAS_102    -0.0790227  -0.0295285  0.0147819
GWAS_105    -0.0845384  0.000685319 0.0108059
GWAS_108    -0.0779536  -0.00380985 0.0101755
GWAS_109    -0.0789908  -0.00534946 0.012742
GWAS_11 0.0152839   0.0185823   -0.0305629
GWAS_110    -0.080786   -0.00255263 0.0131448
    *   第一列样品名称，第二列PC1的值，第三列PC2的值，第四列PC3的值（也就是平时看到的结果图的横纵坐标来源）
*   PCA解释数据结果（特征值）
54.402
32.2402
25.6809
18.0063
13.7968
9.6096
9.46086
9.00158
8.16587
7.60115
    *   这个结果每一个值对应的维度的解释情况，行数与样品数量一致，第一行代表第一维，依次类推；每一行除以所有
https://www.jianshu.com/p/d5a86164e809

1.1.获取bed文件

之前文章里面还有SNP提取，染色体提取，这里再说一次，这是提取我们指定样本的

plink --bfile 文件 --keep ID.txt --make-bed --out 文件2

注意这个ID格式得和fam格式一致，

1.2.plink计算PCA

分析整体可以直接上这个代码

plink --bfile bed --pca 20 --allow-no-sex --out pca

–pca参数用于对数据做pca分析，后面的20是取前20个pca结果的意思，一般不需要取太多，因为最后是用最显著的第一列和第二列作图。

最终得到两个文件：.eigenval 和.eigenvec，我们用.eigenvec文件进行绘图。

需要用最显著的第一列和第二列作图，在linux里面用awk提取包含ID的前四列，多个的话，记得用shell命令很爽的。

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

有一些格式介绍

一个是pca.txt，就是刚刚的文件

一个是对应ID表型文件pheno.txt，表型文件

二，GCTA分析

2.0首先我们得了解一下基本概念，

在此之前，

随机效应：它通常被称为主观变化或主观错误。例如，不同人对同一主题的看法肯定会有所不同（不同），意见是主观的。

MLM和LMM：都是混合线性模型的一种。在MLM中，随机效应是从一个有限的总体中随机抽取的，而在LMM中，随机效应是从一个无限的总体中随机抽取的。总的来说，LMM比MLM更加灵活。

稀疏矩阵稀疏性处理：

将稠密矩阵转换为稀疏矩阵形式几乎总是可以提升处理时间上的效率。

通常，为了处理稀疏性矩阵，我们会构造更为有效的数据结构，压缩稀疏行和列，或者通过PCA、SVD等方法来进行降维。

GCTA: A Tool for Genome-wide Complex Trait Analysis - PMC (nih.gov)

什么是SNP遗传力?-CSDN博客

表型(P)，由基因型(G)和环境因素(E)共同控制，即P = G + E

heritability，遗传力，用来描述表型变异中遗传变异的比例。基因G所占的比例，通过方差来描述遗传变异和表型变异，则遗传力的公式如下（广义遗传力）

SNP heritability，即SNP遗传力，公式如下（1）

用SNP位点来表征样本的遗传变异，在描述SNP位点和表型的关联性时，采用加性模型，将表型y看做是多个位点效应相加的结果

则SNP遗传力，公式如下（2）

这里的SNP位点是属于一个集合的，是部分位点，具体哪些位点取决于两个因素：第一个是检测到的SNP位点数量，芯片，NGS不同平台检测到的位点数不同；第二个是估算SNP遗传力的算法，目前有以下两种算法

GREML(Genomic relatedness matrix REstricted Maximum Likehood)

LDSC(linkage disequilibrium score regression)

我们这里用的GCTA 属于第一种算法，REML（ restricted maximum likelihood ）算法来无偏地估计全部SNP解释的方差 σ2g （继而估计SNP遗传力）。

与单SNP关联检验相对，GCTA中的GREML（genome-based restricted maximum likelihood）方法使用线性混合模型（ linear mixed models ， LMMs），将全部SNP的效应作为随机效应拟合。

模型如下：

y是一个n × 1 的表型向量，n是样本大小；

β 是主成分分析（PCA）固定效应的向量(诸如性别、年龄、以及一个或多个特征）；

u 是SNP效应的向量，服从如下分布：

I 是一个 n × n 的单位矩阵,

ɛ 是残余效应的向量，服从以下分布：

W是标准化的基因型矩阵，其中第ijth 个元素为：

其中 xij 是第j个个体的第i个SNP的参考allele的拷贝数， pi该 allele 的频率。

如果我们定义 A = WW′/N 并且定义 σ2g 为全部SNP解释的方差 , 即 σ2g=Nσ2u, （N为SNP的数量），那么等式 1 将等效于：

其中g是一个n x 1的向量，表示各个个体的全部的遗传效应，服从 g∼N(0,Aσ2g)。A被解释为个体之间的遗传关系矩阵（GRM）。

可以总结一下，就是我们得先估计GRM。如果不去掉有亲缘关系的个体，GRM的估计就会因为其基因组的相似性而引起偏差，最好去除有亲缘关系的个体，并添加PCA作为协变量。进行分析

2.1.1 估计GRM

GRM（遗传亲缘关系矩阵）：在GWAS中亲缘关系矩阵主要用于校正遗传背景，极易形成假阳性（没病，但却检测结果说有病）。GWAS分析要求样本间是相互独立的，在质控阶段，会根据样本间的亲缘关系，剔除亲缘关系较近的样本。

        IBD (Identical by descent)，血缘同源。IBS (Identical by state)，状态同源。

        通常我们会讲，如果一段DNA在两个或者多个个体中均有一致的核苷酸序列，那么可以定义该DNA片段属性为状态同源，即IBS。

        在两个或者多个个体中的一个IBS片段，如果遗传自一个共同的祖先个体，且没有发生重组，那么该片段也是血缘同源的，即IBD。

        一段DNA，如果是IBD，那么肯定也是IBS；如果片段不是IBD，也有可能是IBS，因为不同个体中发生了相同的突变，或者重组没有改变片段等原因所致。

个体间的亲缘关系(relationship)，可以通过以下两种方法估计：

        通过系谱进行计算，利用处于IBD状态的等位基因估计个体间的亲缘关系（通过plink计算IBD距离是最经典的一种）

        通过分子标记，利用处于IBS状态的基因估计个体间共享染色体片段的百分比，估计个体间的相似性（对于两两样本间的亲缘关系，可以用一个矩阵表示，即genetic relationship matrix, 简称GRM）针对亲缘关系大的情况，我们进行过滤，比如设定阈值为0.125, 亲缘关系大于该阈值的样本间就需要剔除其中一个样本。GCTA采用迭代的方式进行剔除，保证剩余样本的个数最大化。代码如下

gcta64 --bfile wenjan --autosome --maf 0.01 --make -grm --out wenjian1 thread-num-10

输入为plink格式的bed/bim/fam文件，–autosome 只用常染色体上的SNP，–maf 0.01 过滤掉maf小于0.01的snp，–make-grm 计算GRM，–out 输出文件的前缀，–thread-num 10 要使用的线程数

计算完成后GRM会存储在以下文件中： wenjain.grm.bin, wenjain.grm.N.bin and wenjain.grm.id

2.1.2可以去除掉有亲缘关系的个体

大家可以去看这个文章一文搞定GCTA软件的学习 - 知乎 (zhihu.com)

gcta64 -grm wenjian --grm-cutoff 0.025 --make-grm --out test_rm025

第二步：利用GCTA-GREML估计SNP遗传力

1

gcta64 --grm test --pheno test.phen --reml --out test --thread-num 10

–pheno 表型文件
–reml 利用reml算法计算snp遗传力

计算结果储存在 test.hsq的文件中，当然我们以可以加入协变量，例如前10个主成分 PCA：

1

gcta64 --grm test --pheno test.phen --reml --qcovar PCA.txt --out test --thread-num 10

gcta64 --bfile wenjain --make-grm --out wenjain1 --thread-num 9#计算GRM

gcta64 --bfile c.G1 --maf 0.01 --make-grm --out wenjian1 --thread-num 9#计算GRM

gcta64 --grm wenjian1 --make-bK-sparse 0.05 --out wenjain_sp_grm
gcta64 --bfile wenian --grm-sparse wenjain_sp_grm --pheno pheno.txt --out mlm_wenjian_gwas --qcovar pca.txt --thread-num 9 --fastGWA-mlm

三.GCTA下载（linux）

弄了快一个星期了，市面上解说GCTA下载的办法的文章。反正解决不了我的问题，为了能在conda里面直接下载软件，我整顿了一下我的conda，查了一些东西

（base）：~$ conda config --add channels conda-forge
（base）：~$ conda config --add channels bioconda

来自 <linux通过conda安装gatk4 - 知乎>

# conda ustc源 ，一个一个的下载
# conda ustc源
conda config --add channels https://mirrors.ustc.edu.cn/anaconda/pkgs/main/
conda config --add channels https://mirrors.ustc.edu.cn/anaconda/pkgs/free/
conda config --add channels https://mirrors.ustc.edu.cn/anaconda/cloud/conda-forge/
conda config --add channels https://mirrors.ustc.edu.cn/anaconda/cloud/msys2/
conda config --add channels https://mirrors.ustc.edu.cn/anaconda/cloud/bioconda/
conda config --add channels https://mirrors.ustc.edu.cn/anaconda/cloud/menpo/
conda config --set show_channel_urls yes

原文链接：https://blog.csdn.net/weixin_44179120/article/details/125844391

之后我们还是得先创造环境，基础环境一般是不下载软件的