#!/bin/bash
#LD PLOT
#author lee

echo"par1 is vcf,par2 is top_chr"
Genome=/public/home/caisl/lee/genome/rice/msu.gff3

while read chr pos
do

pos_up=$((pos-20000))
pos_down=$((pos+20000))
region=$chr":"$pos_up":"$pos_down
LDBlockShow -InVCF $1 -OutPut ${chr}_${pos} -Region $region -InGFF $Genome -OutPng -SeleVar 2
gunzip ${chr}_${pos}.blocks.gz
grep $pos ${chr}_${pos}.blocks|awk '{print}'>${chr}_${pos}_block
start=$(awk '{print $2}' ${chr}_${pos}_block)
end=$(awk '{print $3}' ${chr}_${pos}_block)
awk -v chr=$chr -v start=$start -v end=$end '{if($3=="gene" && $1==chr && $4>=start && $4<=end)print $0}' $Genome>${chr}_${pos}_block_gene;
sed -i 's/|/\t/g' ${chr}_${pos}_block
awk '{for(i=6;i<=NF;i++) a[i,NR]=$i}END{for(i=6;i<=NF;i++) {for(j=1;j<=NR;j++) printf a[i,j] " ";print $1}}' ${chr}_${pos}_block>${chr}_${pos}_block_snps_in_LD
awk '{print $2,$1}' ${chr}_${pos}_block_snps_in_LD>${chr}_${pos}_block_snps_in_LD2
done<$2