Emmax-pipeline

🏡Zhang Lab, Jun 04, 2023 🏡

🔬Multi-omics research center🔬

🔬State Key Laboratory of Crop Stress Adaption and Improvement, HENU, Kaifeng, Henan 🔬

环境配置和软件安装

本流程需要在linux/Unix(MacOS)或者windows WSL环境下配置。建议通过conda完成环境搭建及软件安装。下面将讲解详细步骤。

conda下搭建emmax-gwas分析环境

安装conda

首先到清华镜像站下载miniconda安装包：index of /annconda/miniconda，根据自己操作系统类型选择最新日期对应的安装包。

1.1 通过wget将Miniconda安装包下载到服务器； 1.2 添加执行权限，运行安装脚本，结束的时候注意不要选择随系统启动，以防conda污染环境变量
1.3 创建一个alias快捷启动conda
1.4 替换软件镜像

代码如下：

#| code-fold: true
#| code-summary: "Show the code"
#> 下载到电脑
wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-py39_23.3.1-0-Linux-x86_64.sh

#> --2023-06-14 09:14:37--  https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-py39_23.3.1-0-Linux-x86_64.sh
#> Resolving mirrors.tuna.tsinghua.edu.cn (mirrors.tuna.tsinghua.edu.cn)... 101.6.15.130
#> Connecting to mirrors.tuna.tsinghua.edu.cn (mirrors.tuna.tsinghua.edu.cn)|101.6.15.130|:443... connected.
#> HTTP request sent, awaiting response... 200 OK
#> Length: 63969016 (61M) [application/octet-stream]
#> Saving to: ‘Miniconda3-py39_23.3.1-0-Linux-x86_64.sh’

#> Miniconda3-py39_23.3.1-0-Linux-aa 100%[==========================================================>]  61.00M  53.7MB/s    in 1.1s 2023-06-14 09:14:38 (53.7 MB/s) - ‘
#> Miniconda3-py39_23.3.1-0-Linux-x86_64.sh’ saved [63969016/63969016]
#> 添加可执行权限
chmod +x Miniconda3-py39_23.3.1-0-Linux-x86_64.sh
#> 执行安装脚本
bash Miniconda3-py39_23.3.1-0-Linux-x86_64.sh
#> 下面一路yes+回车就可以了 就是在最后看到是否开机启动的时候选择no
#> Do you wish the installer to initialize Miniconda3
#> by running conda init? [yes|no]
#> [no] >>> no
#> 新建一个放置alias的配置文件
touch ~/.bash_alias
#> 将快捷启动conda的conact写入刚创建的~/.bash_alias
echo 'alias conact="source ~/miniconda3/bin/activate"' >> ~/.bash_alias
#> 将source ~/.bash_alias写到用户配置文件下，每次启动shell bash_alias中的快捷命令自动生效
echo 'source ~/.bash_alias' >> ~/.bashrc\
#> 重新加载下配置文件
source ~/.bash_alias
#> 利用conact启动conda
conact
#> 更换镜像源
vim ~/.condarc
#> 全程英文输入模式！！！
#> 键盘输入i进入编辑模式，左下角能看到状态变为insert
#> 将下列代码粘贴进去
channels:
  - defaults
show_channel_urls: true
default_channels:
  - https://mirrors.bfsu.edu.cn/anaconda/pkgs/main
  - https://mirrors.bfsu.edu.cn/anaconda/pkgs/r
  - https://mirrors.bfsu.edu.cn/anaconda/pkgs/msys2
custom_channels:
  conda-forge: https://mirrors.bfsu.edu.cn/anaconda/cloud
  msys2: https://mirrors.bfsu.edu.cn/anaconda/cloud
  bioconda: https://mirrors.bfsu.edu.cn/anaconda/cloud
  menpo: https://mirrors.bfsu.edu.cn/anaconda/cloud
  pytorch: https://mirrors.bfsu.edu.cn/anaconda/cloud
  pytorch-lts: https://mirrors.bfsu.edu.cn/anaconda/cloud
  simpleitk: https://mirrors.bfsu.edu.cn/anaconda/cloud
#> 按ESC退出编辑模式,然后按：号进入命令执行模式，做下角光标闪烁，输入wq回车
conda clean -i #清理index cache

