基因组复杂性状分析工具"/>
GCTA：全基因组复杂性状分析工具
这是我个人的流程笔记，大概类似这种，主要是怕忘记，最好的解决办法就写笔记。
一.数据获取
plink
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
GCTA分析
2.0首先我们得了解一下基本概念，
（自己总结，可能有误）原文
什么是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 属于第一种算法，总结下面这篇文章，GCTA: A Tool for Genome-wide Complex Trait Analysis - PMC (nih.gov)
与单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）。
我们可以通过利用REML（ restricted maximum likelihood ）算法来无偏地估计全部SNP解释的方差 σ2g （继而估计SNP遗传力）。依赖于从所有SNP估计的GRM。
现在用的最多的是以上内容，有兴趣的可以阅读全文。还有一些关于，他们在五个功能域中开发了GCTA：数据管理，从一组SNP估计GRM，估计单个染色体或整个基因组上所有SNP解释的方差，估计连锁不平衡（LD）结构和模拟。GCTA: A Tool for Genome-wide Complex Trait Analysis - PMC --- GCTA：全基因组复杂性状分析工具 - PMC (nih.gov)
很好可以总结一下，就是我们得先估计GRM，如果不去掉有亲缘关系的个体，GRM的估计就会因为其基因组的相似性而引起偏差，最好去除有亲缘关系的个体，并添加PC作为协变量。进行
除此之外，还得了解一下，
随机效应：它通常被称为主观变化或主观错误。例如，不同人对同一主题的看法肯定会有所不同（不同），意见是主观的。
MLM和LMM：都是混合线性模型的一种。在MLM中，随机效应是从一个有限的总体中随机抽取的，而在LMM中，随机效应是从一个无限的总体中随机抽取的。总的来说，LMM比MLM更加灵活。
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 --grm test --grm-cutoff 0.025 --make-grm --out test_rm025
稀疏矩阵稀疏性处理
将稠密矩阵转换为稀疏矩阵形式几乎总是可以提升处理时间上的效率。
通常，为了处理稀疏性矩阵，我们会构造更为有效的数据结构，压缩稀疏行和列，或者通过PCA、SVD等方法来进行降维。
2.1.1 估计GRM
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可以去除掉有亲缘关系的个体
gcta64 -grm wenjian --grm-cutoff 0.025 --make-grm --out test_rm025
第二步：利用GCTA-GREML估计SNP遗传力
| 1 |
|
- –pheno 表型文件
- –reml 利用reml算法计算snp遗传力
计算结果储存在 test.hsq的文件中，当然我们以可以加入协变量，例如前10个主成分 PCA：
| 1 |
|
gcta64 --bfile wenjain --make-grm --out wenjain1 --thread-num 9#计算GRM
gcta64 --bfile cohort.filtered.noIndel.geno0.02.mind0.02.maf0.05.hwe1e5.het2.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软件的学习 - 知乎 (zhihu)
有一些格式介绍
一个是pca.txt，就是刚刚的文件
一个是对应ID表型文件pheno.txt，表型文件
二.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 /
conda config --add channels /
conda config --add channels /
conda config --add channels /
conda config --add channels /
conda config --add channels /
conda config --set show_channel_urls yes原文链接：
之后我们还是得先创造环境，基础环境一般是不下载软件的
conda create -n gcta
conda activate gcta
conda install -c bioconda gcta
gcta64
三.GCTA分析
bed就是我们一开始的VCF转出的二进制文件，或者即上面指定得到的文件
然后可以画图了。
四.qqman和manhattan的可视化
在可视化之前我们还得整理一下刚刚得到的文件mlm_wenjain_gwas，用awk提取得到这种格式就可以了，NR-1的意思就是从0开始的自然数依次递增要用awk提取的文件格式CHR SNP BP P
awk '{print $1,$2,NR-1,$10}' mlm_P298.weight_gwas.fastGWA > P298weight.txt;
因为我们数据比较多，所以我们只能用qqman来画，电脑带不动啊啊啊。emm
setwd("lunjin")#切到自己的路径
install.packages("qqman")#安装library("qqman")
vignette('qqman') #可以查看，RStudio里面qqman和manhattan自带的help，说的很详细，里面有自带的文件数据应该如何处理的格式，非常方便
results_assoc <- read.table("wenjain", header = T)#先QQ判断性状关联是否被基因组控制
z = qnorm(results_assoc$P/2) #将p值的一半转为正态分布数
lambda = round(median(z^2, na.rm = T)/0.454, 3) #计算lambda值，即基因组膨胀系数，λ越偏离1，就表明存在群体分层的现象越严重，容易有假阳性结果
qq(results_assoc$P, main = "Q-Q plot of GWAS p-values",
xlim = c(0, 7), ylim = c(0,8),
pch = 18, #pch参数用于设置散点的形状
col = "blue4",
cex = 0.5, #cex参数用于设置散点的大小
las = 1)#las参数用于设置标签的方向#当时这两种代码都能跑
manhattan(results_assoc, main = "Manhattan Plot", ylim = c(0, 10), cex = 0.6,
cex.axis = 0.9,
col = c("blue4", "orange3"),
suggestiveline = -log10(1e-05),
genomewideline = -log10(5e-08))manhattan(results_assoc, main = "Manhattan Plot",
ylim = c(0,10), cex = 0.6, #reduce the point size to 60%点
cex.axis = 0.9,# reduce the font size of the axis labels to 90%字
col = c("blue4", "orange3"),
annotatePval = NULL,
annotateTop = T)
全部流程大概是这样的了，在把图片导出来就可以了。
更多推荐
GCTA：全基因组复杂性状分析工具
发布评论