介紹：使用PopLDdecay軟件生成群體的LD衰減圖。從大群體中提取subgroup的vcf作為PopLDdecay的輸入。

#!/bin/bash

SamLst=$1   ##- 每個(gè)樣本名一行
Prefix=$2

gvcf=$3

vcftools  --gzvcf  ${gvcf}  --recode --recode-INFO-all --stdout  --keep  ${SamLst}  > ${Prefix}.snp.vcf

bgzip  ${Prefix}.snp.vcf

tabix -p vcf  ${Prefix}.snp.vcf.gz

/data/mg1/caix/src/biosoft/PopLDdecay/PopLDdecay  -InVCF  ${Prefix}.snp.vcf.gz  -OutStat  ${Prefix}.LD.stat.gz

rm ${Prefix}.snp.vcf.gz  ${Prefix}.snp.vcf.gz.tbi

使用PopLDdecay中程序Plot_MultiPop.pl繪圖

perl  /data/mg1/caix/src/biosoft/PopLDdecay/bin/Plot_MultiPop.pl  -inList  Stat.Sp.Lst  -output  oilseed  -keepR

Stat.Sp.Lst文件格式：

/data/mg1/caix/works/Bra_reSeqProject_2022/Pop_analysis_0411/oilseed/LD_decay/plot/CAturnips.LD.stat.gz CAoilseed
/data/mg1/caix/works/Bra_reSeqProject_2022/Pop_analysis_0411/oilseed/LD_decay/plot/Oilseed.LD.stat.gz   ChineseOilseed
/data/mg1/caix/works/Bra_reSeqProject_2022/Pop_analysis_0411/oilseed/LD_decay/plot/sarson.LD.stat.gz    sarson

結(jié)果圖展示：

image.png

#!/usr/bin/perl -w
use 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;
       } 

    }

}