library(tidyverse);library(data.table)

setwd("~/Desktop/ddRAD_epiRAD_powan/demultiplexed_reads/EpiRAD/Intermediate_files3/08.EpiRAD/read_depth")

my_data <- fread("reads_count_stacks2.txt", header = FALSE)
loci.name <- fread("loci.names.txt", header = FALSE)
id.name <- fread("id.names.txt", header = FALSE)
row.names(my_data) <- loci.name$V1
colnames(my_data) <- id.name$V1


mat <- as.matrix(my_data)
mat2 <- ifelse(mat < 2, 0, mat)

my_data2 <- data.frame(mat2)

lom <- my_data2[,c(1:25,74:116)]
row.names(lom) <- loci.name$V1
lom2 <- lom %>% rownames_to_column("Locus")
lom3 <- lom2[apply(lom, 1, function(x) length(x[x != 0]) >= 22.44),]

fwrite(lom, "../datasets/lom_counts.txt", sep = "\t", row.names = FALSE)
fwrite(lom3, "../datasets/lom_filtered_counts_rep_transformed.txt", sep = "\t", row.names = FALSE)


eck <- my_data[,c(26:73)]
row.names(eck) <- loci.name$V1
eck2 <- eck %>% rownames_to_column("Locus")
eck3 <- eck2[apply(eck, 1, function(x) length(x[x != 0]) >= 15.84),]

fwrite(eck, "../datasets/eck_counts.txt", sep = "\t", row.names = FALSE)
fwrite(eck3, "../datasets/eck_filtered_counts_rep.txt", sep = "\t", row.names = FALSE)


combined <- na.omit(left_join(eck3, lom3, by = "Locus"))
fwrite(combined, "../datasets/combined_counts.txt", sep = "\t", row.names = FALSE)


### create dataset to identify batch loci

batch_effect <- my_data[,c(10:58,79:111)]

row.names(batch_effect) <- loci.name$V1
batch_effect2 <- batch_effect %>% rownames_to_column("Locus")
batch_effect3 <- batch_effect2[apply(batch_effect, 1, function(x) length(x[x != 0]) >= 2),]
fwrite(batch_effect3, "../datasets/batch_data_counts.txt", sep = "\t", row.names = FALSE)



# eck
batch_effect_eck <- my_data[,c(26:58)]
row.names(batch_effect_eck) <- loci.name$V1
batch_effect_eck2 <- batch_effect_eck %>% rownames_to_column("Locus")
batch_effect_eck3 <- batch_effect_eck2[apply(batch_effect_eck, 1, function(x) length(x[x != 0]) >= 2),]
fwrite(batch_effect_eck3, "batch_data_eck.txt", sep = "\t", row.names = FALSE)

# lom
batch_effect_lom <- my_data[,c(10:25,79:111)]
row.names(batch_effect_lom) <- loci.name$V1
batch_effect_lom2 <- batch_effect_lom %>% rownames_to_column("Locus")
batch_effect_lom3 <- batch_effect_lom2[apply(batch_effect_lom, 1, function(x) length(x[x != 0]) >= 2),]
fwrite(batch_effect_lom3, "batch_data_lom.txt", sep = "\t", row.names = FALSE)



#### epiRAD replicates

replicates_eck <- eck3[,c(13,14,21,22,40,41)]
rep_eck_1 <- filter(replicates_eck, ELT08345X == 0 & ELT08345 > 0)
rep_eck_2 <- filter(replicates_eck, ELT08345X > 0 & ELT08345 == 0)

library_batch_loci <- read.table("../datasets/list_of_library_batch_loci_glmFit.txt", header = TRUE)
adapter_batch_loci <- read.table("../datasets/list_of_adapter_batch_loci_glmFit.txt", header = TRUE)

combined1 <- combined[!(combined$Locus %in% library_batch_loci$Locus),]
combined1 <- combined1[!(combined1$Locus %in% adapter_batch_loci$Locus),]