新建一个emmax分析环境

做分析需要下面的软件：

软件	作用
bcftools	操作vcf文件，包括修改染色体prefix，按区间提取，构建索引，过滤snp等
vcftools	和bcftools类似
plink	将vcf转换为plink格式，然后进行vcf文件操作，emmax要求输入文件格式为plink
tassel	基因型文件转换格式，pca，生成phylo文件
fasttree	目标区间快速建立NJ树
admixture	群体结构分析
beagle	缺失基因型填补
tabix	vcf建立索引
gatk	snp calling
samtools	操作sam，bam文件对可以结构变异切割bam文件通过IGV验证
vep	SNP文件注释

代码如下：

#| code-fold: true
#| code-summary: "Show the code"
conda create -c bioconda -n gwas bcftools vcftools plink tassel admixture beagle fasttree tabix samtools gatk -y
#> 安装成功后可以看到以下提醒
# To activate this environment, use
#
#     $ conda activate gwas
#
# To deactivate an active environment, use
#
#     $ conda deactivate
#> 单独给vep创建一个环境
conda create -c bioconda -n ensembl-vep vep -y

安装emmax

emmax有可以直接在unbuntu20.04编译好的软件直接下载，解压使用，不过为了方便，我们还是需要些两个alias方便使用emmax和emmax-kin，或者直接把他们加入环境变量。

代码：

#| code-fold: true
#| code-summary: "Show the code"
wget http://csg.sph.umich.edu//kang/emmax/download/emmax-beta-07Mar2010.tar.gz
# --2023-06-14 11:05:10--  http://csg.sph.umich.edu//kang/emmax/download/emmax-beta-07Mar2010.tar.gz
# Resolving csg.sph.umich.edu (csg.sph.umich.edu)... 141.211.55.24
# Connecting to csg.sph.umich.edu (csg.sph.umich.edu)|141.211.55.24|:80... connected.
# HTTP request sent, awaiting response... 200 OK
# Length: 519928 (508K) [application/x-gzip]
# Saving to: ‘emmax-beta-07Mar2010.tar.gz’

# emmax-beta-07Mar2010.tar.gz          100%[=====================================================================>] 507.74K   416KB/s    in 1.2s

# 2023-06-14 11:05:13 (416 KB/s) - ‘emmax-beta-07Mar2010.tar.gz’ saved [519928/519928]
tar -xvf emmax-beta-07Mar2010.tar.gz
#> 将emmax文件夹放到家目录的software路径下
mkdir ~/software
mv emmax-beta-07Mar2010 ~/software/
#> 创建alias
echo 'alias emmax="~/software/emmax-beta-07Mar2010/emmax"' >> ~/.bash_alias
echo 'alias emmax-kin="~/software/emmax-beta-07Mar2010/emmax-kin"' >> ~/.bash_alias
source ~/.bash_alias
#> 测试
emmax
#> 出现下面提示说明安装完毕
# Usage: emmax [options]
# Required parameters
#   -t [tpedf_prefix] : prefix for tped/tfam files
#   -o [out_prefix]  : output file name prefix
# Likely essential parameters
#   -p [phenof] : 3-column phenotype file with FAMID, INDID at the first two colmns, in the same order of .tfam file. Not required only with -K option  -k [kinf] : n * n matrix containing kinship values in the individual order consistent to [tpedf].tfam file. [tpedf].kinf will be used if not specified
#   -c [covf] : multi-column covariate file with FAMID, INDID at the first two colmns, in the same order of .tfam fileOptional parameters
# Optional parameters
#   -i [in_prefix] : input file name prefix including eigenvectors
#   -d [# digits]  : precision of the output values (default : 5)
#   -s [start index of SNP] : start index of SNP (default : 0)
#   -e [end index of SNP] : end index of SNP (default : #snps)
#   -w : flag for writing eigenvalues/eigenvector files
#   -D [delimiters] : delimter string in quotation marks
#   -P [# heaer cols in tped] : # of column headers in tped file
#   -F [# heaer cols in tfam] : # of column headers in tfam file
# Segmentation fault (core dumped)

下载脚本文件

可以通过下载压缩包和通过git等方式将仓库代码下载到服务器，然后给脚本添加alias..

