跟着Nature学数据分析:plink计算SNP和SV之间的连锁不平衡R方值

论文

Graph pangenome captures missing heritability and empowers tomato breeding
论文链接
https://www.nature.com/articles/s41586-022-04808-9
这篇论文里对应的代码https://github.com/YaoZhou89/TGG
在代码部分并没有找到关于计算ld的代码，论文中也没有找到相关方法的描述。论文中提供了SNP Indel 和 sv数据集。下载下来自己算算试试
数据下载链接http://solomics.agis.org.cn/tomato/ftp/
snp indel 数据集 只下载 chr3的部分
SV数据集的处理

sv的数据集把3号染色体过滤出来
bcftools view 706.sv.vcf.gz -r 3 -O v -o chr3.sv.vcf
自己写一个python脚本修改一些vcf文件里的内容
把id 改成 chr + pos + “_SV”的形式，把INFO列的内容都去掉，把 alt 和 ref 都改成 单碱基的形式 基因型只保留前三个字符
python 20240524_01.py chr3.sv.vcf chr3.sv.edited.vcf
20240524_01.py脚本的内容
import sys

fw = open(sys.argv[2],’w’)

with open(sys.argv[1],’r’) as fr:
for line in fr:
if line.startswith(“#”):
fw.write(line)
else:
temp_list = line.strip().split(“\t”)
new_list = [None]*len(temp_list)

        new_list[0] = temp_list[0]
        new_list[1] = temp_list[1]
        new_list[2] = temp_list[0] + "_" + temp_list[1] + "_SV"
        new_list[3] = "A"
        new_list[4] = "T"
        new_list[5] = temp_list[5]
        new_list[6] = temp_list[6]
        new_list[7] = "."
        new_list[8] = "GT"

        new_list[9:] = ['./.' if j.count(".") > 0 else j for j in [i[0:3] for i in temp_list[9:]]]

        fw.write("%s\n"%('\t'.join(new_list)))

fw.close()
这个vcf文件里不知道为啥会有很多 .:1 .:0这种基因型，如果是这种同意替换成./.
根据缺失率对数据进行过滤
vcftools –vcf chr3.sv.edited.vcf –max-missing 0.8 –recode –recode-INFO-all –out chr3.sv.edited.filter
还剩12688个位点
做一个基因型填充
beagle gt=chr3.sv.edited.filter.recode.vcf out=chr3.sv.edited.filter.impute nthreads=24
SNP INDEL数据集的处理

只保留二等位的位点，缺失率
time vcftools –gzvcf chr3.vcf.gz –maf 0.05 –min-alleles 2 –max-alleles 2 –max-missing 0.8 –recode –recode-INFO-all –out chr3.snp.filter
549904 out of a possible 2827014 Sites
15Min
基因型填充
time java -Xmx48g -jar ~/anaconda3/envs/syri/share/beagle-5.2_21Apr21.304-0/beagle.jar gt=chr3.snp.filter.recode.vcf out=chr3.snp.filter.impute nthreads=24
7min
把填充后的id改一下
bgzip -d chr3.snp.filter.impute.vcf.gz
python 20240524_02.py chr3.snp.filter.impute.vcf chr3.snp.filter.impute.edited.vcf
脚本内容
import sys

fw = open(sys.argv[2],’w’)

with open(sys.argv[1],’r’) as fr:
for line in fr:
if line.startswith(“#”):
fw.write(line)
else:
temp_list = line.strip().split(“\t”)

        if len(temp_list[3]) == len(temp_list[4]):
            temp_list[2] = temp_list[0] + "_" + temp_list[1] + "_SNP"
        else:
            temp_list[2] = temp_list[0] + "_" + temp_list[1] + "_INDEL"
            temp_list[3] = "A"
            temp_list[4] = "T"

        fw.write("%s\n"%('\t'.join(temp_list)))

fw.close()
把SNP INDEL和SV数据集合并

bgzip chr3.snp.filter.impute.edited.vcf
tabix chr3.snp.filter.impute.edited.vcf.gz

tabix chr3.sv.edited.filter.impute.vcf.gz

vcfcat chr3.sv.edited.filter.impute.vcf.gz chr3.snp.filter.impute.edited.vcf.gz > merged.sv.snp.vcf

vcfsort merged.sv.snp.vcf > merged.sv.snp.sorted.vcf

计算ld R2

参考链接
https://speciationgenomics.github.io/ld_decay/
这里介绍的还挺详细的
plink –vcf merged.sv.snp.sorted.vcf \
–double-id \
–allow-extra-chr \
–maf 0.05 \
–geno 0.1 \
–mind 0.5 \
–thin 0.4 \
-r2 gz \
–ld-window 100 \
–ld-window-kb 1000 \
–ld-window-r2 0 \
–make-bed \
–out tomato.chr3.ld
每个参数都是什么意思在上面的链接里都有介绍(这个计算起来非常快)
利用输出数据作图

R语言代码
library(data.table)
library(tidyverse)
dat.ld<-fread(“tomato.chr3.ld.ld.gz”) dat.ld%>%filter(abs(BP_A-BP_B)<=50000) -> dat.ld
save(dat.ld,file = “tomato.chr3.ld.Rdata”)

image.png
作图代码
library(tidyverse)
library(patchwork)
load(“tomato.chr3.ld.Rdata”)

dat.ld %>% head()

dat.ld.new<-dat.ld %>%
mutate(group=paste0(str_extract(SNP_A,pattern = “[A-Za-z]+$”),
“_”,
str_extract(SNP_B,pattern = “[A-Za-z]+$”)))

dat.ld.new %>%
mutate(new_group=case_when(
group == “SNP_SV”|group==”SV_SNP” ~ “SNP_SV”,
group == “INDEL_SV”|group==”SV_INDEL” ~ “INDEL_SV”,
group == “SNP_SNP” ~ “SNP_SNP”,
group == “SV_SV” ~ “SV_SV”,
TRUE ~ “OTHERS”
)) -> dat.ld.new

p1<-ggplot(data=dat.ld.new %>% filter(new_group==”SNP_SV”),
aes(x=R2))+
geom_histogram(aes(y=after_stat(count / sum(count))),
bins = 150,
alpha=0.8,
color=”grey”,
fill=”#009f73″)+
ylim(NA,0.025)+
theme_bw(base_size = 15)+
theme(panel.grid = element_blank())+
labs(title = “SNP_SV”)

p2<-ggplot(data=dat.ld.new %>% filter(new_group==”INDEL_SV”),
aes(x=R2))+
geom_histogram(aes(y=after_stat(count / sum(count))),
bins = 150,
alpha=0.8,
color=”grey”,
fill=”#56b4e8″)+
ylim(NA,0.025)+
theme_bw(base_size = 15)+
theme(panel.grid = element_blank())+
labs(title = “INDEL_SV”)

p3<-ggplot(data=dat.ld.new %>% filter(new_group==”SV_SV”),
aes(x=R2))+
geom_histogram(aes(y=after_stat(count / sum(count))),
bins = 150,
alpha=0.8,
color=”grey”,
fill=”#d55e00″)+
ylim(NA,0.025)+
theme_bw(base_size = 15)+
theme(panel.grid = element_blank())+
labs(title = “SV_SV”)

p4<-ggplot(data=dat.ld.new %>% filter(new_group==”SNP_SNP”),
aes(x=R2))+
geom_histogram(aes(y=after_stat(count / sum(count))),
bins = 150,
alpha=0.8,
color=”grey”,
fill=”#0072b1″)+
ylim(NA,0.025)+
theme_bw(base_size = 15)+
theme(panel.grid = element_blank())+
labs(title = “SNP_SNP”)

p1+p2+p3+p4+
plot_layout(ncol = 2,nrow = 2)

image.png
SNP 和SNP 的R2和论文中的图的分布还是挺像的，SNP和SV的分布还是不一样的，如果用上所有染色体的数据可能还会有变化