🏡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
完成环境搭建及软件安装。下面将讲解详细步骤。
首先到清华镜像站下载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
做分析需要下面的软件:
软件 | 作用 |
---|---|
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有可以直接在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
至此,上游分析环境搭建工作完毕。
首先最原始的基因型文件应该是经过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格式
通过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 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 => #曼哈顿图
首先可以通过Morin_Flower_miss_result_filter.xlsx 文件中连续高信号出现染色体位置划分范围。 其次可以通过用IGV看Morin_Flower_miss_result_filter.gwas 划定范围
通过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 &