ASMC(ascertained sequentiallyMarkovian coalescent)是2018年Alkes L. Price實(shí)驗(yàn)室開(kāi)發(fā)的算法，相關(guān)文章發(fā)表在2018年的Nature Genetics上制市。
該方法利用 ASMC允粤，分析SNP芯片和 WGS數(shù)據(jù)集萍聊，目的是快速計(jì)算基因組位點(diǎn)的溯祖時(shí)間问芬，從而可以檢測(cè)正選擇和背景選擇的特征。
1. Installation
首先需要配置一下 C++ 庫(kù)寿桨，基本上conda都可以裝此衅。
C++ compiler (C++17 or later)
CMake (3.15 or later)
Boost (1.62 or later)
Eigen (3.3.4 or later)
{fmt}
range-v3
OpenMP
zlib
2. 編譯ASMC、FastSMC和二進(jìn)制轉(zhuǎn)換可執(zhí)行文件
# Get the source
git clone --recurse-submodules https://github.com/PalamaraLab/ASMC
cd ASMC
# Create a build directory
mkdir build && cd build
# Configure and build
# On first run, CMake will build the required dependencies
cmake -DASMC_NO_PYTHON=TRUE ..
cmake --build . --parallel 4
#Tool to compute decoding quantities.
pip install asmc-preparedecoding
3. 輸入輸出文件格式
（1）輸入文件
Genetic map (.map/map.gz)
可以使用plink來(lái)將vcf轉(zhuǎn)換成map格式
plink --vcf input.vcf --recode --out output --double-id --autosome-num 42
Demographic history (.demo)
該文件包含了群體過(guò)去的有效規(guī)模的分段大小亭螟。包含兩列：
TimeStart PopulationSize
其中 TimeStart是第一代人口的大小為 PopulationSize挡鞍。 第一行代表第 0 代，可以使用PSMC/MSMC/SMC++獲得該文件數(shù)據(jù)预烙。
Time discretization (.disc)
在 ASMC 的輸入中提供的時(shí)間間隔列表墨微。每行包含一個(gè)數(shù)字，表示以連續(xù)代測(cè)量的時(shí)間扁掸，并從 0.0 代開(kāi)始翘县。
0.0
30.0
60.0
90.0
120.0
150.0
180.0
... <lines omitted>
79855.6
96263.0
124311.7
此文件中定義的間隔為：{0.0-30.0, 30.0-60.0, ..., 96263.0-124311.7, 124311.7-Infinity}。
（2）輸出文件
Decodingquantities(.decodingQuantities.gz)
這個(gè)文件是ASMCprepareDecoding的輸出文件谴分，是ASMC的輸入文件锈麸，它用于執(zhí)行溯祖時(shí)間的有效推斷，不需要了解該文件的內(nèi)容牺蹄。
Time discretization intervals (.intervalsInfo)
*.intervalsInfo文件是 ASMCprepareDecoding的輸出文件忘伞，作為ASMC的輸入文件。它包含分析中使用的離散時(shí)間間隔數(shù)相對(duì)應(yīng)的行數(shù)。 以下是文件格式：
IntervalStart ExpectedCoalescenceTime IntervalEnd
IntervalStart 和 IntervalEnd 表示每個(gè)離散時(shí)間間隔的開(kāi)始/結(jié)束氓奈，ExpectedCoalescenceTime 是已推斷在此時(shí)間間隔內(nèi)溯祖的一對(duì)個(gè)體的預(yù)期溯祖時(shí)間匿刮，并且取決于人口統(tǒng)計(jì)模型。
成對(duì)后驗(yàn)合并概率之和.sumOverPairs.gz
ASMC的輸出寫(xiě)入 *.{00,01,11}.sumOverPairs.gz文件探颈。 每個(gè)文件都包含一個(gè)大小為 SxD 的矩陣，其中 S 是數(shù)據(jù)中的位點(diǎn)數(shù)训措，D 是分析中使用的離散時(shí)間間隔數(shù)伪节。 每個(gè)矩陣的第 [i,j] 個(gè)條目包含 SNP i 和離散合并時(shí)間 j 處所有樣本的后驗(yàn)溯祖概率之和。 輸出使用 {00,01,11} 后綴分解每個(gè)位點(diǎn)攜帶不同等位基因的樣本的合并事件绩鸣。 具體來(lái)說(shuō)：
*.00.sumOverPairs.gz 中矩陣的第i行包含在位點(diǎn)i處為純合0的所有樣本對(duì)的后驗(yàn)概率之和怀大。
*.01.sumOverPairs.gz 中矩陣的第 i 行包含在位點(diǎn) i 雜合的所有樣本對(duì)的后驗(yàn)概率之和。
*.11.sumOverPairs.gz 中矩陣的第 i 行包含在位點(diǎn) i 處為純合 1 的所有樣本對(duì)的后驗(yàn)概率之和呀闻。
4. 運(yùn)行
（1）Creating decoding quantities
首先需要使用已知的群體歷史（*.demo文件化借，離散Ne值），設(shè)定的離散時(shí)間間隔（.disc文件)捡多，以及單群體的等位基因頻率文件（ *frq文件蓖康，可以用vcftools計(jì)算得到），計(jì)算解碼文件垒手。
在python里面計(jì)算：
import pathlib
from asmc.preparedecoding import *
files_dir = (pathlib.Path('..') / 'test' / 'regression').resolve()
demo_file = str(files_dir / 'input_CEU.demo')
disc_file = str(files_dir / 'input_30-100-2000.disc')
freq_file = str(files_dir / 'input_UKBB.frq')
# Calculate with a discretization file
dq = prepare_decoding(
demography=demo_file,
discretization=disc_file,
frequencies=freq_file,
samples=50, # Use a larger number (300 is suggested) for real analysis
)
# Or calculate with pre-defined quantiles, and a user-defined number of additional quantiles
dq = prepare_decoding(
demography=demo_file,
discretization=[[30.0, 15], [100.0, 15], 39],
frequencies=freq_file,
samples=50, # Use a larger number (300 is suggested) for real analysis
)
# Or use with built-in demographies
dq = prepare_decoding(
demography='CEU',
discretization=[[30.0, 15], [100.0, 15], 39],
frequencies=freq_file,
samples=50, # Use a larger number (300 is suggested) for real analysis
)
# Or use with built-in frequency information from UKBB
dq = prepare_decoding(
demography='CEU',
discretization=[[30.0, 15], [100.0, 15], 39],
frequencies='UKBB',
samples=50, # Use a larger number (300 is suggested) for real analysis
)
#寫(xiě)出
# This will create a file `files_dir/output.decodingQuantities.gz`
dq.save_decoding_quantities(str(files_dir / 'output'))
# This will create a file `files_dir/output.csfs`
dq.save_csfs(str(files_dir / 'output'))
# You may also save other files, which may be of use:
dq.save_intervals(str(files_dir / 'output'))
dq.save_discretization(str(files_dir / 'output'))
save_demography(str(files_dir), 'CEU')
（2）運(yùn)行ASMC
./build/ASMC_exe \
--decodingQuantFile ASMC_data/decoding_quantities/30-100-2000_CEU.decodingQuantities.gz \
--inFileRoot ASMC_data/examples/asmc/exampleFile.n300 \
--majorMinorPosteriorSums \
--mode sequence
參數(shù)解讀：
必須參數(shù)：
--inFileRoot arg haps|haps|hap.gz|haps.gz 和 sample|samples 文件的前綴
--decodingQuantFile arg 解碼文件
以下選其一：
--posteriorSums 后驗(yàn)分布總和的輸出文件蒜焊。
--majorMinorPosteriorSums 后驗(yàn)分布總和的輸出文件，分為主要/次要等位基因科贬。
可選參數(shù)：
--outFileRoot arg 后驗(yàn)分布總和的輸出文件前綴（默認(rèn)值：--hapsFileRoot 參數(shù)）
--jobs int (=0) 并行完成的任務(wù)數(shù)
--jobInd int (=0) 任務(wù)索引 (1..jobs)
--mode string (=array) 解碼模式泳梆。可選 {sequence, array} 榜掌。
--compress (=false) 壓縮為二進(jìn)制（無(wú) CSFS）
--useAncestral (=false) 假設(shè)祖先等位基因編碼為 1 （否則將假設(shè) 1 = 次要等位）
--skipCSFSdistance float (=0.0) 兩個(gè) CSFS 排放之間的遺傳距離（以 cM 為單位）
5. 結(jié)果可視化
gunzip -c Alineage.1-1.00.sumOverPairs.gz | \
awk '{for (i=5; i<=NF; i++) {c[i]+=$i; t+=$i;}} END{for (i=5; i<=NF; i++) print c[i]/t}' \
> norm.txt
factor=1.5
cat /ASMC/input/gene_regions.txt | \
tr '-' '\t' | \
awk -v factor=$factor '{
l=$3-$2;
chr=$1;
from=$2-factor*l;
to=$3+factor*l;
file="region"NR".chr"chr"."from"-"to;
print "python3 plotPosteriorHeatMap.py -i /ASMC/decoding_quantities/Alineage.intervalsInfo -t 1 18 -n norm.txt -o " file " -f "$1".txt.gz" " -p " ($2-2*l) " " ($3+2*l)
}' | \
while read l; do
$l;
done
其中优妙，gene_regions.txt文件表示想要展示的基因組區(qū)域，第一列是染色體ID憎账。如下格式：
1 1.42207-11.4235
1 33.7288-43.7335
1 37.104-47.1191
plotPosteriorHeatMap.py是ASMC自帶的腳本套硼。
6. Density of recent coalescence (DRC) statistic
如何通過(guò)計(jì)算ASMC來(lái)確定最近受選擇區(qū)域呢？作者沒(méi)有給出具體的代碼鼠哥，只給了分析思路熟菲。
可以通過(guò)將 sumOverPairs 矩陣中與所需時(shí)間間隔相對(duì)應(yīng)的條目相加，并在一個(gè)窗口內(nèi)對(duì) SNP 進(jìn)行平均來(lái)獲得DRC 統(tǒng)計(jì)數(shù)據(jù)朴恳。
關(guān)于分析 DRC 統(tǒng)計(jì)的幾點(diǎn)說(shuō)明：使用 Gamma 分布作為 DRC 統(tǒng)計(jì)的零模型是一種近似值抄罕，需要進(jìn)一步的工作來(lái)推導(dǎo)一個(gè)精確的零模型。DRC 不能大于 1于颖，因此在 Gamma 模型下為非常大的 DRC 值獲得的近似 p 值將非常写艋摺（假陽(yáng)性）。
作者發(fā)現(xiàn) Gamma 近似在某些條件下是合理的，例如最近的有效群體很大做入，這使得 DRC 在中性地區(qū)變小冒晰，但這種近似在所有人群（例如孤立人群）或時(shí)間間隔比較大時(shí)都不能很好地工作。
如果您決定使用這種方法來(lái)檢測(cè)候選區(qū)域以進(jìn)行選擇竟块，作者建議檢查 Gamma 近似值是否合理壶运。由于特別大的 DRC 值可能會(huì)導(dǎo)致人為地減小 p 值±嗣兀或者可以因此僅使用 Gamma 近似來(lái)確定近似的顯著性閾值蒋情，然后輸出大于此閾值的原始 DRC 值，而不是輸出近似的 p 值耸携。
最后棵癣，作者建議檢查以這種方式檢測(cè)到的候選區(qū)域是否與基因組的問(wèn)題區(qū)域不重疊，例如重組率非常高/低或 LD 的區(qū)域夺衍，以及可能影響的大結(jié)構(gòu)變異狈谊，例如倒位重組。
參考：
- P. Palamara, J. Terhorst, Y. Song, A. Price. High-throughput inference of pairwise coalescence times identifies signals of selection and enriched disease heritability. Nature Genetics, 2018.
- J. Nait Saada, G. Kalantzis, D. Shyr, F. Cooper, M. Robinson, A. Gusev, P. F. Palamara. Identity-by-descent detection across 487,409 British samples reveals fine-scale evolutionary history and trait associations. Nature Communications, 2020.
- Spence, J.P. and Song, Y.S. Inference and analysis of population-specific fine-scale recombination maps across 26 diverse human populations. Science Advances, Vol. 5, No. 10, eaaw9206 (2019)
- ASMC github
