|
大家好，我是鄧飛。 LD衰減圖，可以形象的查看群體LD衰減的情況。LD衰減是由于連鎖不平衡所致，LD衰減速度在不同物種或者不同亞種中差異不同，通常用LD衰減到一般的距離來作為群體的衰減距離（還有其它計算方法），如果LD衰減很快，則在進(jìn)行GWAS分析時需要更多的位點(diǎn)才能達(dá)到一定的精度。（計算群體GWAS分析所需要的最少SNP個數(shù)） 另外，LD衰減也可以反應(yīng)群體受選擇的情況，一般來說，野生群體比馴化改良群體LD衰減快，異花授粉比自花授粉植物L(fēng)D衰減快。 之前寫過推文教程（LD衰減圖繪制--PopLDdecay） 
出的圖是上面這個樣子的，如果人為查看的衰減一半的距離，大約是100kb左右，如何更科學(xué)的計算呢？ 網(wǎng)上看到了一個perl腳本可以根據(jù)PopLDdecay的結(jié)果自動計算衰減一半的距離：（https://www.jianshu.com/p/8205dbcb3839） 代碼如下：calculate_LDlength.pl #!/usr/bin/perl -wuse strict;
my $in0 = $ARGV[0]; ##- sarson.LD.stat.gz
open IN0, "gzip -dc $in0 | ";<IN0>;my $firstLine = <IN0>;chomp($firstLine);my @firstLine = split(/\t/, $firstLine);my $max = $firstLine[1];close IN0;
my %dis2Value = ();open IN1, "gzip -dc $in0 | ";<IN1>;while(<IN1>){ chomp; my @temp = split(/\t/, $_); $dis2Value{$temp[0]} = $temp[1];}close IN1;
my $halfValue = $max/2;
for my $key1(sort {$a<=>$b} keys %dis2Value){
my $next = $key1 + 1; if(exists $dis2Value{$next}) {
my $currentValue = $dis2Value{$key1}; my $nextValue = $dis2Value{$next}; if($currentValue >= $halfValue && $nextValue < $halfValue){ print "Processing ", $in0, "\n"; print "max LD: r2: ", $max, "\n"; print "half LD: r2: ", $halfValue, "\t", "LD length: ", $key1, "\n"; last; }
}
}
例如PopLDdecay生成的文件為：LDdecay.stat.gz，用上面的程序處理這個文件，就能得到衰減一半的距離，調(diào)用方法： $ perl calculate_LDlength.pl LDdecay.stat.gz Processing LDdecay.stat.gzmax LD: r2: 0.8551half LD: r2: 0.42755 LD length: 67
可以看到LD衰減一半的值是0.427，對應(yīng)的距離是67kb。 以上。
|
