1、轉(zhuǎn)換文件格式

從 VCF 文件轉(zhuǎn)換為 relate 軟件可以讀取的 hap 格式：

PATH_TO_RELATE/bin/RelateFileFormats \

? ? ? ? ? ? ? ? --mode ConvertFromVcf \

? ? ? ? ? ? ? ? --haps example.haps \

? ? ? ? ? ? ? ? --sample example.sample \

? ? ? ? ? ? ? ? -i example

--haps? 指定輸出的.haps文件路徑。

--sample? 指定輸出的.sample文件路徑。

-i, --input? 輸入文件的基名（不帶擴(kuò)展名），并優(yōu)先選擇非壓縮文件（如無壓縮文件則使用壓縮文件）。

--chr? 可選參數(shù)：指定染色體索引（整數(shù)類型），將作為HAPS文件的第一列。默認(rèn)值：0。

輸出文件為 example.haps 和 example.sample。

# 這里的 VCF 文件需要進(jìn)行 phase。

java -Xmx163840m -jar? ?beagle.18May20.d20.jar? ?\

? ? ? ?gt=chr01.vcf? ?\

? ? ? ?out=chr01.phased? ?\

? ? ? ?window=8? #默認(rèn)40，縮小可以提高計(jì)算速度，增加可以提高估計(jì)的準(zhǔn)確性 \

? ? ? ?ne=100000? #默認(rèn)10e6 \

? ? ? ?nthreads=8

2、估計(jì)譜系進(jìn)化歷史

首先需要進(jìn)行一次簡(jiǎn)單的計(jì)算，后續(xù)所有計(jì)算都需要在此基礎(chǔ)上運(yùn)行，運(yùn)算時(shí)間較短。

PATH_TO_RELATE/bin/Relate \

? ? ? --mode All \

? ? ? -m 1.25e-8 \

? ? ? -N 30000 \

? ? ? --haps data/example.haps \

? ? ? --sample data/example.sample \

? ? ? --map data/genetic_map.txt \

? ? ? --annot data/example.annot \

? ? ? --seed 1 \

? ? ? -o example

輸出文件分別為 example.anc 和 example.mut。

--mode? 指定Relate的運(yùn)行模式。模式"All"將執(zhí)行算法的所有階段，其他模式可用于單獨(dú)執(zhí)行特定階段（參見文檔）。此功能在并行化運(yùn)行Relate時(shí)非常有用。

-m, --mutation_rate? ? 每代每個(gè)堿基的突變率。

-N, --effective_size? 單倍型的有效群體大小。（指單倍型而非個(gè)體！若需計(jì)算個(gè)體的有效群體大小，需將單倍型的有效群體大小乘以2）

example.haps 和 example.sample 為上一步驟生成的文件；

genetic_map.txt 則需要單獨(dú)生成，具體格式如下：

將三列中的第i個(gè)元素分別記為p[i]、r[i]、rdist[i]，通過以下公式計(jì)算：r[i]=( (rdist[i+1]?rdist[i])/(p[i+1]?p[i]) )×1e6。

example.annot 是可選項(xiàng)，如果想放，按如下方式生成：

PATH_TO_RELATE/bin/RelateFileFormats \

? ? ? ? ? ? ? ? --mode GenerateSNPAnnotations \

? ? ? ? ? ? ? ? --haps example.haps \

? ? ? ? ? ? ? ? --sample example.sample \

? ? ? ? ? ? ? ? --poplabels example.poplabels \

? ? ? ? ? ? ? ? --ancestor ancestor.fa \

? ? ? ? ? ? ? ? -o example

example.poplabels 文件指定樣本所在群體信息，按如下格式整理：

ancestor.fa 是可選項(xiàng)，如果要放，需要是和數(shù)據(jù)集比對(duì)后的 fasta 格式。

3、估計(jì)有效群體大小

PATH_TO_RELATE/scripts/EstimatePopulationSize/EstimatePopulationSize.sh \

? ? ? ?-i example \

? ? ? ?-m 1.25e-8 \

? ? ? ?--poplabels data/example.poplabels \

