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

用户7010445

发布于 2024-05-27 20:30:13

1720

发布于 2024-05-27 20:30:13

文章被收录于专栏：小明的数据分析笔记本小明的数据分析笔记本

论文

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")

作图代码

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)