大家好,我是鄧飛。
LD衰減圖,可以形象的查看群體LD衰減的情況。LD衰減是由于連鎖不平衡所致,LD衰減速度在不同物種或者不同亞種中差異不同,通常用LD衰減到一般的距離來(lái)作為群體的衰減距離(還有其它計(jì)算方法),如果LD衰減很快,則在進(jìn)行GWAS分析時(shí)需要更多的位點(diǎn)才能達(dá)到一定的精度。(計(jì)算群體GWAS分析所需要的最少SNP個(gè)數(shù))
另外,LD衰減也可以反應(yīng)群體受選擇的情況,一般來(lái)說(shuō),野生群體比馴化改良群體LD衰減快,異花授粉比自花授粉植物L(fēng)D衰減快。
LD圖分為單個(gè)群體和多個(gè)群體的圖,今天使用軟件PopLDdecay來(lái)實(shí)現(xiàn)一下。這個(gè)軟件需要Linux系統(tǒng)。
目標(biāo)圖:
基因型數(shù)據(jù)是vcf數(shù)據(jù),plink格式的數(shù)據(jù)可以轉(zhuǎn)化為vcf數(shù)據(jù)。
$ ls hebing.ped hebing.map
hebing.map hebing.ped
上面示例中用的是plink數(shù)據(jù),這里轉(zhuǎn)化為vcf格式的數(shù)據(jù):
plink --file hebing --recode vcf-iid --allow-extra-chr --chr-set 40 --out file
結(jié)果:
$ ls file.vcf
file.vcf
官網(wǎng):https://github.com/BGI-shenzhen/PopLDdecay
安裝方法:
git clone https://github.com/hewm2008/PopLDdecay.git
cd PopLDdecay; chmod 755 configure; ./configure;
make;
mv PopLDdecay bin/; # [rm *.o]
測(cè)試:
$ PopLDdecay
Usage: PopLDDecay -InVCF <in.vcf.gz> -OutStat <out.stat>
-InVCF <str> Input SNP VCF Format
-InGenotype <str> Input SNP Genotype Format
-OutStat <str> OutPut Stat Dist ~ r^2/D' File
-SubPop <str> SubGroup Sample File List[ALLsample]
-MaxDist <int> Max Distance (kb) between two SNP [300]
-MAF <float> Min minor allele frequency filter [0.005]
-Het <float> Max ratio of het allele filter [0.88]
-Miss <float> Max ratio of miss allele filter [0.25]
-EHH <str> To Run EHH Region decay set StartSite [NA]
-OutFilterSNP OutPut the final SNP to calculate
-OutType <int> 1: R^2 result 2: R^2 & D' result 3:PairWise LD Out[1]
See the Help for more OutType [1-8] details
-help Show more help [hewm2008 v3.40]
顯示安裝完成
~/bin/PopLDdecay -InVCF file.vcf -OutStat LDdecay
~/bin/Plot_OnePop.pl -inFile LDdecay.stat.gz --output re1 -bin1 10 -bin2 100
這里面,我用的是SNP數(shù)據(jù)位點(diǎn)較少,所以曲線不平滑,可以修改bin2=500,在進(jìn)行繪圖:
~/bin/Plot_OnePop.pl -inFile LDdecay.stat.gz --output re2 -bin1 10 -bin2 500
如果要運(yùn)行A,B,C三個(gè)品種的LD圖,先要把ID整理為三個(gè)文件:
文件:
$ ls *txt
A.txt B.txt C.txt total_name.txt
分別針對(duì)于每個(gè)品種,計(jì)算:
~/bin/PopLDdecay -InVCF file.vcf -OutStat A.stat.gz -SubPop A.txt
~/bin/PopLDdecay -InVCF file.vcf -OutStat B.stat.gz -SubPop B.txt
~/bin/PopLDdecay -InVCF file.vcf -OutStat C.stat.gz -SubPop C.txt
結(jié)果文件:
$ ls *stat.gz
A.stat.gz B.stat.gz C.stat.gz
生成multi.list,格式為兩列,第一列為文件名稱,第二列為品種名稱,內(nèi)容為:
$ cat multi.list
A.stat.gz A
B.stat.gz B
C.stat.gz C
作圖:
perl ~/bin/Plot_MultiPop.pl -inList multi.list --output re3 -bin1 10 -bin2 500
聯(lián)系客服