? ? ? ?--seed 1 \

? ? ? -o example_popsize

較上一步驟，該步驟需要的時(shí)間較久。

這里 -i 的輸入文件為上一步驟生成的 .anc 和 .mut 文件，無需后綴。

此外，有一些可選項(xiàng)可供選擇：

--noanc? ? 可選參數(shù)：指定是否在最后重新推斷分支長(zhǎng)度。0：重新推斷，1：不重新推斷。默認(rèn)值：0。

--pop_of_interest? ? 可選參數(shù)：指定需要估計(jì)種群大小的群體。使用.poplabels文件的第二列作為群體標(biāo)簽?？梢灾付ǘ鄠€(gè)群體，用逗號(hào)分隔（例如pop1,pop2）。使用此選項(xiàng)時(shí)，建議（但非必需）在運(yùn)行Relate時(shí)添加--annot。如果未使用--annot運(yùn)行Relate，也可后續(xù)補(bǔ)充（詳見文檔）。

--first_chr& --last_chr? ? ?可選參數(shù)：允許多個(gè)輸入文件，假設(shè)文件名按染色體索引編號(hào)。示例：-i example --first_chr 1 --last_chr 2 會(huì)假設(shè)輸入文件為 example_chr1.anc, example_chr1.mut, example_chr2.anc, example_chr2.mut。

--chr? ? 可選參數(shù)：包含染色體名稱的文件（每行一個(gè)名稱）。此選項(xiàng)會(huì)覆蓋--first_chr和--last_chr的設(shè)置。示例：-i example --chr chr.txt 會(huì)假設(shè)輸入文件為example_chr*.anc,example_chr*.mut。

--threads? ? 可選參數(shù)：最大線程數(shù)。默認(rèn)值：1。

--bins? ? 可選參數(shù)：指定分箱區(qū)間。格式為lower,upper,stepsize，對(duì)應(yīng)函數(shù)c(0,10^seq(lower, upper, stepsize))。

--years_per_gen? ? 可選參數(shù)：每代的年數(shù)。默認(rèn)值：28.0。

--num_iter? ? 可選參數(shù)：算法迭代次數(shù)。默認(rèn)值：5。

--threshold? ? 可選參數(shù)：[0,1]范圍內(nèi)的浮點(diǎn)數(shù)，表示按突變映射數(shù)量排序后需丟棄的樹的比例。此參數(shù)有兩個(gè)作用：1) 移除低質(zhì)量/無信息量的樹；2) 加速推斷。值越大，丟棄的樹越多。默認(rèn)值：0.5。

--seed? ? 可選參數(shù)：用于估計(jì)分支長(zhǎng)度的MCMC算法中的隨機(jī)數(shù)生成器種子。

輸出文件：

example_popsize.pdf????展示估計(jì)的種群大小的圖表。

example_popsize.anc.gz, example_popsize.mut.gz, example_popsize.dist

如果指定--noanc 1，則不會(huì)輸出這些文件。包含更新分支長(zhǎng)度后的譜系信息。無論threshold參數(shù)值如何，均會(huì)保留所有樹的信息。若指定了--noanc 1，可通過獨(dú)立運(yùn)行 ReEstimateBranchLengths.sh 腳本重新估計(jì)分支長(zhǎng)度。

example_popsize.coal????包含將所有樣本視為單一群體時(shí)的共祖率（coalescence rates）和交叉共祖率（cross-coalescence rates）。

example_popsize.pairwise.coal 和 .bin? ? 共祖率文件及其對(duì)應(yīng)的二進(jìn)制文件，記錄樣本對(duì)之間的共祖率?？赏ㄟ^ RelateCoalescentRates --mode FinalizePopulationSize 提取任意 .poplabels文件的共祖率。

example_popsize_avg.rate????包含時(shí)間平均的突變率。

4、重新估計(jì)分支長(zhǎng)度

基于上述估計(jì)的有效群體大小重新估計(jì)分支長(zhǎng)度。

PATH_TO_RELATE/scripts/SampleBranchLengths/ReEstimateBranchLengths.sh \

