最近开始对重测序数据进行分析，针对构建系统发育树做记录：
进入环境：
source activate bininfo
文件准备：SNP的.vcf文件
软件准备：bgzip（压缩为gz格式）、bcftools（合并文件）、plink（plink2）（数据过滤）、phylip（构树）、tassel（转换格式）。
①压缩文件：公司返回的vcf文件是单独的，即一个物种有一个snp.vcf，我们需要合成一个总vcf文件进行构树。另外有时公司返回的文件已经压缩snp.vcf.gz，需要用bzip解压再压缩。可在centos下文件管理解压。
bgzip -c -f xxx1.vcf > xxx1.vcf.gz
bgzip -c -f xxx2.vcf > xxx2.vcf.gz
②合并文件
bcftools index xxx1.vcf.gz
bcftools index xxx2.vcf.gz
bcftools merge xxx1.vcf.gz xxx2.vcf.gz -Oz -o merge.vcf.gz
输出了合并文件merge.vcf.gz，之后解压得到最终的合并文件merge.vcf
③plink（原文//www.greatytc.com/p/44993691145f）
过滤SNP
plink --vcf merge.vcf --geno 0.1 --maf 0.01 --out all.missing_maf --recode vcf-iid --allow-extra-chr --set-missing-var-ids @:# --keep-allele-order
如果提示memory报错，则在后面加--memory 报错数字。
plink --vcf merge.vcf --geno 0.1 --maf 0.01 --out all.missing_maf --recode vcf-iid --allow-extra-chr --set-missing-var-ids @:# --keep-allele-order --memory 217364096
如果pink持续报错，度娘不出来解决办法，使用plink2试试。
输出文件all.missing_maf
生成保留位点文件
plink --vcf all.missing_maf.vcf --indep-pairwise 50 10 0.2 --out tmp.ld --allow-extra-chr --set-missing-var-ids @:#
过滤LD
plink --vcf all.missing_maf.vcf --make-bed --extract tmp.ld.prune.in --out all.LDfilter --recode vcf-iid --keep-allele-order --allow-extra-chr --set-missing-var-ids @:#
输出文件all.LDfilter.vcf
④构树
将vcf转为phylip格式文件
run_pipeline.pl -Xms1G -Xmx5G -importGuess all.LDfilter.vcf -ExportPlugin -saveAs sequences.phy -format Phylip_Inter
生成dnadist需要的配置文件
echo -e "sequences.phy\nY" > dnadist.cfg
运行dnadist生成距离矩阵文件
dnadist < dnadist.cfg >dnadist.log
生成neighbor程序需要的配置文件
mv outfile infile.dist
echo -e "infile.dist\nY" > neighbor.cfg
构建nj树
neighbor < neighbor.cfg >nj.log
整理结果格式
less infile.dist | tr '\n' '|'| sed 's/| / /g' | tr '|' '\n' >infile.dist.table
less outtree | tr '\n' ' '|sed 's/ //g' > outtree.nwk
输出文件outtree.nwk，用MEGA进行优化。
2022.10.30 by ZXG
