From 72b81e3c8abdd6fa84ad6fe72f30b7d5e306e1be Mon Sep 17 00:00:00 2001 From: ixxmu Date: Tue, 28 May 2024 10:57:08 +0000 Subject: [PATCH] A custom message for the commit --- ...Dblock\347\232\204\350\204\232\346\234\254.md" | 15 --------------- 1 file changed, 15 deletions(-) delete mode 100644 "docs/2024-05/\344\270\200\344\270\252\346\211\271\351\207\217\345\201\232LDblock\347\232\204\350\204\232\346\234\254.md" diff --git "a/docs/2024-05/\344\270\200\344\270\252\346\211\271\351\207\217\345\201\232LDblock\347\232\204\350\204\232\346\234\254.md" "b/docs/2024-05/\344\270\200\344\270\252\346\211\271\351\207\217\345\201\232LDblock\347\232\204\350\204\232\346\234\254.md" deleted file mode 100644 index c1bb2436..00000000 --- "a/docs/2024-05/\344\270\200\344\270\252\346\211\271\351\207\217\345\201\232LDblock\347\232\204\350\204\232\346\234\254.md" +++ /dev/null @@ -1,15 +0,0 @@ ---- -title: "一个批量做LDblock的脚本" -date: 2024-05-28T10:56:05Z -draft: ["false"] -tags: [ - "fetched", - "育种数据分析之放飞自我" -] -categories: ["Acdemic"] ---- -一个批量做LDblock的脚本 by 育种数据分析之放飞自我 ------- -

大家好，我是邓飞。

今天，有小伙伴在星球询问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 1fi
vcf=$1;chr=$2;pos=$3;qujian=$4st=$(( $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 2807198827 1435474132 1419428142 1419427351 17550492610 336195821 24513624410 1340346861 540492201 2852738455 182405365 832334878 1232729356 9326360510 14417231610 35982628 47665935 2008132768 1605363001 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 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 7 143547413 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 2 141942814 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 2 141942735 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 1 175504926 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 10 33619582 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 1 245136244 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 10 134034686 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 1 54049220 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 1 285273845 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 5 18240536 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 5 83233487 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 8 123272935 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 6 93263605 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 10 144172316 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 10 3598262 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 8 4766593 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 5 200813276 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 8 160536300 100000bash plot_ld_block_vcf.sh ../geno/genotype.vcf 1 286271514 100000

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

bash temp.sh

静等结果就可以了。

666

推荐阅读：


想要更好的学习和交流，快来加入飞哥的知识星球，这是一个生物统计+数量遗传学+GWAS+GS的社区，在这里你可以向飞哥提问、帮你指定学习计划、跟着飞哥一起做实战项目，冲冲冲。点击这里加入吧：飞哥的学习圈子


1，快来领取 | 飞哥的GWAS分析教程


2，飞哥汇总 | 入门数据分析资源推荐


3，数量遗传学，分享几本书的电子版


4，R语言学习看最新版的电子书不香嘛？


5，书籍及配套代码领取--统计遗传分析导论

-
-原文链接