? ? ? ? ? ? ? ? -i example \

? ? ? ? ? ? ? ? -o example_sub \

? ? ? ? ? ? ? ? -m 1.25e-8 \

? ? ? ? ? ? ? ? --coal popsize.coal \

? ? ? ? ? ? ? ? --first_bp 100000 \

? ? ? ? ? ? ? ? --last_bp 200000 \

? ? ? ? ? ? ? ? --seed 1

-i, --input? ? 輸入文件名（不包含擴(kuò)展名），應(yīng)為.anc和.mut文件的名稱。需要包含輸入haps/sample文件中的所有 SNP。如果未包含所有 SNP，請(qǐng)參見--dist選項(xiàng)。

-o, --output? ? 輸出文件的基名（不包含擴(kuò)展名）。

-m, --mu? ? 突變率，例如1.25e-8。

--coal? ? 估計(jì)的共祖率，將所有單倍型視為一個(gè)群體。

--threads? ? 可選參數(shù)：最大線程數(shù)。默認(rèn)值為 1。

--first_bp? ? 可選參數(shù)：起始?jí)A基對(duì)位置（bp）。

--last_bp? ? 可選參數(shù)：結(jié)束堿基對(duì)位置（bp）。

--dist? ? 可選參數(shù)：包含 SNP 堿基對(duì)位置及其相鄰 SNP 間距離的文件名（文件格式）。例如，可通過RelateExtract --mode AncMutForSubregion獲取。注意：如果.anc/.mut文件僅覆蓋輸入haps/sample文件的子區(qū)域，則必須提供此文件。

--seed? ? 可選參數(shù)：用于分支長(zhǎng)度估計(jì)的隨機(jī)種子。

輸出文件為 example_sub.anc/mut.

5、估計(jì)突變率

計(jì)算隨時(shí)間變化的突變率：

PATH_TO_RELATE/bin/RelateMutationRate \

? ? ? ? ? ? ? ? --mode Avg\

? ? ? ? ? ? ? ? -i example \

? ? ? ? ? ? ? ? -o example

-i, --input? ? 輸入文件名（不包含擴(kuò)展名），應(yīng)為.anc和.mut文件的名稱。

-o, --output? ? 輸出文件的基名（不包含擴(kuò)展名）。

--first_chr & --last_chr? ? 可選參數(shù)：允許多個(gè)輸入文件，假設(shè)文件名按染色體索引編號(hào)。示例：-i example --first_chr 1 --last_chr 2會(huì)假設(shè)輸入文件為example_chr1.anc,example_chr1.mut,example_chr2.anc,example_chr2.mut。

--chr? ? 可選參數(shù)：包含染色體名稱的文件（每行一個(gè)名稱）。此選項(xiàng)會(huì)覆蓋--first_chr和--last_chr的設(shè)置。示例：-i example --chr chr.txt會(huì)假設(shè)輸入文件為example_chr*.anc,example_chr*.mut。

--bins? ? 可選參數(shù)：指定分箱區(qū)間。格式為 lower, upper, stepsize，對(duì)應(yīng)函數(shù)c(0,10^seq(lower, upper, stepsize))。

--years_per_gen? ? 可選參數(shù)：每代的年數(shù)。默認(rèn)值：28.0。

--dist? ? 可選參數(shù)：指定.dist文件。該文件存儲(chǔ) SNP 之間的堿基對(duì)距離（單位：bp）。如果未指定，則使用.mut文件中提供的距離。僅在通過某些功能（如此函數(shù)）移除部分 SNP（及對(duì)應(yīng)的樹）時(shí)需要此選項(xiàng)。.dist文件可通過此功能生成。注意：如果與--first_chr或--last_chr同時(shí)指定，請(qǐng)?zhí)峁┎粠U(kuò)展名的文件名。

輸出文件為：example_avg.rate.

6、檢測(cè)正選擇

計(jì)算譜系選擇證據(jù)的p值。建議先估計(jì)種群大小歷史，并設(shè)置 --threshold 0。此操作會(huì)基于估計(jì)的種群大小歷史更新所有樹的分支長(zhǎng)度。

