library(tidyverse)
library(cowplot)
### load depth data
temp = list.files(path = "metasembly_100bp/",pattern="*2R_meta-contig.txt",full.names = TRUE)
##do not load the recombinant SD-ZI
temp = temp[c(1:6,8:16)]
temp
myfiles = lapply(temp, function(x){read.delim(gzfile(x,'rt'),sep = '\t',header = FALSE)})
nms = gsub('metasembly_100bp//','',temp)
nms = gsub('.2R_meta-contig.txt', "",nms)
nms
names(myfiles) = nms
datos = do.call(rbind,myfiles)
head(datos)
colnames(datos) = c("contig","start",'stop','dp')
datos$line = sub("\\..*", "", rownames(datos))
head(datos)
datos$geno = ifelse(grepl('ERR',datos$line),'iso-1','sd-zi')
head(datos)
####cut to Rsp major locus
w = 180000
### average depth chrom 2 for normalization
ave.line <- read.table("metasembly_100bp/ave.depth.chrom2.txt", quote="\"", comment.char="")
colnames(ave.line) = c("line","ave.line")
datos = left_join(datos,ave.line) %>% mutate("norm.dp" = dp/ave.line)
# p = ggplot(filter(datos,stop<w),aes(x = start,y = norm.dp,color=geno))+
# geom_point()+
# facet_grid(.~contig,scales = 'free')+
# coord_cartesian(ylim=c(0,10))
#
# p
#### Load repeats location near Rsp major
repeat_bed = read.csv('meta-assembly_r6-PBcR-BLASR.V2.rsp_maj.txt',sep = '\t',header = FALSE,comment.char = "#")
head(repeat_bed)
repeat_bed = repeat_bed %>% filter(grepl(pattern = 'rsp',ignore.case = TRUE,x = V9)) %>%
select('rsp'=V9,'contig'=V1,"start"=V4,"stop"=V5) %>%
mutate('loc'='Rsp-major')
colnames(repeat_bed)=c('rsp','contig','start','stop','loc')
head(repeat_bed)
repeat_bed$start = as.numeric(repeat_bed$start)
repeat_bed$stop = as.numeric(repeat_bed$stop)
str(repeat_bed)
### define colors
cols = c("iso-1"="darkslategrey",'sd-zi'="darkorange")
# ggplot(datos,aes(x=start))+
# geom_rect(data=repeat_bed, aes(xmin = start, xmax = stop), ymin = -Inf, ymax = Inf, alpha = 0.1,fill="dodgerblue")+
# geom_point(aes(y = norm.dp,color=geno,fill=geno), alpha=.8,size=1)+
# geom_line(aes(y = norm.dp,color=geno,fill=geno))+
# scale_color_manual(values = cols)+
# scale_fill_manual(values=cols)+
# coord_cartesian(ylim=c(0,12))+
# theme_minimal()+
# ylab('Relative Depth')+
# #xlab("utg_2021")+
# theme(legend.position = 'none')
#### get average across all SD lines (exclude B6-Z1138 because low coverage in 2R)
sumdat = datos %>%
filter(line != 'B6-Z1138') %>%
group_by(contig,start,stop,geno) %>%
summarise(mean.norm.dp = mean(norm.dp),stdev=sd(norm.dp))
head(sumdat)
### plot depth in windows
dp = ggplot(filter(sumdat,stop<w),aes(x=start))+
geom_rect(data=repeat_bed, aes(xmin = start, xmax = stop), ymin = -Inf, ymax = Inf, alpha = 0.4,fill="gray")+
geom_ribbon(aes(ymax=mean.norm.dp+stdev,ymin=mean.norm.dp-stdev,fill=geno),alpha=.5)+
geom_point(aes(y = mean.norm.dp,color=geno,fill=geno), alpha=.8,size=1)+
#geom_line(aes(y = mean.norm.dp,color=geno,fill=geno))+
scale_color_manual(values = cols,labels=c('Iso-1','SD-Mal'))+
scale_fill_manual(values=cols,labels=c('Iso-1','SD-Mal'))+
#coord_cartesian(ylim=c(0,15))+
ylab('Relative Depth')+
#xlab("utg_2021")+
theme_minimal(base_size = 20)+
#theme(legend.position = c(0.8, 0.2))+
theme(legend.position = c(0.95,0.95),
legend.title = element_blank(),
legend.text = element_text(face='italic'))
dp
#################plot repeat coordinates from gff
source('meta_gff_plot.R')
require(cowplot)
plot_grid(dp,pgff,ncol = 1,align = "v")
dpm = dp +
theme(axis.title.x = element_blank(),
axis.text.x = element_blank())
dpm
pgff
pgff = pgff+
xlab(label = 'Rsp Major')+
theme(axis.title.y = element_blank(),
axis.title.x = element_text(face = 'italic'),
axis.text.y = element_text(face = 'italic'))
pgff
plot_grid(dpm,pgff,ncol = 1,align = "v",rel_heights = c(5,4))