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