(本文软件和操作均基于Linux,我的SNP数量1400多万,个体数346个)
在群体遗传学中,一个物种内,多个材料的基因型文件可以构建系统发育树,来推断群体材料的分化关系。同时,可以使用群体结构图(Structure或者fastStructure),展示群体的分层现象。Structure软件运行较慢,而fastStructure则速度占优,因而考虑使用这个软件。同时R包LEA也可以计算和绘制群体结构图。
群体材料系统发育树的构建方法和软件非常多元,可以使用SNPhylo、MEGA X、fastTree、IQtree和TASSEL等,本文选择SNPhylo这个软件,它可以直接输入VCF文件格式或HapMap格式,非常方便。
1.1 fastStructure安装步骤(python2.7)fastStructure是python软件编写的软件,依赖较多,官网编译下载复杂,使用conda或者docker更方便,二者优先考虑conda,docker虽然轻量,它的image通常较大。
## 创建conda环境,下载fastStructure(Python2.7)conda create -n faststructure_py2 python=2.7mamba install -c bioconda faststructure
1.2 fastStructure使用步骤fastStructure使用bed格式基因型作为输入(plink生成):
## 将下面plink格式(map/ped)转为bed格式## map和ped文件可以使用vcftools、plink或者TASSEL均能生成-rwxrwxrwx 1 root root 19525144502 Feb 4 00:32 chrall.plink.ped-rwxrwxrwx 1 root root 368083368 Feb 4 00:04 chrall.plink.map## 生成bed格式基因型plink --file chrall.plink --make-bed --out chrallfas## 基因型结果如下-rwxrwxrwx 1 root root 1112 Feb 4 13:25 chrallfas.log-rwxrwxrwx 1 root root 424514400 Feb 4 13:25 chrallfas.bim-rwxrwxrwx 1 root root 7776 Feb 4 13:24 chrallfas.fam-rwxrwxrwx 1 root root 1227374949 Feb 4 13:24 chrallfas.bed-rwxrwxrwx 1 root root 4662 Feb 4 13:24 chrallfas.nosex
然后运行下面的主体计算命令:
for i in {2..10}donohup python /home/cfc424/bin/miniconda3/envs/faststructure_py2/bin/structure.py -K ${i} --input=chrallfas --output=genotypes_output &done## 筛选最合适的K群体个数和绘制图片python chooseK.py --input=genotypes_output ## 假定K=3python distruct.py -K 3 --input=genotypes_output --output=genotypes_output.svg
实际绘图时,可以根据材料的先验知识,选择合适的K值绘图。
我在绘图步骤,遇到报错:
## 运行python distruct.py -K 3 --input=genotypes_output --output=genotypes_output.svg ## 报错Tkinter.py, line 1819, in __init__ self.tk = _tkinter.create(screenName, baseName, className, interactive, wantobjects, useTk, sync, use)_tkinter.TclError: no display name and no $DISPLAY environment variable
查询发现这是distruct.py有问题,文件header处增加下面两行即可:
import matplotlib as mpl ## add line onempl.use('svg') ## add line two
修改后distruct.py的header如下所示(注意导入顺序):
import numpy as npimport matplotlib as mpl ## add line onempl.use('svg') ## add line two import matplotlib.pyplot as plotimport colorsysimport getoptimport sys, pdb
保存后,再运行:
## K=3时,structure.py结果如下-rwxrwxrwx 1 root root 819 Feb 6 02:04 genotypes_output.3.log-rwxrwxrwx 1 root root 391M Feb 6 02:05 genotypes_output.3.meanP-rwxrwxrwx 1 root root 9.8K Feb 6 02:04 genotypes_output.3.meanQ## 通过上面的结果绘图python /home/cfc424/bin/miniconda3/envs/faststructure_py2/bin/distruct.py -K 3 --input=genotypes_output --output=genotypes_output3.svg ## 输出结果genotypes_output3.svg## 使用Python3软件包 cairosvg将svg格式转换为pdf或者png等## 注意该软件的Python3和fastStructure的Python2版本不同## 需要在Python3环境中操作(我的conda base为Python3)## 此外,使用Adobe Illustrator导入svg文件调节后导出图片mamba install cairosvgcairosvg genotypes_output3.svg -o genotypes_output3.pdffor i in {2..10};do cairosvg genotypes_output${i}.svg -o genotypes_output${i}.png ;done
*** Note:
fastStructure生成的图片后面调节颜色是个问题,特别是多个K值对应时,这个我还没弄。
SNPhylo是计算群体进化树的软件,依赖软件多,使用conda安装即可:
conda create -n snphylo_r4.0_py2 r=4.1 python=2.7conda activate snphylo_r4.0_py2
在github release中下载SNPhylo的安装包,地址如下:
https://github.com/thlee/SNPhylo/releases/download/20180901/SNPhylo.20180901.zip
解压文件,进入配置命令,检查软件依赖软件:
bash ./setup.shR ## R4.1python ## python2.7muscle ## 比对软件dnaml 缺少R包:"getopt", "gdsfmt","SNPRelate","ape","phangorn"
提示环境中R包不存在,下载这些R包:
if (!requireNamespace("BiocManager", quietly = TRUE));install.packages("BiocManager") BiocManager::install(c("getopt", "gdsfmt","SNPRelate","ape","phangorn"))
重新运行配置脚本,软件安装成功:
bash ./setup.sh## 最终提示SNPHYLO is successfully installed!!!
2.1 SNPhylo使用 SNPhylo软件接受VCF文件、HapMap文件和Simple SNP文件作为输入。
使用VCF文件的命令例子如下:
~/bin/SNPhylo-master/snphylo.sh -v ./Chrall.filter.maf0.05_maxmissing0.8.vcf -p 5 -c 5
SNPhylo软件输出结果如下,其中包括snphylo.output.ml.tree为Newwick格式的tree文件,可以将该文件修改结尾为.nwk,导入软件MEGA X或者ITOL(网络平台)或者Evolview(网络平台)以及R包ggtree进行调节展示。
如果要求不高,该软件已经提供PNG图片snphylo.output.ml.png。
-rwxrwxrwx 1 root root 88116 Feb 4 09:09 snphylo.output.ml.png-rwxrwxrwx 1 root root 8664 Feb 4 09:09 snphylo.output.ml.tree-rwxrwxrwx 1 root root 94320 Feb 4 09:09 snphylo.output.ml.txt-rwxrwxrwx 1 root root 3682994 Feb 4 02:46 snphylo.output.phylip.txt-rwxrwxrwx 1 root root 3347459 Feb 4 02:46 snphylo.output.fasta-rwxrwxrwx 1 root root 121234 Feb 4 02:44 snphylo.output.id.txt-rwxrwxrwx 1 root root 344191155 Feb 4 02:43 snphylo.output.gds-rwxrwxrwx 1 root root 5212597003 Feb 4 00:50 snphylo.output.filtered.vcf
以上,对于我的数据集(SNP 1400多万,个体数346个),fastStructure(K=2 …K=10)共运行67个小时(3天不到,K越大,速度越慢),SNPhylo运行7.5小时(服务器500G运行内存,逻辑cpu 8个)。
两个的软件具体参数可以查阅官网文档。
参考
https://rajanil.github.io/fastStructure/
https://github.com/thlee/SNPhylo
https://cloud.tencent.com/developer/news/311649
https://www.jianshu.com/p/64db4635e511
https://github.com/rajanil/fastStructure/issues/36