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% 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