require(dplyr)
require(stringr)
require(tidyverse)
require(ggplot2)
####Read data
snps <- read.delim("recmin/maf2_338_variants.simple.vcf")
str(snps)
in2rmal = c(18774475,14591034,15616195,8855601)
####Binary table/matrix
makebin = function(x) {ifelse(x==snps$A2.C7.Z1125,0,1)}
bin = select(snps,-MELREF, -CHROM,-num_A,-num_T,-num_G,-num_C, pos = POS)
head(bin)
bin = apply(bin[c(2:ncol(bin))],2,makebin)
bin = as.data.frame(bin)
bin$pos = snps$POS
head(bin)
pos = data.frame("pos1" = c(in2rmal[1],bin$pos))
pos$pos2 = c(bin$pos,in2rmal[2])
# ####Prepare and export for RecMin
# bin = as.matrix(bin)
# tbin = t(bin)
# colnames(tbin) = tbin[10,]
# tbin = tbin[-10,]
#
# write.table(tbin,'bin_data_refZI125.txt',col.names = TRUE,row.names = FALSE,sep = '')
####Visualize haplotype blocks
hdat = gather(bin, key = 'ind',value = 'nt',1:9)
head(hdat)
str(hdat)
hdat$ind = sub('.*\\.', '', hdat$ind)
hdat$ind = sub('Z1',"ZI",hdat$ind)
order_tree <- read_delim("recmin/tree.txt",
"\t", escape_double = FALSE, col_names = FALSE, trim_ws = TRUE)
colnames(order_tree) = c("ind","phylo")
order_tree = arrange(order_tree,desc(phylo))
hdat1 = left_join(hdat,order_tree,'ind') %>%
arrange(desc(phylo)) %>%
mutate(ind2 = factor(ind, levels=order_tree$ind))
hplot = ggplot(hdat1,aes(x = pos,y = ind2,color=nt,fill=nt ))+
geom_line(size=5,alpha=.8)+
xlim(c(min(in2rmal),max(in2rmal)))+
xlab("2R")+
ylab("SD")+
#geom_vline(xintercept = in2rmal,linetype='dotted')+
scale_color_gradient(low='darkorange',high = 'darkorange4')+
theme_minimal()+
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.line.x =element_blank(),
axis.text.x=element_blank(),
plot.margin = margin(t = 0,r = 0,b = 3,l = 0)
)
hplot
###Add marks where recombination happened
loc <- read.csv("recmin/recmin_output.txt",header = FALSE)
str(loc)
loc$V3 = (loc$V1+loc$V2)/2
rm = ggplot(loc,aes(x = V1))+
geom_point(y = 1,shape=2,size = 3) +
xlim(c(min(in2rmal),max(in2rmal)))+
ylim(c(0.9,1.1))+
ylab("Rm")+
xlab("2R")+
theme_minimal()+
#geom_vline(xintercept = in2rmal,linetype='dotted')+
theme(axis.line.y = element_blank(),
axis.text.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank())
rm
require(cowplot)
bin_hapl = plot_grid(hplot,rm,nrow=2,align = "v", rel_heights = c(8,2))
bin_hapl