当前位置：首页 > news >正文

GWAS_LD

news 来源：原创 2025/4/21 16:56:23

局部LDblock 绘图

1. 查看显著位点附近基因情况

链接pvalue显著位点文件

 ln -s ~/yiyaoran/GWAS/my_GWAS_J/P3.GWAS/01.tassel/mlm_output.manht_figure.sigSite.out .
#也可以自己筛选awk '$2 == 9 && $4 <= 0.000028481' mlm_output.manht_input>368_GWAS.snp

snp两端延伸50000bp作为绘图窗口, 注意chr名称问题

awk '$3 ~ /[0-9]/{print $2"\t"$3-1"\t"$3}' 368_GWAS.snp | \
awk '{if($1<10){print "Chr0" $0 }else{print "Chr" $0 }}' > ./GWAS.snp.bed#准备gene_length文件  注意染色体名称
cp  ~/yiyaoran/01.ref/chrsize.txt .
awk 'NR > 1 {printf "Chr%02d\t%d\n", $1, $2}' chrsize.txt > new_chrsize.txt
#"Chr%02d\t%d\n"：
#Chr 是固定的前缀。
#%02d 表示将第一个字段（$1）格式化为两位数字，不足两位时前面补零（例如 1 变成 01）。
#\t 是制表符，用于分隔列。
#%d 表示第二个字段（$2）。
#\n 是换行符。#bedtools slop 用于扩展或收缩 BED 文件中的每个区间。
bedtools slop -b 50000 -i GWAS.snp.bed -g new_chrsize.txt > GWAS.snp.50000.bed#bedtools merge 用于合并重叠或相邻的区间
bedtools merge -i GWAS.snp.50000.bed > GWAS.snp.merge.bed

准备基因位置信息文件gff3或gtf文件

 awk '$3=="gene"' Zea_mays.Zm-B73-REFERENCE-NAM-5.0.60.gtf >gene.gtf

提取窗口内基因信息以及与显著snp最近的基因

注意bed文件的染色体名称和gtf文件的染色体名称是否一致

#修改染色体名称chr09为9
awk '{gsub(/^Chr/, ""); gsub(/^0+/, ""); print}' GWAS.snp.merge.bed > GWAS.snp.merge.fixed.bed
#提取窗口内的基因信息：
bedtools intersect -wo -a GWAS.snp.merge.bed -b gene.gtf > GWAS.snp.merge.gene
#bedtools intersect：用于找出两个文件中重叠的区间。
#-wo：输出重叠的区间，并在输出中包含重叠的长度（overlap size）#修改染色体名称
awk '{gsub(/^Chr/, ""); gsub(/^0+/, ""); print}' GWAS.snp.bed >GWAS.snp.fixed.bed
#找出与显著 SNP 最近的基因：
bedtools closest -a GWAS.snp.fixed.bed -b gene.sorted.gtf -d >GWAS.snp.gene
#bedtools closest：用于找出每个区间最近的非重叠区间。
#-d：在输出中包含最近区间的距离。