PATH_TO_RELATE/scripts/DetectSelection/DetectSelection.sh \

? ? ? ? ? ? ? ? -i example \

? ? ? ? ? ? ? ? -o example_selection \

? ? ? ? ? ? ? ? -m 1.25e-8 \

? ? ? ? ? ? ? ? --first_bp 1000000 \

? ? ? ? ? ? ? ? --last_bp 1000000 \

? ? ? ? ? ? ? ? --years_per_gen 28 \

? ? ? ? ? ? ? ? --coal popsize.coal

-i, --input? ? 輸入文件名（不包含擴(kuò)展名），應(yīng)為.anc和.mut文件的名稱。需要包含haps/sample輸入文件中的所有 SNP。

-o, --output? ? 輸出文件的基名（不包含擴(kuò)展名）。

-m, --mu? ? 突變率，例如1.25e-8。

--first_bp? ? 可選參數(shù)：起始?jí)A基對(duì)位置（bp）。

--last_bp? ? 可選參數(shù)：結(jié)束堿基對(duì)位置（bp）。

--years_per_gen? ? 可選參數(shù)：每代的年數(shù)。默認(rèn)值：28.0。

--coal? ? 可選參數(shù)：將所有單倍型視為一個(gè)群體時(shí)的共祖率估計(jì)。如果.anc/.mut文件是通過EstimatePopulationSize.sh腳本獲得的，則不需要此參數(shù)；否則強(qiáng)烈建議提供。

--seed? ? 可選參數(shù)：用于分支長(zhǎng)度估計(jì)的隨機(jī)種子。僅在指定--coal時(shí)使用。

輸出文件

.freq 文件：記錄衍生等位基因在genN到gen1世代（文件頭中指定）的頻率。

.lin 文件：記錄樹中譜系的數(shù)量，覆蓋genN到gen1世代（文件頭中指定），以及突變頻率為2時(shí)的譜系數(shù)量。

.sele 文件：記錄選擇證據(jù)的 log10 p 值，覆蓋genN到gen1世代（文件頭中指定），以及突變頻率為2時(shí)的 p 值。

獲得選擇的 P 值：

PATH_TO_RELATE/bin/RelateSelection \

? ? ? ? ? ? ? ? --mode Frequency \

? ? ? ? ? ? ? ? -i example \

? ? ? ? ? ? ? ? -o example

輸出文件為：example.freq, example.lin.

將上面的兩個(gè)文件作為輸入文件（example.freq, example.lin），計(jì)算每個(gè)SNP位點(diǎn)的P值：

PATH_TO_RELATE/bin/RelateSelection \

? ? ? ? ? ? ? ? --mode Selection \

? ? ? ? ? ? ? ? -i example \

? ? ? ? ? ? ? ? -o example

輸出文件為：example.sele。

7、推斷樣本支長(zhǎng)

PATH_TO_RELATE/scripts/SampleBranchLengths/SampleBranchLengths.sh \

? ? ? ? ? ? ? ? -i example \

? ? ? ? ? ? ? ? -o example_sub \

? ? ? ? ? ? ? ? -m 1.25e-8 \

? ? ? ? ? ? ? ? --coal popsize.coal \

? ? ? ? ? ? ? ? --format a \

? ? ? ? ? ? ? ? --num_samples 5 \

? ? ? ? ? ? ? ? --first_bp 100000 \

? ? ? ? ? ? ? ? --last_bp 200000 \

? ? ? ? ? ? ? ? --seed 1

-i, --input? ? 輸入文件名（不包含擴(kuò)展名），應(yīng)為.anc和.mut文件的名稱。需要包含輸入haps/sample文件中的所有 SNP。若未包含所有 SNP，請(qǐng)參考--dist選項(xiàng)。

-o, --output? ? 輸出文件的基名（不包含擴(kuò)展名）。

-m, --mu? ? 突變率，例如1.25e-8。

--coal? ? 估計(jì)共祖率（coalescent rate），將所有單倍型視為一個(gè)群體。

