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

2024-05-27 20:30:13 浏览数 (2)

论文

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号染色体过滤出来

代码语言：javascript复制
bcftools view 706.sv.vcf.gz -r 3 -O v -o chr3.sv.vcf

自己写一个python脚本修改一些vcf文件里的内容

把id 改成 chr pos "_SV”的形式，把INFO列的内容都去掉，把 alt 和 ref 都改成 单碱基的形式 基因型只保留前三个字符

代码语言：javascript复制
python 20240524_01.py chr3.sv.vcf chr3.sv.edited.vcf

20240524_01.py脚本的内容

代码语言：javascript复制
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("%sn"%('t'.join(new_list)))
            
fw.close()

这个vcf文件里不知道为啥会有很多 .:1 .:0这种基因型，如果是这种同意替换成./.

根据缺失率对数据进行过滤

代码语言：javascript复制
vcftools --vcf chr3.sv.edited.vcf --max-missing 0.8 --recode --recode-INFO-all --out chr3.sv.edited.filter

还剩12688个位点

做一个基因型填充

代码语言：javascript复制
beagle gt=chr3.sv.edited.filter.recode.vcf out=chr3.sv.edited.filter.impute nthreads=24

SNP INDEL数据集的处理

只保留二等位的位点，缺失率

代码语言：javascript复制
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

基因型填充

代码语言：javascript复制
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改一下

代码语言：javascript复制
bgzip -d chr3.snp.filter.impute.vcf.gz
python 20240524_02.py chr3.snp.filter.impute.vcf chr3.snp.filter.impute.edited.vcf

脚本内容

代码语言：javascript复制
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("%sn"%('t'.join(temp_list)))
            
fw.close()

把SNP INDEL和SV数据集合并

代码语言：javascript复制
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/

这里介绍的还挺详细的

代码语言：javascript复制
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语言代码

代码语言：javascript复制
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")

作图代码

代码语言：javascript复制
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)

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

数据分析 list 脚本 论文 数据

0 人点赞