步骤	脚本名称	alias	作用
协变量文件格式化	MakeEmmaxCov.R	emmax-cov	修改协变量文件格式
表型数据文件格式化	split_phenotype.R	emmax-pheno-split	表型分割
开始关联分析及下游文件生成	run_emmax_pipeline.sh	run_emmax_pipeline	数据准备好后执行该脚本会将下游分析一键完成
开始关联分析及下游文件生成	EMMAX_DownStream.R	emmax-down	已经被集成到run_emmax-pipeline
区域单倍型划分	emmax-haplotype.sh	emmax-haplotype	根据候选区间进行单倍型划分
区域单倍型热图	emmax-haploheatmap.R	emmax-haploheatmap	倍集成到emmax-haplotype

代码：

#| code-fold: true
#| code-summary: "Show the code"
cd # 返回家目录
mkdir src && cd src # 创建一个src路径 然后进入该文件夹
git clone https://github.com/ShawnWx2019/emmax-pipeline.git # 通过git将emmax-pipeline下载到本地
echo "alias emmax-pca='Rscript ~/src/emmax-pipeline/code/PCA_var_proportion.R'" >> ~/.bash_alias
echo "alias emmax-cov='Rscript ~/src/emmax-pipeline/code/MakeEmmaxCov.R'" >> ~/.bash_alias
echo "alias emmax-down='Rscript ~/src/emmax-pipeline/code/EMMAX_DownStream.R'" >> ~/.bash_alias
echo "alias emmax-pheno-split='Rscript ~/src/emmax-pipeline/code/split_phenotype.R'" >> ~/.bash_alias
echo "alias run_emmax_pipeline='bash ~/src/emmax-pipeline/code/run_emmax_pipeline.sh'" >> ~/.bash_alias
echo "alias emmax-haplotype='bash ~/src/emmax-pipeline/code/emmax-haplotype.sh'" >> ~/.bash_alias
echo "alias emmax-haploheatmap='Rscript ~/src/emmax-pipeline/code/emmax-haploheatmap.R'" >> ~/.bash_alias
source ~/.bash_alias

至此，上游分析环境搭建工作完毕。

Emmax上游分析

基因型文件过滤

首先最原始的基因型文件应该是经过GATK call snp后进行质控后的.vcf文件, 可以直接使用vcftools进行进一步过滤，对于过滤参数的标准可以参考文献中对应物种的过滤方式。例如在@duResequencing243Diploid2018 文章中，陆地棉的过滤条件为：

only high-quality SNPs (coverage depth ≥3 and ≤50, mapping quality ≥20, missing ratio of samples within population ≤10% (1,980,539 SNPs) or 20% (3,665,030 SNPs), and minor allele frequency (MAF) ≥0.05) were retained for subsequent analyses. SNPs with a miss ratio ≤10% were used in PCA/phylogenetictree/structure analyses, whereas SNPs with a miss ratio ≤20% were used in the rest of the analyses.

第一步: 首先emmax要求染色体为纯数字，如果你发现GATK生成的vcf文件中染色体为"chr01"等并非纯数字时需要通过bcftools修改一下染色体。

#| code-fold: true
#| code-summary: "Show the code"
##> GATK 一般不会给snp加ID，这里用bcftools对vcf文件中snp命名，方式为chr:pos
bcftools annotate --set-id '%CHROM:%POS' raw.vcf.gz
#> 手动写一个配置文件，左边为原始的
vim chr_cn.txt
i ## 进入编辑模式
chr1    1
chr2    2
chr3    3
chr4    4
chr5    5

ESC #退出编辑模式
wq # 保存
## 使用bcftools 修改染色体名称
bcftools annotate --rename-chrs chr_cn.txt raw.vcf.gz | bgzip -c >raw_rename.vcf.gz
## 给新生成的vcf文件建立索引
bcftools index -t raw_rename.vcf.gz

第二步: 用plink对vcf文件过滤，生成的Gh_383.maf0.05.int0.8.vcf

#| code-fold: true
#| code-summary: "Show the code"
plink --vcf raw_rename.vcf.gz --maf 0.05 --geno 0.2 --recode vcf-fid --out Gh_383.maf0.05.int0.8
#> 生成emmax分析所需要的plink文件
plink --vcf Gh_383.maf0.05.int0.8.vcf --recode 12 transpose --output-missing-genotype 0 --out Gh_383.maf0.05.int0.9