--format? ? 輸出文件格式：'a': anc/mut 格式（默認(rèn)）。'n': newick/sites 格式（供 CLUES 使用）。'b': 二進(jìn)制格式（供 CLUES 和 PALM 使用）注意：anc/mut 格式會(huì)修改以支持存儲(chǔ)同一分支的多個(gè)長(zhǎng)度。

--num_samples? ? 分支長(zhǎng)度采樣次數(shù)。需為整數(shù)且 ≥1。

--num_proposals? ? 可選參數(shù)：采樣間隔的 MCMC 提議次數(shù)。默認(rèn)值為max(10*N, 1000)（N 為單倍型數(shù)量）。若設(shè)為0，所有樣本將與輸入一致。

--first_bp? ? 可選參數(shù)：起始?jí)A基對(duì)位置（bp）。

--last_bp? ? 可選參數(shù)：結(jié)束堿基對(duì)位置（bp）。

--dist? ? 可選參數(shù)：包含 SNP 的堿基對(duì)位置及相鄰 SNP 間距的文件名（文件格式）。例如，通過RelateExtract --mode AncMutForSubregion生成。注意：若.anc/.mut文件僅覆蓋輸入haps/sample文件的子區(qū)域，則必須提供此文件。

--seed? ? 可選參數(shù)：分支長(zhǎng)度估計(jì)的隨機(jī)種子。

輸出文件為 example_sub.anc/mut 或 example_sub.newick/.sites 或 example_sub.timeb.

該步驟的輸出文件可以用于后續(xù)的 CLUES 或 PALM 分析，我后面有時(shí)間會(huì)寫一下如何利用 CLUES 推斷位點(diǎn)的受選擇時(shí)間。

詳見：CLUES 推斷位點(diǎn)的受選擇時(shí)間（進(jìn)化軌跡）

8、繪制樹圖

PATH_TO_RELATE/scripts/TreeView/TreeView.sh \

? ? ? ? ? ? ? ? --haps data/example.haps \

? ? ? ? ? ? ? ? --sample data/example.sample \

? ? ? ? ? ? ? ? --anc example.anc \

? ? ? ? ? ? ? ? --mut example.mut \

? ? ? ? ? ? ? ? --poplabels data/example.poplabels \

? ? ? ? ? ? ? ? --bp_of_interest 3000000 \

? ? ? ? ? ? ? ? --years_per_gen 28 \

? ? ? ? ? ? ? ? -o example

基于Relate計(jì)算的結(jié)果便可計(jì)算。

PATH_TO_RELATE/scripts/TreeView/TreeViewMutation.sh \

? ? ? ? ? ? ? ? ? --haps data/example.haps

? ? ? ? ? ? ? ? ? --sample data/example.sample \

? ? ? ? ? ? ? ? ? --anc example.anc \

? ? ? ? ? ? ? ? ? --mut example.mut \

? ? ? ? ? ? ? ? ? --poplabels data/example.poplabels \

? ? ? ? ? ? ? ? ? --bp_of_interest 3000000 \

? ? ? ? ? ? ? ? ? --years_per_gen 28 \

? ? ? ? ? ? ? ? ? -o example

在樹圖中，高亮顯示攜帶衍生等位的樣本。

PATH_TOcowplot/TreeView/TreeViewSamples.sh \

? ? ? ? ? ? ? ? --haps data/example.haps

? ? ? ? ? ? ? ? --sample data/example.sample \

? ? ? ? ? ? ? ? --anc example.anc \

? ? ? ? ? ? ? ? --mut example.mut \

? ? ? ? ? ? ? ? --dist example.dist \

? ? ? ? ? ? ? ? --poplabels data/example.poplabels \

? ? ? ? ? ? ? ? --bp_of_interest 3000000 \

? ? ? ? ? ? ? ? --years_per_gen 28 \

? ? ? ? ? ? ? ? -o example

這個(gè)要用 SampleBranchLengths 的輸出文件，在每個(gè)節(jié)點(diǎn)添加置信區(qū)間。