一个批量做LDblock的脚本

邓飞

发布于 2024-05-29 14:50:06

1000

发布于 2024-05-29 14:50:06

文章被收录于专栏：育种数据分析之放飞自我育种数据分析之放飞自我

大家好，我是邓飞。

今天，有小伙伴在星球询问GWAS分析的显著性位点，如何批量绘制LDblock，之前写过LDblock的安装和使用说明：（LDblock绘制连锁不平衡和单体型图）

问题有3个核心点：

1，snp文件总是有^M符号，怎么去除

2，如何批量对显著位点进行LDblock分析

3，有没有现成代码

今天就搞定这个问题。

首先，看一下正常LDblock分析的代码：

LDBlockShow -InVCF hebing1.vcf -Region 3:400986543:401186543 -OutPng -SeleVar 1 --OutPut re3

有vcf文件，染色体，物理位置区间，结果文件名称。

来一个脚本：

$ cat ~/bin/plot_ld_plot_vcf.sh 
#!/bin/bash

if [ $# != 4 ]
then
    echo $0 name.vcf chronome_number location 20000; 第一个是vcf或者vcf.gz文件，第二个是染色体位置，第三个是物理位置，第四个是上下游区间
    exit 1
fi

vcf=$1;chr=$2;pos=$3;qujian=$4
st=$(( $pos - $qujian))
en=$(( $pos + qujian))

ot=re_${chr}_${pos}
echo $chr $pos $pos >temp.txt
~/bin/LDBlockShow -InVCF $vcf -Region ${chr}:${st}:${en} -SeleVar 1 -BlockType 2 -NoShowLDist 10000000 -OutPut ${ot}
~/bin/ShowLDSVG -InPreFix ${ot} -SpeSNPName temp.txt -OutPut ${ot}_a

上面的脚本，给定vcf文件，染色体，物理位置，上下游区间，就会自动分析。

如何批量处理呢？

比如snp.txt文件，如下图所示：

1 280719882
7 143547413
2 141942814
2 141942735
1 175504926
10 33619582
1 245136244
10 134034686
1 54049220
1 285273845
5 18240536
5 83233487
8 123272935
6 93263605
10 144172316
10 3598262
8 4766593
5 200813276
8 160536300
1 286271514

上下游区间是100kb，那么先将每个snp生成一个命令行，然后运行就行了。

awk走起：

awk '{print "bash plot_ld_block_vcf.sh ../geno/genotype.vcf", $1,$2,100000}' snp.txt >temp.sh

生成的脚本文件：

$ cat temp.sh 
bash plot_ld_block_vcf.sh ../geno/genotype.vcf 1 280719882 100000
bash plot_ld_block_vcf.sh ../geno/genotype.vcf 7 143547413 100000
bash plot_ld_block_vcf.sh ../geno/genotype.vcf 2 141942814 100000
bash plot_ld_block_vcf.sh ../geno/genotype.vcf 2 141942735 100000
bash plot_ld_block_vcf.sh ../geno/genotype.vcf 1 175504926 100000
bash plot_ld_block_vcf.sh ../geno/genotype.vcf 10 33619582 100000
bash plot_ld_block_vcf.sh ../geno/genotype.vcf 1 245136244 100000
bash plot_ld_block_vcf.sh ../geno/genotype.vcf 10 134034686 100000
bash plot_ld_block_vcf.sh ../geno/genotype.vcf 1 54049220 100000
bash plot_ld_block_vcf.sh ../geno/genotype.vcf 1 285273845 100000
bash plot_ld_block_vcf.sh ../geno/genotype.vcf 5 18240536 100000
bash plot_ld_block_vcf.sh ../geno/genotype.vcf 5 83233487 100000
bash plot_ld_block_vcf.sh ../geno/genotype.vcf 8 123272935 100000
bash plot_ld_block_vcf.sh ../geno/genotype.vcf 6 93263605 100000
bash plot_ld_block_vcf.sh ../geno/genotype.vcf 10 144172316 100000
bash plot_ld_block_vcf.sh ../geno/genotype.vcf 10 3598262 100000
bash plot_ld_block_vcf.sh ../geno/genotype.vcf 8 4766593 100000
bash plot_ld_block_vcf.sh ../geno/genotype.vcf 5 200813276 100000
bash plot_ld_block_vcf.sh ../geno/genotype.vcf 8 160536300 100000
bash plot_ld_block_vcf.sh ../geno/genotype.vcf 1 286271514 100000

然后运行temp.sh文件就可以了：