第三步: 根据LD对SNP筛选，构建运行群体结构分析的基因型文件

#| code-fold: true
#| code-summary: "Show the code"
#> 过滤LD
plink --vcf Gh_383.maf0.05.int0.8.vcf --indep-pairwise 50 10 0.2 --out Gh_383.maf0.05.int0.8
#> 提取筛选SNP id
plink --vcf Gh_383.maf0.05.int0.8.vcf --make-bed --extract Gh_383.maf0.05.int0.8.prune.in --out Gh_383.maf0.05.int0.8.prune.in
#> 将筛选后的in结果转换为vcf
plink --bfile Gh_383.maf0.05.int0.8.prune.in --recode vcf-fid --out Gh_383.maf0.05.int0.8.prune.in
#> 转换为admixture需要的输入文件格式
plink -bfile Gh_383.maf0.05.int0.8.prune.in --recode 12 --out Gh_383.maf0.05.int0.8.prune.in

第四步: 利用admixture进行群体结构分析.

通过下面代码进行admixture群体结构分析，根据结果选择合适的K值，然后选择对应的群体结构协变量文件

#| code-fold: true
#| code-summary: "Show the code"
for i in {1..20}
do
admixture --cv Gh_383.maf0.05.int0.8.prune.in.ped ${i} -j48 >> log.txt
done
##> 运行结束后通过下面代码查看最佳分群
$ grep CV log.txt
##> emmax-cov进行格式转换

emmax-cov -n Gh_383.maf0.05.int0.8 \
          -q admixture.Q3
          -o Cov_Q3_emmax.cov
          -k 3

第五步: 利用emmax-kin生成亲缘关系矩阵

#| code-fold: true
#| code-summary: "Show the code"
emmax-kin Gh_383.maf0.05.int0.8 -v -d 10
#> 运行完毕Gh_383.maf0.05.int0.8.BN.kinf 文件

至此运行emmax的基因型文件，协变量文件已经具备了。

(vep) shawn @ bio-Super-Server: 01.raw_g $ tree
.
├── Cov_P3_emmax.cov # 协变量文件 PCA
├── Cov_Q3_emmax.cov # 协变量文件 群体结构
├── Gh_383.maf0.05.int0.8.bed # 基因型文件plink格式 二进制
├── Gh_383.maf0.05.int0.8.bim # 基因型文件plink格式 二进制
├── Gh_383.maf0.05.int0.8.fam # 基因型文件plink格式 二进制
├── Gh_383.maf0.05.int0.8.hBN.kinf #亲缘关系矩阵
├── Gh_383.maf0.05.int0.8.log # 基因型文件plink格式 二进制
├── Gh_383.maf0.05.int0.8.nosex # 基因型文件plink格式
├── Gh_383.maf0.05.int0.8.tfam # 基因型文件plink格式
├── Gh_383.maf0.05.int0.8.tped # 基因型文件plink格式
├── Gh_383.maf0.05.int0.8.vcf # 基因型文件 vcf格式

SNP注释文件准备

通过vep软件对snp文件进行注释

代码：

# 格式化
grep -v "#" CRI_Gh_v2.gff3 | sort -k1,1 -k4,4n -k5,5n -t$'\t' | bgzip -c > data.gff.gz
# tabix索引
tabix -p gff data.gff.gz
# 注释
vep -i  raw.vcf.gz  --gff data.gff.gz --fasta CRI_Gh_v2.fa --fork 4 --tab -o Gh_383_anno.xls Gh_383_anno.html

表型数据准备

首先按照下表准备好输入表型文件traits.txt

sample	trait1	trait2	trait3
s_001	212	123	NA
s_002	221	1223	2
s_004	NA	123	3

然后运行

emmax-pheon-split -t Gh_383.maf0.05.int0.8.tfam -p traits.txt

然后我们就得到了trait1.txt trait2.txt trait3.txt可以用于emmax表型分析的文件

至此，所有运行下游分析的文件都准备完毕

Emmax一步流程

首先我们看下帮助文档

-------------------------------

Emmax pipeline . From raw to Manhattan plot

-------------------------------

Usage:

run_emmax_pipline \
 -t [tpedf_prefix] \
 -o [out_prefix] \
 -p [phenof] \
 -k [kinf] \
 -a [anno] \
 -c [covf] (option) \
 -i [img_type] (option) \
 -s [point_size] (option) \
 -w [point_color] (option) \
 -r [rerun] (option) \


-------------------------------

Author Shawn Wang (shawnwang2016@126.com)
Date Web Apr 26, 2023
Version V.0.0.0.99 beta

-------------------------------

Required parameters:

[-t]:tpedf_prefix,   prefix for tped/tfam files
[-o]:out_prefix,  output file name prefix
[-p]:phenof,  3-column phenotype file with FAMID, INDID at the first two colmns, in the same order of .tfam file.
[-k]:kinf,  n * n matrix containing kinship values in the individual order consistent to [tpedf].tfam file. [tpedf].kinf will be used if not specified.
[-a]:anno,  SNP annotation file by vep.

Optional parameters:

[-c]:covf,  multi-column covariate file with FAMID, INDID at the first two colmns, in the same order of .tfam fileOptional parameters.
[-i]:img_type,  output img file type, jpg or pdf.
[-s]:point_size,  point size of manhattan plot, default: cut_off 0_4_6 ==> 0.3_0.5_0.5.
[-w]:point_color,  point color of manhattan plot, default: cut_off 4_6 ==> green_red.
[-r]:rerun,  If you need to rerun, please use this parameter. 1, 2 or 12

-------------------------------

必选参数 -t plink 生成的转置的emmax文件 -o 输出文件的前缀 -p 表型文件 -k 亲缘关系矩阵 -a 注释文件
可选参数 -c 协变量文件 -i 输出图片文件格式 -s 曼哈顿图点的大小 -w 显著位点点的颜色（-logp 4，6） -r 断点重新运行

run_emmax_pipeline \
  -t Gh_383.maf0.05.int0.8 \
  -o Morin_raw \
  -p Morin.txt \
  -k Gh_383.maf0.05.int0.8.hBN.kinf \
  -a Gh_383_anno.xls \
  -c Cov_Q3_emmax.cov

大概等待10分钟左右就可以出结果，基因型文件较小的话更快。

运行完毕得到的文件如下：

── LDFilebyChr => #我们根据显染色体上显著位点的个数提取了前5个chr的基因型文件数据用来做LDblock
│   ├── A01_genotype.txt
│   ├── A11_genotype.txt
│   └── D07_genotype.txt
├── Morin_Flower_miss.log
├── Morin_Flower_miss.ps => # emmax 输出的包含位点p值的结果文件
├── Morin_Flower_miss.reml
├── Morin_Flower_miss_result_filter.gwas => #格式化后的结果文件，用与使用IGV看结果
├── Morin_Flower_miss_result_filter.xlsx => #Excel文件，包含显著位点的snp信息，包括了p值和snp注释
├── Morin_Flower_miss_result_genotype_all.txt
├── QQplot.Morin_Flower_miss.jpg => #QQ 图
└── Rectangular-Manhattan.Morin_Flower_miss.jpg => #曼哈顿图

Emmax 下游分析

确定LDblock画图范围

首先可以通过Morin_Flower_miss_result_filter.xlsx 文件中连续高信号出现染色体位置划分范围。其次可以通过用IGV看Morin_Flower_miss_result_filter.gwas 划定范围

进行LDblock可视化

通过LDblock进行可视化

经过测试，报错的原因主要是染色体名字不符，许多软件在做gwas分析的时候都会改变chr的前缀，导致基因型文件中染色体标志和gff文件中的不符合，此时需要用bcftools重新修改基因型文件。

LDBlockShow -InGFF  ~/GeekCloud/01.Database/Cotton/AD1/Gh.gff -Cutline 6 -OutPdf -SeleVar 2 -TopSite -InGenotype A11_Gh_383_miss.geno --OutPut Morin_A11_BGLU --InGWAS A11_genotype.txt -Region A11:91108756:91221656 &

License

ShawnWx2019/emmax-pipeline

ErrorLooks like something went wrong!

About

Topics

Resources

License

Uh oh!

Stars

Watchers

Forks

Releases

Packages 0

Languages

Error
Looks like something went wrong!

Packages