最近開(kāi)始對(duì)重測(cè)序數(shù)據(jù)進(jìn)行分析，針對(duì)構(gòu)建系統(tǒng)發(fā)育樹(shù)做記錄：
進(jìn)入環(huán)境：
source activate bininfo
文件準(zhǔn)備：SNP的.vcf文件
軟件準(zhǔn)備：bgzip（壓縮為gz格式）、bcftools（合并文件）、plink（plink2）（數(shù)據(jù)過(guò)濾）、phylip（構(gòu)樹(shù)）、tassel（轉(zhuǎn)換格式）。
①壓縮文件：公司返回的vcf文件是單獨(dú)的，即一個(gè)物種有一個(gè)snp.vcf，我們需要合成一個(gè)總vcf文件進(jìn)行構(gòu)樹(shù)。另外有時(shí)公司返回的文件已經(jīng)壓縮snp.vcf.gz，需要用bzip解壓再壓縮?？稍赾entos下文件管理解壓。
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（原文http://www.itdecent.cn/p/44993691145f）
過(guò)濾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報(bào)錯(cuò)，則在后面加--memory 報(bào)錯(cuò)數(shù)字。
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持續(xù)報(bào)錯(cuò)，度娘不出來(lái)解決辦法，使用plink2試試。
輸出文件all.missing_maf
生成保留位點(diǎn)文件
plink --vcf all.missing_maf.vcf --indep-pairwise 50 10 0.2 --out tmp.ld --allow-extra-chr --set-missing-var-ids @:#
過(guò)濾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
④構(gòu)樹(shù)
將vcf轉(zhuǎn)為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
運(yùn)行dnadist生成距離矩陣文件
dnadist < dnadist.cfg >dnadist.log
生成neighbor程序需要的配置文件
mv outfile infile.dist
echo -e "infile.dist\nY" > neighbor.cfg
構(gòu)建nj樹(shù)
neighbor < neighbor.cfg >nj.log
整理結(jié)果格式
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進(jìn)行優(yōu)化。
2022.10.30 by ZXG
