生信：種群歷史有效群體大小推斷 SMC++(二)

2023-07-21 00:00 作者:國靚 0人讀過 | 我要投稿

smc++ plot?繪圖效果

cluster文件，兩列第一列樣本名稱第二列g(shù)roup名稱

代碼塊1

cat?../00.raw/cluster?|?cut?-f?2?|?sort?-u?|?while?read?-r?a;?do ????awk?-v?var="$a"?'$2?==?var?{?sp?=?sp?$1?","?}?END?{?print?var?":"?sp?}'?../00.raw/cluster?|?sed?'s/,$//'?|?awk?-v?var="$a"?'$0?!~?var"$"'?>?"${a}.list" done

目標(biāo)：生成樣品信息文件，格式為?pop:sp1,sp2?的形式。
結(jié)果：對于輸入文件?../00.raw/cluster，根據(jù)第2列的唯一值進(jìn)行循環(huán)。使用?awk?命令提取符合條件的樣品信息，并使用?sed?命令去掉末尾的逗號，并使用第2列的值作為文件名保存。

代碼塊2

cat?../00.raw/cluster?|?cut?-f?2?|?sort?-u?|?while?read?-r?a;?do ????awk?-v?var="$a"?'$2?==?var?{print?$1}'?../00.raw/cluster?|?head?-n?4?>?"${a}.sample.list" done

目標(biāo)：生成Composite likelihood樣品列表文件，一般選取2-10個樣品。
結(jié)果：對于輸入文件?../00.raw/cluster，根據(jù)第2列的唯一值進(jìn)行循環(huán)。使用?awk?命令提取符合條件的樣品名，并使用?head?命令限制數(shù)量為4個，并保存到以第2列的值為文件名的文件中。

代碼塊3

bgzip?-c?/data/Work/22.selective_sweeps/00.raw/all_snps_imputation.eff.reh.chr.vcf?>?../00.raw/all_snps_imputation.eff.reh.chr.vcf.gz tabix?../00.raw/all_snps_imputation.eff.reh.chr.vcf.gz

目標(biāo)：將 VCF 文件進(jìn)行壓縮，并生成索引文件。
結(jié)果：將輸入的 VCF 文件?/data/zhougl/Work/22.selective_sweeps/00.raw/all_snps_imputation.eff.reh.chr.vcf?使用?bgzip?命令進(jìn)行壓縮，并將壓縮文件保存為?../00.raw/all_snps_imputation.eff.reh.chr.vcf.gz。然后使用?tabix?命令生成該壓縮文件的索引。

代碼塊4

cat?../00.raw/cluster?|?cut?-f?2?|?sort?-u?|?while?read?-r?a;?do ????mkdir?"${a}_vcf2smc" done

目標(biāo)：創(chuàng)建文件夾以準(zhǔn)備進(jìn)行 VCF 轉(zhuǎn) SMC 輸入格式的操作。
結(jié)果：對于輸入文件?../00.raw/cluster，根據(jù)第2列的唯一值進(jìn)行循環(huán)。使用?mkdir?命令創(chuàng)建名為?${a}_vcf2smc?的文件夾。

代碼塊5

awk?'BEGIN?{for?(i=1;?i<=12;?i++)?{printf?"%d\t%s%02d\n",?i,?"chr",?i}}'?|?cut?-f?2?>?chr.list

目標(biāo)：生成染色體列表文件。
結(jié)果：使用?awk?命令生成以?chr?開頭的數(shù)字編號的染色體列表，并將結(jié)果保存到?chr.list?文件中。

代碼塊6

cp?tmp.sh?G2_vcf2smc.sh cat?chr.list?|?while?read?chr;?do ????cat?G2.sample.list?|?while?read?sample;?do ????????echo?"smc++?vcf2smc?-m?$mask?-d?$sample?$sample?$vcf?G2_vcf2smc/$sample.$chr.smc.gz?$chr?`cat?G2.list`?"?>>?G2_vcf2smc.sh ????done done sbatch?G2_vcf2smc.sh

目標(biāo)：執(zhí)行?smc++ vcf2smc?命令將 VCF 文件轉(zhuǎn)換為 SMC 輸入格式。
結(jié)果：首先復(fù)制?tmp.sh?文件為?G2_vcf2smc.sh。然后，對于每個染色體和樣品，使用?smc++ vcf2smc?命令將 VCF 文件轉(zhuǎn)換為 SMC 輸入格式，并保存到對應(yīng)的輸出文件中。最后，使用?sbatch?命令提交?G2_vcf2smc.sh?腳本進(jìn)行執(zhí)行。

代碼塊7

mkdir?G2_analysis smc++?estimate?--cores?10?--knots?10?--spline?cubic?-o?G2_analysis?6.5e-9?G2_vcf2smc/*.smc.gz

目標(biāo)：進(jìn)行種群歷史的有效群體大小估計。
結(jié)果：創(chuàng)建名為?G2_analysis?的文件夾。然后，使用?smc++ estimate?命令進(jìn)行種群歷史的有效群體大小估計，指定核心數(shù)為 10，節(jié)點(diǎn)數(shù)為 10，樣條函數(shù)為 cubic，并將結(jié)果輸出到?G2_analysis?文件夾中，使用突變速率參數(shù)為?6.5e-9，輸入文件為?G2_vcf2smc/*.smc.gz。

代碼塊8

smc++?plot?-g?1?--linear?G2_analysis.plot.pdf?G2_analysis/model.final.json

目標(biāo)：繪制種群歷史估計結(jié)果的圖表。
結(jié)果：使用?smc++ plot?命令，指定世代時間比例為 1，繪制線性圖表，并將結(jié)果輸出到?G2_analysis.plot.pdf?文件中，輸入文件為?G2_analysis/model.final.json。

國靚的文案

主要分享：1：基于R語言、python和shell的數(shù)據(jù)分析，數(shù)據(jù)可視化及生物信息分析；2：零基礎(chǔ)學(xué)英語從How are you開始學(xué)；3：運(yùn)動健康。（能量是守恒的，喜歡是互相的，關(guān)注我，世界上就多了一個愛你的人！）