```r
#This script generates stacked barplots comparing TE content in the entire genome to that of piRNA clusters
#load libraries
library(dplyr)   
library(ggplot2)
library(tidyverse)
library(gtools)

my.outprefix <- paste(Sys.Date(),"Cluster",sep="_")

#read in the outfile from repeatmasker
out_file <- read.table("data/Amil/GCF_013753865.1_Amil_v2.1_genomic.fna.out", comment.char = "", header = F, fill = T, skip = 2)

#retains only 1)chromosome 2)start 3)stop 4)TE name 5)Family Name and 6)percent divergence
bed_file <- out_file[,c(5:7,11,2)]

#Assess how the families breakdown for the whole genome
#remove empty lines 
bed_file <- bed_file %>%
  drop_na(.)

#output bed file
write.table(bed_file, file = "output/18-Apul-piRNA/2_TE_piRNA_cluster_overlap/APUL_repeatmasker_bedfile.bed", sep = "\t", col.names = F, row.names = F, quote = F)

#retain only major family name and combine unclear and Unknown
bed_file$V11 <- factor(bed_file$V11)
bed_file$V11 <- gsub("/.*","",bed_file$V11)
bed_file$V11 <- gsub("unclear","Unknown",bed_file$V11)

#make a table of the families of TEs. How frequently do the different TE families appear in the genome?
table_of_TEs <- as.data.frame(table(bed_file$V11))

#remove simple repeats and unify "unknown" and "unclear"
table_of_TEs <- table_of_TEs[(!grepl("Low_complexity", table_of_TEs$Var1)),]
table_of_TEs <- table_of_TEs[(!grepl("Simple_repeat", table_of_TEs$Var1)),]

#create matrix with TE family frequency in the genome
cols.g <- rainbow(nrow(table_of_TEs))
table_of_TEs$percent = round(100*table_of_TEs$Freq/sum(table_of_TEs$Freq), digits = 1)
table_of_TEs$label = paste(table_of_TEs$Var1," (", table_of_TEs$percent,"%)", sep = "")
```