一、Linux系统和生信分析软件安装
1、Linux系统安装配置
ubuntu20.04 版本；默认安装。
安装向日葵远程控制
sudo dpkg -i SunloginClient_15.2.0.63064_amd64.deb
sudo apt-get install -f
sudo dpkg -i SunloginClient_15.2.0.63064_amd64.deb
设置远程控制长期密码。
ubuntu24.04 版本；默认安装。
安装向日葵远程控制
wget http://archive.ubuntu.com/ubuntu/pool/universe/g/gconf/gconf2-common_3.2.6-7ubuntu2_all.deb
wget http://archive.ubuntu.com/ubuntu/pool/universe/g/gconf/libgconf-2-4_3.2.6-7ubuntu2_amd64.deb
sudo dpkg -i gconf2-common_3.2.6-7ubuntu2_all.deb
sudo dpkg -i libgconf-2-4_3.2.6-7ubuntu2_amd64.deb
sudo apt --fix-broken install
https://sunlogin.oray.com/download/linux?type=personal&ici=sunlogin_navigation
sudo dpkg -i SunloginClient_15.2.0.63064_amd64.deb
sudo apt --fix-broken install
#######
如果其它主机连接本机黑屏可禁用WaylandEnable
vim /etc/gdm3/custom.conf
grep WaylandEnable /etc/gdm3/custom.conf
WaylandEnable=false
:wq！
或者注销当前账户，在登陆时选择在登录界面： 在用户名下方或密码框旁边找到一个齿轮图标。 点击齿轮图标，选择 “Ubuntu on Xorg”
2、Win系统 下载并安装MobaXterm
https://download.mobatek.net/2502024121622306/MobaXterm_Installer_v25.0.zip
3、GWAS分析相关软件安装
[1] 下载Miniconda3 https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/
最新版Miniconda3-py39_24.11.1-0-Linux-x86_64.sh
#拷贝到Linux子系统的root目录或新建用户名的Home目录
[2] 安装Miniconda
bash Miniconda3-py39_24.11.1-0-Linux-x86_64.sh
配置conda
source ~/.bashrc
Conda升级 conda update -n base -c defaults conda
[3] 使用conda安装plink
conda install -c bioconda plink
[4] 使用conda安装bedtools
conda install -c bioconda bedtools
[5] 使用conda安装admixture
conda install -c bioconda admixture
[6] 使用conda安装gcta
conda install -c bioconda gcta
[7] 安装vcftools
conda install -c bioconda vcftools
[8] 安装gemma
conda install -c bioconda gemma
[8] 安装Emmax
https://csg.sph.umich.edu//kang/emmax/download/index.html
https://csg.sph.umich.edu//kang/emmax/download/emmax-intel-binary-20120210.tar.gz
[9] EIGENSTRAT
从https://github.com/DReichLab/EIG下载
直接拷贝安装软件文件夹
[10] 安装xpclr
conda create -n xpclr_env python=2.7
conda activate xpclr_env # 更换为xpclr运行环境，在使用XPCLR前，需要再次运行。使用结束后，输入conda deactivate，或输入source ~/.bashrc切换回来
conda install -c bioconda xpclr
[11] 使用conda安装bwa
conda install -c bioconda bwa
[12] 使用conda安装blast
conda install -c bioconda blast
二、GWAS分析命令及解释
1. 过滤VCF文件
使用plink过滤低质量和低频率的SNP并直接生成VCF文件
复制代码
plink --vcf test.vcf --maf 0.05 --geno 0.1 --recode vcf-iid --out test.maf0.05.int0.9 --allow-extra-chr
1. Plink：这是PLINK软件的命令行工具，用于处理基因数据。
2. --vcf test.vcf：指定输入文件为test.vcf。
3. --maf 0.05：设置MAF阈值为0.05。过滤掉MAF小于0.05的位点（SNPs）。
4. --geno 0.1：设置基因型缺失率的阈值为0.1。过滤掉缺失率大于10%的位点。
5. --recode vcf-iid：将过滤后的数据重新编码为VCF格式，并保留样本IID。
6. --out test.maf0.05.int0.9：指定输出文件的前缀为test.maf0.05.int0.9。
7. --allow-extra-chr：允许非人类标准染色体名称。
可选：将VCF数据转换为hapmap格式
plink --vcf test.vcf --maf 0.05 --geno 0.1 --recode hapmap --allow-extra-chr --out rice_filtered_hapmap
2. LD修剪（去除连锁不平衡）
复制代码
plink --vcf test.maf0.05.int0.9.vcf --indep-pairwise 50 10 0.2 --out test.maf0.05.int0.9 --allow-extra-chr
1. --indep-pairwise 50 10 0.2：窗口50 SNP，步长5，保留r² < 0.2的SNP。
2. --out test.maf0.05.int0.9：指定输出文件的前缀为test.maf0.05.int0.9
输出文件：test.maf0.05.int0.9.prune.in：包含保留的SNP ID。
test.maf0.05.int0.9.prune.out：包含被去除的SNP ID。
3. 提取筛选结果
复制代码
plink --vcf test.maf0.05.int0.9.vcf --make-bed --extract test.maf0.05.int0.9.prune.in --out test.maf0.05.int0.9.prune.in --allow-extra-chr
1. --make-bed：将输入的VCF文件转换为PLINK的BED格式文件。
2. --extract test.maf0.05.int0.9.prune.in：
指定一个包含要提取的SNP ID的文件test.maf0.05.int0.9.prune.in。
作用：这个文件通常是由前面的LD剪枝步骤生成的，包含保留的SNP ID。通过--extract选项，PLINK只会提取这些指定的SNP，从而生成一个更小的、只包含独立SNP的数据集。
3. --out test.maf0.05.int0.9.prune.in
指定输出文件的前缀为test.maf0.05.int0.9.prune.in。
PLINK会生成三个文件：
test.maf0.05.int0.9.prune.in.bed：基因型数据的二进制文件。
test.maf0.05.int0.9.prune.in.bim：包含SNP信息的文件（染色体、位置、SNP ID等）。
test.maf0.05.int0.9.prune.in.fam：包含样本信息的文件（样本ID、家族ID等）。
4. 将筛选后的数据转换为vcf格式
复制代码
plink --bfile test.maf0.05.int0.9.prune.in --recode vcf-iid --out test.maf0.05.int0.9.prune.in --allow-extra-chr
1. --bfile test.maf0.05.int0.9.prune.in：指定输入文件的前缀为test.maf0.05.int0.9.prune.in。PLINK的--bfile选项用于指定BED格式文件的前缀。PLINK会自动查找与该前缀匹配的三个文件：
test.maf0.05.int0.9.prune.in.bed：基因型数据的二进制文件。
test.maf0.05.int0.9.prune.in.bim：包含SNP信息的文件。
test.maf0.05.int0.9.prune.in.fam：包含样本信息的文件。
5. 将筛选后的文件转换为admixture数据格式用于群体结构分析
复制代码
plink -bfile test.maf0.05.int0.9.prune.in --recode 12 --out test.maf0.05.int0.9.prune.in --allow-extra-chr
1. --recode 12：将输入的BED文件重新编码为PED格式，并保留基因型数据的完整信息。
2. --recode：指定输出文件的格式。这里选择PED格式。12：表示输出的基因型数据将使用数字编码（1和2）来表示等位基因。具体编码规则如下：1：表示第一个等位基因。2：表示第二个等位基因。0：表示缺失基因型。这种编码方式便于后续的统计分析和其他软件的兼容性。
6. 将基因型文件转换为EMMAX格式
复制代码
plink --vcf test.maf0.05.int0.9.vcf --recode 12 transpose --output-missing-genotype 0 --out test.maf0.05.int0.9 --autosome-num 90 --allow-extra-chr
1. --recode：指定输出文件的格式。这里选择PED格式。
12：表示输出的基因型数据将使用数字编码（1和2）来表示等位基因。具体编码规则如下：
1：表示第一个等位基因。
2：表示第二个等位基因。
0：表示缺失基因型（由--output-missing-genotype 0指定）。
2. transpose：对输出的PED文件进行转置。在转置后的文件中，行表示SNP，列表示样本。这与常规的PED文件格式（行表示样本，列表示SNP）相反。转置后的文件更适合某些分析工具（如PCA分析）得到.tped和.tfam文件。
3. --output-missing-genotype 0：指定缺失基因型的编码方式为0。而不是默认的-9
7. 亲缘关系计算（亲缘关系矩阵）
复制代码
chmod +x emmax-kin 或者 Chmod +x emmax-kin-intel64
./emmax-kin test.maf0.05.int0.9 -v -h -d 10 -o test.maf0.05.int0.9.BN.kinf
## 如果你下载的emmax软件是intel这个版本的，这里需要注意，命令有点不同（通常使用下边的指令）
./emmax-kin-intel64 test.maf0.05.int0.9 -v -d 10 -o test.maf0.05.int0.9.BN.kinf
1. test.maf0.05.int0.9：输入文件的前缀。通常情况下，这个前缀对应于一个.tped文件和一个.tfam文件。
.tped：转置后的基因型文件，行表示SNP，列表示样本。
.tfam：样本信息文件，包含样本的家族ID（FID）和个体ID（IID）。
2. -v：启用详细模式（Verbose Mode），在运行过程中输出更多的调试信息和进度信息。这有助于监控程序的运行状态和排查问题。
3. -d 10：指定计算亲缘关系矩阵时使用的参数。具体来说，-d 参数用于设置亲缘关系矩阵计算的精度（decimal places）。作用：在计算亲缘关系矩阵时，保留小数点后10位数字。这可以提高计算的精度，但可能会增加计算时间和存储需求。
4. -o test.maf0.05.int0.9.BN.kinf：指定输出文件的名称为test.maf0.05.int0.9.BN.kinf。作用：输出的亲缘关系矩阵文件将以指定的名称保存，格式为.kinf。
输出文件
5. test.maf0.05.int0.9.BN.kinf：这是一个亲缘关系矩阵文件，包含样本之间的亲缘关系系数。格式通常是一个对称矩阵，矩阵的行和列分别代表样本，矩阵中的每个元素表示两个样本之间的亲缘关系系数。
test.maf0.05.int0.9.BN.kinf 文件展示
R语言绘制亲缘关系热图（结果可视化）
为test.maf0.05.int0.9.BN.kinf 文件添加列名和行名，然后进行以下步骤（行名和列名均为品种名称。）
library(pheatmap)
kinship_matrix <- read.table("JSK.BN.kinf", header = TRUE, row.names = 1)
p1 <- pheatmap(kinship_matrix, #要绘制热图的矩阵
color = colorRampPalette(c('blue','white','red'))(100), #热图色块颜色是从蓝到红分为100个等级
border_color = "NA", #热图中每个色块的边框颜色，NA表示无边框
scale = "row", #按行进行归一化，"column"表示按列，"none"表示不进行归一化
cluster_rows = FALSE, #是否对行进行聚类
cluster_cols = FALSE, #是否对列进行聚类
legend = TRUE, #是否显示图例
legend_breaks = c(-1, 0, 1), #设置图例的断点
legend_labels = c("low","","heigh"), #设置图例断点处的标签
show_rownames = FALSE, #是否显示行名
show_colnames = FALSE, #是否显示列名
fontsize = 8 #字体大小，可以通过fontsize_row、fontsize_col参数分别设置行列名的字体大小
)
pdf("热图1.pdf",width = 12,height = 8)
p1
8. 群体结构文件 首先admixture群体结构分析（Q矩阵获取）
复制代码
for i in {1..15}
do
admixture --cv test.maf0.05.int0.9.prune.in.ped ${i} -j48 >> log.txt
done
1. for i in {1..15}：这是一个循环语句，i的值从1到15，每次循环i的值增加1。
作用：循环运行admixture命令，分别指定不同的群体数量（从1到15可以自行调节）。
2. do：开始循环体的标记。
admixture --cv test.maf0.05.int0.9.prune.in.ped ${i} -j48 >> log.txt：
这是循环体中的命令，用于运行admixture软件。
3. 参数解释：
4. admixture：调用admixture软件，用于群体结构分析。
5. --cv：启用交叉验证（Cross-Validation），用于评估不同群体数量（K值）的模型拟合优度。
6. test.maf0.05.int0.9.prune.in.ped：输入文件，通常是PLINK格式的PED文件，包含基因型数据。
7. ${i}：当前循环的群体数量（K值）。在每次循环中，i的值从1到15，分别表示不同的群体数量。
8. -j48：指定使用48个线程进行计算，加速处理过程。
9. >> log.txt：将每次运行的输出结果追加到log.txt文件中。这样可以记录每个K值的交叉验证结果，便于后续分析。
10. done结束循环体的标记。
11. 生成文件名：
test.maf0.05.int0.9.prune.in.1.Q
test.maf0.05.int0.9.prune.in.1.P
8.1查看最佳分群（选择CV error最小的的K对应的Q矩阵，水稻研究选择的K值通常不宜太大）
复制指令
grep CV log.txt
K=5时CV值最小，所以最佳分群为5 K=5时候的Q矩阵用于后续分析。
admixture产生的Q矩阵格式如下以K=5为案例(前6行)。
因为文件相对较小，手动添加表头信息，以及删除最后一列，保存5.Q用于后续分析,之后如下：
复制代码
Vi test.maf0.05.int0.9.prune.in.5.Q
i（进入编辑模式）
添加文件信息
ESC（退出编辑模式）
：wq！（保存并退出文件）
8.2表型数据（一次输入一个表型数据）
可以直接使用excel进行转化，然后导出成为txt格式即可
9. 协变量文件
R语言内完成
复制代码
tfam <- read.table("test.maf0.05.int0.9.tfam", header = F, stringsAsFactors = F)
admix <- read.table("14.Q", header = T, check.names = F, stringsAsFactors = F, skip = 1)
admix <- admix[match(tfam$V1, admix$`<Trait>`), ]
admix <- cbind(admix[,1], admix[,1], rep(1, nrow(admix)), admix[,-1])
write.table(admix, file = "test.emmax.cov.txt", col.names = F, row.names = F, sep = "\t", quote = F)
协变量文件格式
10. GWAS分析
复制代码
emmax-intel64 -v -d 10 -t test.maf0.05.int0.9 -o test_emmax_cov -p test_trait_emmax.txt –k test.maf0.05.int0.9.BN.kinf -c test.cov.emmax.txt
1. emmax-intel64：这是 EMMAX 程序的可执行文件，用于进行全基因组关联分析（GWAS）。
2. –v：打印详细信息，用于调试和记录运行过程。
3. -d 10：指定亲缘关系矩阵的维度（或缩放因子）。通常用于调整亲缘关系矩阵的计算精度。
4. -t test.maf0.05.int0.9：指定输入的基因型数据前缀。test.maf0.05.int0.9 是经过过滤的基因型数据文件前缀，通常为 .tped 和 .tfam 格式。
5. -o test_emmax_cov：指定输出文件的前缀。结果文件将以 test_emmax_cov 为前缀生成。
6. -p test_trait_emmax.txt：指定表型数据文件。表型数据文件需要与基因型数据的样本顺序一致。
7. -k test.maf0.05.int0.9.BN.kinf：指定亲缘关系矩阵文件。这是由 emmax-kin-intel64 生成的文件，用于校正群体结构和亲缘关系。
8. -c test.cov.emmax.txt：指定协变量文件。协变量文件用于校正环境因素或其他非遗传因素对表型的影响
11. 曼哈顿图及QQ图绘制
12. 注意事项
文件一致性：确保 .tped、.tfam 和表型文件中的样本顺序完全一致。