檢測(cè)基因組選擇信號(hào)的方法有很多種，其中 XP-CLR 方法是常用的一種。XP-CLR 是陳華老師、Nick Patterson 和 David Reich 在 2010 年發(fā)表的方法，全稱叫 the cross-population composite likelihood ratio test（跨群體復(fù)合似然比檢驗(yàn)），是一種是基于選擇掃蕩（selective sweeep）的似然方法。
選擇掃蕩可以增加群體之間的遺傳分化，導(dǎo)致等位基因頻率偏離中性條件下的預(yù)期值。XP-CLR 利用了兩個(gè)群體之間的多基因座等位基因頻率差異（multilocus allele frequency differentiation）建立模型，使用布朗運(yùn)動(dòng)來(lái)模擬中性下的遺傳漂移，并使用確定性模型來(lái)近似地對(duì)附近的單核苷酸多態(tài)性（SNPs）進(jìn)行選擇性掃描。
python版本XP-CLR
python版本XP-CLR，具體可看 https://github.com/hardingnj/xpclr
1 安裝
lone this git repository into your working directory and cd.
python setup.py install
Also available via conda via the bioconda channel:
conda create -n xpclr -c bioconda xpclr
2 輸入文件
該軟件可接受輸入格式為:
(1) hdf5，通過(guò)scikit-allel 將vcf轉(zhuǎn)換成該格式，
(2) vcf格式
(3) 或者像之前xpclr中的3種類型格式
具體可以看xpclr-master/fixture/下的示例文件
- 兩個(gè)群體的gen 基因型文件，其組成如下,缺失可以用9代替
1 1 0 1 9 9
0 1 1 0 0 1
1 1 0 0 1 0
每2列為一個(gè)樣本的基因型，以此類推。
- snp, 位置信息，其格式如下
rs1081440 9 0.000109 36587 C T
rs9408752 9 0.000938 91857 A G
rs2810979 9 0.001323 152695 G A
每一列以此為：SNP編號(hào)(自行定義即可)，染色體，遺傳距離(可簡(jiǎn)單的物理距離/100000000)，SNP位置，Ref, Alt
3 運(yùn)行
基本參數(shù)說(shuō)明
xpclr -h
usage: xpclr [-h] --out OUT [--format FORMAT] [--input INPUT]
[--gdistkey GDISTKEY] [--samplesA SAMPLESA] [--samplesB SAMPLESB]
[--rrate RRATE] [--map MAP] [--popA POPA] [--popB POPB] --chr
CHROM [--ld LDCUTOFF] [--phased] [--verbose VERBOSE]
[--maxsnps MAXSNPS] [--minsnps MINSNPS] [--size SIZE]
[--start START] [--stop STOP] [--step STEP]
--out: 輸出文件
--format： 輸入格式，vcf，hdf5，zarr，txt（對(duì)應(yīng)2種基因型，和一個(gè)snp位點(diǎn)文件）
--input：輸入vcf，或者h(yuǎn)df5， 選擇txt時(shí)，不選擇，
--samplesA： 樣本A名稱文件； txt時(shí)，不選擇
--samplesB：樣本B名稱文件； txt時(shí)，不選擇
--map: 基因位點(diǎn)文件
--popA： 樣本A基因型文件，txt時(shí)使用
--popB： 樣本B基因型文件，txt時(shí)使用
--chr： 染色體，和輸入文件一致即可，每次分析一個(gè)
--ld： LD值，進(jìn)行權(quán)重分析
--maxsnps： 一個(gè)窗口最大的SNP
--minsnps： 一個(gè)窗口最小的SNP
--size： 窗口大小
--step： 步長(zhǎng)
其中，群體A為reference，群體B為目標(biāo)群體。
運(yùn)行示例文件
## 輸入VCF文件
xpclr --out test -Sa samplesA.txt -Sb samplesB.txt \
-I small.vcf.gz -C 3L --ld 0.95 --phased --maxsnps 600 \
--size 200000 --step 20000
### 輸入txt文件
xpclr --format txt --out test --map mapfile.snp \
--popA genotype1.geno --popB genotype2.geno \
--chr 1 --ld 0.95 --phased --maxsnps 600 \
--size 200000 --step 20000
結(jié)果
結(jié)果文件中倒數(shù)兩列分別為xpclr標(biāo)準(zhǔn)化值，以及xpclr的值
原版XP-CLR
原版XP-CLR多年沒(méi)有更新
1 安裝
軟件可以從此處下載https://reich.hms.harvard.edu/software
wget https://reich.hms.harvard.edu/sites/reich.hms.harvard.edu/files/inline-files/XPCLR.tar
tar XPCLR.tar
在bin目錄下有執(zhí)行程序，以及3個(gè)示例文件
如果 XPCLR 無(wú)法在你的系統(tǒng)中運(yùn)行，則需要自己用 src 的源碼編譯：
cd src
make
make install
2 參數(shù)說(shuō)明
./XPCLR -h
Usage:
XPCLR -xpclr hapmapInput1 hapmapInput2 mapInput outFile \
-w1 gWin(Morgan) snpWin gridSize(bp) chrN -p corrLevel
# -xpclr: 后面接是兩個(gè)群體的 .geno 文件（genofile1 和 genofile2）、 .snp 文件（mapfile）、輸出文件（outputFile）
# -w1: 后面接的參數(shù)依次為：gWin 是 Morgan 為單位的window size (一般可以設(shè)為 0.005)；snpWin 代表一個(gè) window 中最大的 SNP 數(shù)量；gridSize 是 bp 為單位的兩個(gè) grid points 的間距；chrN 是染色體數(shù)
# p: -p1 代表 phased 的數(shù)據(jù)，-p0 代表未 phased; 如果存在2個(gè)SNP，A/G, G/C，不能明確A，G或者A，C為同一染色體，則是unphased
#corrLevel：加權(quán)復(fù)合似然比檢驗(yàn)中的 criterion，一般可設(shè)為0.95
3 運(yùn)行示例數(shù)據(jù)
/data/pub/shehb/soft/XPCLR/bin/XPCLR -xpclr CEU.9 YRI.9 9.xpclr.b36.snp new -w1 0.005 200 1000 9 -p0 0.95
結(jié)果文件沒(méi)一列依次為：chr, grid, ofSNPs_in_window, physical_pos, genetic_pos, XPCLR_score max_s.
4 其它說(shuō)明
如果在運(yùn)行python版本的過(guò)程當(dāng)中，出現(xiàn)如下錯(cuò)誤
No permission to write in the specified directory: {0}".format(outdir
則找到對(duì)應(yīng)xpclr文件，加入
fn = args.out
fn = os.path.abspath(args.out)
outdir = os.path.dirname(fn)
- snp位置信息文件的制作
過(guò)濾后的vcf文件
zcat in.vcf.gz | awk '($1=="MC01"){print $1":"$2"\tMC01\t"$2/100000000"\t"$2"\t"$4"\t"$5}' |grep -v '#' >MC01.map
- 兩個(gè)genotype文件的制作
首先準(zhǔn)備群體名稱文件，相同的兩列
head genoe1.txt
MB MB
MF MF
S001 S001
S002 S002
S003 S003
S004 S004
S005 S005
S006 S006
然后利用plink進(jìn)行調(diào)去相應(yīng)序列
plink --vcf in.vcf.gz --keep genoe1.txt --chr MC01 \
--out MC01.out --recode 01 transpose \
-output-missing-genotype 9 --allow-extra-chr \
--set-missing-var-ids @:# --keep-allele-order
## -output-missing-genotype； 不符合就定義為9
得到一個(gè).tped文件，然后調(diào)取對(duì)應(yīng)基因型即可
cut -d " " -f 5- MC01.out.tped >popA_MC01
