rm()
gc()
##           used (Mb) gc trigger (Mb) max used (Mb)
## Ncells  522230 27.9    1181252 63.1   621665 33.3
## Vcells 1003983  7.7    8388608 64.0  1601566 12.3
if(!require(data.table)) install.packages("data.table")
## Loading required package: data.table
library(data.table)
if(!require(tidyverse)) install.packages("tidyverse")
## Loading required package: tidyverse
## -- Attaching packages ---------------------------------- tidyverse 1.2.1 --
## <U+221A> ggplot2 3.2.1     <U+221A> purrr   0.3.2
## <U+221A> tibble  2.1.3     <U+221A> dplyr   0.8.3
## <U+221A> tidyr   0.8.3     <U+221A> stringr 1.4.0
## <U+221A> readr   1.3.1     <U+221A> forcats 0.4.0
## -- Conflicts ------------------------------------- tidyverse_conflicts() --
## x dplyr::between()   masks data.table::between()
## x dplyr::filter()    masks stats::filter()
## x dplyr::first()     masks data.table::first()
## x dplyr::lag()       masks stats::lag()
## x dplyr::last()      masks data.table::last()
## x purrr::transpose() masks data.table::transpose()
library(tidyverse)
if(!require(grr)) install.packages("grr")
## Loading required package: grr
## 
## Attaching package: 'grr'
## The following object is masked from 'package:dplyr':
## 
##     matches
## The following object is masked from 'package:tidyr':
## 
##     extract
library(grr)

'%ni%' = Negate('%in%')

# getwd()

PW363

IDs, evg headers and seq lenghts

tr = read.table(file = "../output/intermediate/PW363/PW363.tr_lengths.tsv", 
               header = TRUE, 
               sep = "\t", 
               quote = NULL,
               dec = ".", 
               stringsAsFactors = FALSE,
               comment.char = "@")

cds = read.table(file = "../output/intermediate/PW363/PW363.cds_lengths.tsv", 
               header = TRUE, 
               sep = "\t", 
               quote = NULL,
               dec = ".", 
               stringsAsFactors = FALSE,
               comment.char = "@")

colnames(tr)[1] = "name"
colnames(cds)[1] = "name"

head(tr)
##                                                                                                  name
## 1   CLCdnPW10_1 evgclass=main,okay,match:TPW_comp43052_c0_seq1,pct:100/100/.; aalen=134,88%,partial5;
## 2            CLCdnPW10_56 evgclass=main,okay,match:VdnPW4_13959,pct:100/96/.; aalen=280,79%,complete;
## 3          CLCdnPW10_79 evgclass=main,okay,match:CLCdnPW27_111,pct:100/100/.; aalen=850,95%,complete;
## 4             CLCdnPW10_158 evgclass=main,okay,match:VdnPW2_559,pct:100/97/.; aalen=481,84%,complete;
## 5 CLCdnPW10_162 evgclass=main,okay,match:TPW_comp60167_c0_seq1,pct:100/100/.; aalen=200,86%,complete;
## 6           CLCdnPW10_203 evgclass=main,okay,match:VdnPW7_1580,pct:100/100/.; aalen=706,90%,complete;
##   seq qual length
## 1  NA   NA    456
## 2  NA   NA   1066
## 3  NA   NA   2668
## 4  NA   NA   1718
## 5  NA   NA    698
## 6  NA   NA   2337
head(cds)
##                                                                                                                                              name
## 1   CLCdnPW10_1 type=CDS; aalen=134,88%,partial5; clen=456; strand=-; offs=456-52;  evgclass=main,okay,match:TPW_comp43052_c0_seq1,pct:100/100/.;
## 2           CLCdnPW10_56 type=CDS; aalen=280,79%,complete; clen=1066; strand=+; offs=84-926;  evgclass=main,okay,match:VdnPW4_13959,pct:100/96/.;
## 3        CLCdnPW10_79 type=CDS; aalen=850,95%,complete; clen=2668; strand=+; offs=10-2562;  evgclass=main,okay,match:CLCdnPW27_111,pct:100/100/.;
## 4          CLCdnPW10_158 type=CDS; aalen=481,84%,complete; clen=1718; strand=+; offs=126-1571;  evgclass=main,okay,match:VdnPW2_559,pct:100/97/.;
## 5 CLCdnPW10_162 type=CDS; aalen=200,86%,complete; clen=698; strand=+; offs=61-663;  evgclass=main,okay,match:TPW_comp60167_c0_seq1,pct:100/100/.;
## 6        CLCdnPW10_203 type=CDS; aalen=706,90%,complete; clen=2337; strand=-; offs=2252-132;  evgclass=main,okay,match:VdnPW7_1580,pct:100/100/.;
##   seq qual length
## 1  NA   NA    405
## 2  NA   NA    843
## 3  NA   NA   2553
## 4  NA   NA   1446
## 5  NA   NA    603
## 6  NA   NA   2121
tr = tr[, -c(2,3)]
cds = cds[, -c(2,3)]

tmp = strsplit(tr$name, " ")
tmp = data.table::transpose(tmp)
length(tmp)
## [1] 3
ID = tmp[[1]]
evgHead1 = tmp[[2]]
evgHead2 = tmp[[3]]
## gsub("[;\"]"
evgHead2 = sapply(1:length(evgHead2), function(x) gsub("aalen=", "", evgHead2[x]))
evgHead2 = sapply(1:length(evgHead2), function(x) gsub(";", "", evgHead2[x]))
temp = strsplit(evgHead2, ",")
temp = data.table::transpose(temp)
evgAAlen = temp[[1]]
evgAAperc = temp[[2]]
evgAAcomplete = temp[[3]]

evgHead1 = sapply(1:length(evgHead1), function(x) gsub("evgclass=", "", evgHead1[x]))
evgHead1 = sapply(1:length(evgHead1), function(x) gsub(";", "", evgHead1[x]))
temp = strsplit(evgHead1, ",")
temp = data.table::transpose(temp)
class = temp[[1]]
pass = temp[[2]]
match = temp[[3]]
pct = temp[[4]]

tr.2 = cbind(ID, class, pass, match, pct, evgAAlen, evgAAperc, evgAAcomplete, tr$length)
colnames(tr.2)[dim(tr.2)[2]] = "length"

tmp = strsplit(cds$name, " ")
tmp = data.table::transpose(tmp)
length(tmp)
## [1] 8
ID = tmp[[1]]
evgHead1 = tmp[[8]]
evgHead2 = tmp[[3]]

evgHead2 = sapply(1:length(evgHead2), function(x) gsub("aalen=", "", evgHead2[x]))
evgHead2 = sapply(1:length(evgHead2), function(x) gsub(";", "", evgHead2[x]))
temp = strsplit(evgHead2, ",")
temp = data.table::transpose(temp)
evgAAlen = temp[[1]]
evgAAperc = temp[[2]]
evgAAcomplete = temp[[3]]

evgHead1 = sapply(1:length(evgHead1), function(x) gsub("evgclass=", "", evgHead1[x]))
evgHead1 = sapply(1:length(evgHead1), function(x) gsub(";", "", evgHead1[x]))
temp = strsplit(evgHead1, ",")
temp = data.table::transpose(temp)
class = temp[[1]]
pass = temp[[2]]
match = temp[[3]]
pct = temp[[4]]

clen = gsub("clen=", "", tmp[[4]])
strand = gsub("strand=", "", tmp[[5]])
offs = gsub("offs=", "", tmp[[6]])
clen = gsub(";", "", clen)
strand = gsub(";", "", strand)
offs = gsub(";", "", offs)

# dim(cbind(clen, strand, offs))

cds.2 = cbind(ID, class, pass, match, pct, evgAAlen, evgAAperc, evgAAcomplete, clen, strand, offs, cds$length)
colnames(cds.2)[dim(cds.2)[2]] = "length"

utrorfs

dim(cds.2) - dim(tr.2)
## [1] 3066    3
utrorf = cds.2[grep("utrorf", cds.2[,1])]
tr.utrorf = gsub("utrorf", "", utrorf)
ind = match(tr.utrorf, tr.2[,1])

cat("missing in .tr: ", sum(is.na(ind)), "sequences while existing in .cds as utrorf \n")
## missing in .tr:  2608 sequences while existing in .cds as utrorf
ind2 = which(is.na(ind))
ind3 = match(utrorf[ind2],cds.2[,1])

mymat = matrix(NA, length(ind2), dim(tr.2)[2])
mymat[,1] = tr.utrorf[ind2]
colnames(mymat) = colnames(tr.2)
mymat = as.data.frame(mymat)
mymat$class = rep("lost.dropped", length(ind2))
mymat$pass = rep("as.utrorfs.cds", length(ind2))

# tr.2 = rbind(tr.2, mymat) # shift to utrorf file


uniqueness = gsub("utrorf", "", cds.2[,1])
n_occur <- data.frame(table(uniqueness))
duplicates = n_occur[n_occur$Freq > 1,]
# dim(duplicates)
## partial match
## duplicatedCDS = sapply(1:dim(duplicates)[1], function(x) grep(as.character(duplicates[x,1]), cds.2[,1]))
## duplicatedCDSTable = cds.2[unlist(duplicatedCDS),]
duplicatedCDSTable = rbind (cds.2[match(duplicates[,1],cds.2[,1]),], 
                            cds.2[match(paste0(duplicates[,1], "utrorf"),cds.2[,1]),])
duplicatedCDSTable = as.data.frame(duplicatedCDSTable)
duplicatedCDSTable = duplicatedCDSTable[with(duplicatedCDSTable, order(ID)), ]


cat(".tr with 2nd .cds: ", dim(duplicatedCDSTable)[1]/2, " \n")
## .tr with 2nd .cds:  458
write.table(duplicatedCDSTable, file = "../output/PW363/PW363_secondaryORFTable.tsv", 
            append = FALSE, quote = FALSE, sep = "\t",
            eol = "\n", na = "NA", dec = ".", row.names = FALSE,
            col.names = TRUE, qmethod = c("escape", "double"),
            fileEncoding = "")
colnames(tr.2)[2:dim(tr.2)[2]] = paste0("tr.", colnames(tr.2)[2:dim(tr.2)[2]])
colnames(cds.2)[2:dim(cds.2)[2]] = paste0("cds.", colnames(cds.2)[2:dim(cds.2)[2]])

# http://www.datasciencemadesimple.com/join-in-r-merge-in-r/

tr.3 = merge(tr.2, cds.2, by="ID", all.x =TRUE)


utrorf = cds.2[grep("utrorf", cds.2[,1])]
tr.utrorf = gsub("utrorf", "", utrorf)


ind4 = match(tr.utrorf, tr.2[,1])
ind5 = match(utrorf, cds.2[,1])
# colnames(cds.2[ind5,])  
# colnames(tr.2[ind4,])
tmp = (cds.2[ind5,])  
# dim(tr.2[ind4,]) 
colnames(tmp) = colnames(cds.2)
colnames(tmp)[1] = "ID2"
tmp = cbind(tmp, gsub("utrorf", "", tmp[,1]))
colnames(tmp)[dim(tmp)[2]] = "ID"

tr.3_utrorf = merge(tr.2, tmp, by="ID", all.y =TRUE)
colnames(tr.3_utrorf)[1] = "ID3"
colnames(tr.3_utrorf)[dim(tr.2)[2]+1] = "ID"
colnames(tr.3_utrorf)[1] = "ID2"
# dim(tr.3_utrorf)

colnames(tr.3_utrorf)
##  [1] "ID2"               "tr.class"          "tr.pass"          
##  [4] "tr.match"          "tr.pct"            "tr.evgAAlen"      
##  [7] "tr.evgAAperc"      "tr.evgAAcomplete"  "tr.length"        
## [10] "ID"                "cds.class"         "cds.pass"         
## [13] "cds.match"         "cds.pct"           "cds.evgAAlen"     
## [16] "cds.evgAAperc"     "cds.evgAAcomplete" "cds.clen"         
## [19] "cds.strand"        "cds.offs"          "cds.length"
colnames(mymat)
## [1] "ID"            "class"         "pass"          "match"        
## [5] "pct"           "evgAAlen"      "evgAAperc"     "evgAAcomplete"
## [9] "length"
mymat2 = matrix(NA, dim(mymat)[1], dim(tr.3_utrorf)[2])
colnames(mymat2) = colnames(tr.3_utrorf)
mymat2 = as.data.frame(mymat2)
mymat2$ID2 = mymat$ID
mymat2$tr.class = mymat$class
mymat2$tr.pass = mymat$pass
mymat2$tr.length = mymat$length
mymat2$ID = mymat$ID
i <- sapply(tr.3_utrorf, is.factor) 
tr.3_utrorf[i] <- lapply(tr.3_utrorf[i], as.character) 
i <- sapply(mymat2, is.factor) 
mymat2[i] <- lapply(mymat2[i], as.character) 
dim(mymat2)
## [1] 2608   21
dim(tr.3_utrorf)
## [1] 3066   21
sum(mymat2$ID2 %in% tr.3_utrorf$ID2)
## [1] 2608
ind = match(mymat2$ID2, tr.3_utrorf$ID2)
# names(tr.3_utrorf)
tr.3_utrorf$tr.class[ind] = mymat2$tr.class # shift to utrorf file; already in
tr.3_utrorf$tr.pass[ind] = mymat2$tr.pass # shift to utrorf file; already in

sum(duplicatedCDSTable$ID %in% tr.3_utrorf$ID2)
## [1] 458
sum(duplicatedCDSTable$ID %in% tr.3_utrorf$ID)
## [1] 458
tmp = tr.3[which(tr.3$ID %in% duplicatedCDSTable$ID ),]
tmp$ID2 = tmp$ID
ind = match(colnames(tr.3_utrorf), colnames(tmp))
tmp = tmp[,ind]
removeAdd = match(tmp$ID, tr.3$ID)

tr.3_utrorf = rbind(tmp, tr.3_utrorf) # shift to utrorf file

tr.3 = tr.3[-removeAdd,] # shift to utrorf file

lenghts

hist(as.numeric(tr.3$tr.length), 
     breaks = seq(0, max(as.numeric(tr.3$tr.length), na.rm = TRUE) + 50, 50),
     main = "PW363 evigene transcripts length", xlab = "transcript length",
     xlim=c(0, max(300,max(as.numeric(tr.3$tr.length), na.rm = TRUE))),
     col = "grey90")

hist(as.numeric(tr.3$cds.length), 
     breaks = seq(0, max(as.numeric(tr.3$cds.length), na.rm = TRUE) + 50, 50),
     main = "PW363 evigene coding sequences length", xlab = "CDS length",
     xlim=c(0, max(300,max(as.numeric(tr.3$cds.length), na.rm = TRUE))),
     col = "grey45", las=2)
abline(v=180, col = "red")
abline(v=200, col = "blue")
abline(v=300, col = "green")
axis(1, at=c(180), labels=c(180), las=2, col = c("red"), cex.axis=0.75, col.axis = c("red"))
axis(1, at=c(200), labels=c(200), las=2, col = c("blue"), cex.axis=0.75, col.axis = c("blue"))
axis(1, at=c(300), labels=c(300), las=2, col = c("green"), cex.axis=0.75, col.axis = c("green"))

hist(as.numeric(tr.3_utrorf$tr.length), 
     breaks = seq(0, max(as.numeric(tr.3_utrorf$tr.length), na.rm = TRUE) + 50, 50),
      main = "PW363 evigene transcripts with secondary ORF length", xlab = "transcript length",
     xlim=c(0, max(300,max(as.numeric(tr.3_utrorf$tr.length), na.rm = TRUE))),
     col = "grey90")

hist(as.numeric(tr.3_utrorf$cds.length), 
     breaks = seq(0, max(as.numeric(tr.3_utrorf$cds.length), na.rm = TRUE) + 50, 50),
     main = "PW363 evigene secondary ORF coding sequences length", xlab = "CDS length",
     xlim=c(0, max(300,max(as.numeric(tr.3_utrorf$cds.length), na.rm = TRUE))),
     col = "grey45", las=2)
abline(v=180, col = "red")
abline(v=200, col = "blue")
abline(v=300, col = "green")
axis(1, at=c(180), labels=c(180), las=2, col = c("red"), cex.axis=0.75, col.axis = c("red"))
axis(1, at=c(200), labels=c(200), las=2, col = c("blue"), cex.axis=0.75, col.axis = c("blue"))
axis(1, at=c(300), labels=c(300), las=2, col = c("green"), cex.axis=0.75, col.axis = c("green"))

Match Annot

MA = read.table(file = "../output/intermediate/PW363/PW363.tr_DM_STARlong.Aligned.out.sorted.sam.matchAnnot.parsed.txt", 
               header = FALSE, 
               sep = " ", 
               quote = NULL,
               dec = ".", 
               stringsAsFactors = FALSE,
               na.strings = "",
               comment.char = "@")

# dim(MA)
head(MA)
##                V1             V2                 V3 V4 V5
## 1 CLCdnPW36_76306 no_genes_found               <NA> NA NA
## 2 CLCdnPW12_58062 no_genes_found               <NA> NA NA
## 3 CLCdnPW12_58062 no_genes_found               <NA> NA NA
## 4  CLCdnPW12_1591 no_genes_found               <NA> NA NA
## 5  CLCdnPW12_1591 no_genes_found               <NA> NA NA
## 6   SdnPW1_197082 Sotub03g026420 Sotub03g026420.1.1  1  0
colnames(MA) = c("ID", "DMgeneID", "DMtrID", "exon_match", "match_score")
MA = MA[order(MA$ID),] 
duplRow = which(duplicated(MA) | duplicated(MA[nrow(MA):1, ])[nrow(MA):1])
MA[duplRow[1:10],]
##                   ID             DMgeneID               DMtrID exon_match
## 14321 CLCdnPW1_10007       Sotub02g012770   Sotub02g012770.1.1          1
## 14322 CLCdnPW1_10007       Sotub02g012770   Sotub02g012770.1.1          1
## 14881 CLCdnPW1_10013       Sotub02g014920   Sotub02g014920.1.1          1
## 14882 CLCdnPW1_10013       Sotub02g014920   Sotub02g014920.1.1          1
## 85784 CLCdnPW1_10051 PGSC0003DMG400015125 PGSC0003DMT400039120          1
## 85785 CLCdnPW1_10051 PGSC0003DMG400015125 PGSC0003DMT400039120          1
## 91683 CLCdnPW1_10061       no_genes_found                 <NA>         NA
## 91684 CLCdnPW1_10061       no_genes_found                 <NA>         NA
## 39373 CLCdnPW1_10107       no_genes_found                 <NA>         NA
## 39374 CLCdnPW1_10107       no_genes_found                 <NA>         NA
##       match_score
## 14321           0
## 14322           0
## 14881           5
## 14882           5
## 85784           5
## 85785           5
## 91683          NA
## 91684          NA
## 39373          NA
## 39374          NA
MA = MA[!duplicated(MA),]

hist(MA$match_score, main = 'All Match Annot scores', xlab = "MA score")

data <- data.frame(cbind(seq(1,dim(MA)[1],1),MA))
colnames(data)[1] = "enumerate"
setDT(data)
DT5 = data[ , .(match.score = paste(match_score,
                         collapse=",")), by = ID]

no.gene_score.0 = 
    unlist(sapply(1:length(DT5$match.score), function(i) {
      aa = as.numeric(unlist(strsplit(DT5$match.score[i], ",")))
      if(all(is.na(aa)) | max(aa, na.rm=TRUE) == 0) return(i)
    }))
ind = unlist(matches((unlist(DT5[no.gene_score.0 ,1])), MA[,1],  list =TRUE))
no.gene_score.0 = MA[ind ,]
unique(no.gene_score.0[,5])
## [1]  0 NA
MA.2 = MA[-which(MA$match_score == 0),]
MA = MA.2[-which(is.na(MA.2$match_score)),]

# only ones with max score
# uniqueID = unique(MA$ID)
keep = NULL
# to slow!!!
# keep = unlist(
#   sapply(1:length(uniqueID), function(i) {
#     if (!(i %% 100)) print(i)
#       # print(uniqueID[i])
#       ind = matches(uniqueID[i], MA[,1], all.x=FALSE,all.y=FALSE)[,2]
#       mmax = max(MA[ind,5])
#       ind2 = matches(mmax, MA[ind,5], all.x=FALSE,all.y=FALSE)[,2]
#       c(keep,(ind[ind2]))
#     }
#   )
# )

data <- data.frame(cbind(seq(1,dim(MA)[1],1),MA))
colnames(data)[1] = "enumerate"

dd = data %>% group_by(ID) %>% summarize(maxID = paste(enumerate[which(match_score == max(match_score))],  
                                                       collapse=","))



keep = as.numeric(unlist(strsplit(dd$maxID, ",")))
dim(MA)[1] - length(keep)
## [1] 192
MA = MA[keep,]
hist(MA$match_score, main = 'Selected Match Annot scores', xlab = "MA score")

setDT(MA)
setDT(no.gene_score.0)
# dim(MA)
DT1 = MA[ , .(DM.geneID = paste(DMgeneID,
                         collapse="; ")), by = ID]
DT2 = MA[ , .(DM.trID = paste(DMtrID,
                         collapse="; ")), by = ID]
DT3 = MA[ , .(exon.match = paste(exon_match,
                         collapse="; ")), by = ID]
DT4 = MA[ , .(match.score = paste(match_score,
                         collapse="; ")), by = ID]
DT1.1 = no.gene_score.0[ , .(DM.geneID = paste(DMgeneID,
                         collapse="; ")), by = ID]
DT2.1 = no.gene_score.0[ , .(DM.trID = paste(DMtrID,
                         collapse="; ")), by = ID]
DT3.1 = no.gene_score.0[ , .(exon.match = paste(exon_match,
                         collapse="; ")), by = ID]
DT4.1 = no.gene_score.0[ , .(match.score = paste(match_score,
                         collapse="; ")), by = ID]

# dim(MA)

myVec = unique(MA[,1])
myVec.1 = unique(no.gene_score.0[,1])
# dim(myVec)
# dim(myVec.1)
# dim(DT1)
# dim(DT1.1)

myVec = merge(myVec, DT1, by="ID", all.x =TRUE)
myVec = merge(myVec, DT2, by="ID", all.x =TRUE)
myVec = merge(myVec, DT3, by="ID", all.x =TRUE)
myVec = merge(myVec, DT4, by="ID", all.x =TRUE)
myVec.1 = merge(myVec.1, DT1.1, by="ID", all.x =TRUE)
myVec.1 = merge(myVec.1, DT2.1, by="ID", all.x =TRUE)
myVec.1 = merge(myVec.1, DT3.1, by="ID", all.x =TRUE)
myVec.1 = merge(myVec.1, DT4.1, by="ID", all.x =TRUE)

# dim(myVec)
# dim(myVec.1)

no.gene_score.0 = merge(tr.3, myVec.1, all.y =TRUE)

# print(dim(merge(tr.3, myVec, by="ID", all.x =TRUE)))

tr.3 = merge(tr.3, myVec, by="ID", all.x =TRUE)
colnames(myVec)[1] = "ID2"
tr.3_utrorf = merge(tr.3_utrorf, myVec, by="ID2", all.x =TRUE)
# dim(tr.3_utrorf)

IPR

IPR = read.table(file = "../output/intermediate/PW363/PW363_IPS_filtered_aggregated_filtered.tsv", 
               header = TRUE, 
               sep = "\t", 
               quote = NULL,
               dec = ".", 
               stringsAsFactors = FALSE,
               comment.char = "@")

# dim(IPR)
colnames(IPR)[1] = c("ID")

# print(dim(merge(tr.3, IPR, by="ID", all.x =TRUE)))

tr.3 = merge(tr.3, IPR, by="ID", all.x = TRUE)
tr.3_utrorf = merge(tr.3_utrorf, IPR, by="ID", all.x = TRUE)
# dim(tr.3_utrorf)
no.gene_score.0 = merge(no.gene_score.0, IPR, all.x = TRUE)

colnames(tr.3)[which(colnames(tr.3) == "Analysis")] = "IPS_Analysis"
colnames(tr.3_utrorf)[which(colnames(tr.3_utrorf) == "Analysis")] = "IPS_Analysis"
colnames(no.gene_score.0)[which(colnames(no.gene_score.0) == "Analysis")] = "IPS_Analysis"

colnames(tr.3)[which(colnames(tr.3) == "Signature_Accession")] = "IPS_Signature_Accession"
colnames(tr.3_utrorf)[which(colnames(tr.3_utrorf) == "Signature_Accession")] = "IPS_Signature_Accession"
colnames(no.gene_score.0)[which(colnames(no.gene_score.0) == "Signature_Accession")] = "IPS_Signature_Accession"

vecscreen

vs = read.table(file = "../output/intermediate/PW363/PW363_vecscreen.tsv", 
               header = TRUE, 
               sep = "\t", 
               quote = NULL,
               dec = ".", 
               stringsAsFactors = FALSE,
               comment.char = "@")

# dim(vs)
colnames(vs)[1] = "ID"
length(unique(vs$ID))
## [1] 7900
setDT(vs)
# dim(vs)
DT1 = vs[ , .(coverage = paste(paste0(Matching_vector._starting_with_uv., " [",
                                       Lower_end_of_the_alignment_in_the_vector, "-",
                                       Upper_end_of_the_alignment_in_the_vector, "] ",
                                      coverage, "%"),
                         collapse="; ")), by = ID]
DT2 = vs[ , .(vectorEvidence = paste(paste0(The_strength_of_this_vecscreen_match_, ", ",
                                       The_strength_of_the_strongest_vecscreen_match_for_this_query, ", ",
                                       Whether_there_is_any_dangling_part_.called_.Suspect._by_vecscreen._at_either_end_of_the_query, ", ",
                                       the_classification_of_the_match),
                         collapse="; ")), by = ID]
DT3 = vs[ , .(blastnVector = paste(blastn_desc, collapse="; ")), by = ID]


myVec = unique(vs[,1])
# dim(myVec)
# dim(DT1)

myVec = merge(myVec, DT1, by="ID", all.x =TRUE)
myVec = merge(myVec, DT2, by="ID", all.x =TRUE)
myVec = merge(myVec, DT3, by="ID", all.x =TRUE)



# print(dim(merge(tr.3, myVec, by="ID", all.x =TRUE)))

tr.3 = merge(tr.3, myVec, by="ID", all.x =TRUE)
no.gene_score.0 = merge(no.gene_score.0, myVec, all.x = TRUE)
colnames(myVec)[1] = "ID2"
tr.3_utrorf = merge(tr.3_utrorf, myVec, by="ID2", all.x =TRUE)
# dim(tr.3_utrorf)

colnames(tr.3)[which(colnames(tr.3) == "coverage")] = "VecScreen_coverage"
colnames(tr.3_utrorf)[which(colnames(tr.3_utrorf) == "coverage")] = "VecScreen_coverage"
colnames(no.gene_score.0)[which(colnames(no.gene_score.0) == "coverage")] = "VecScreen_coverage"

DIAMOND BLASTX

(myfiles = list.files(path = "../output/intermediate/PW363/", pattern = "ENCHformat_top1.tsv"))
## [1] "01_PW363.cds_solanumTuberosum.out.ENCHformat_top1.tsv"
## [2] "02_PW363.tr_solanumTuberosum.out.ENCHformat_top1.tsv" 
## [3] "03_PW363.cds_Solanaceae.out.ENCHformat_top1.tsv"      
## [4] "04_PW363.tr_Solanaceae.out.ENCHformat_top1.tsv"       
## [5] "05_PW363.cds_SP_TrEMBL_plants.out.ENCHformat_top1.tsv"
## [6] "06_PW363.tr_SP_TrEMBL_plants.out.ENCHformat_top1.tsv" 
## [7] "07_PW363.cds_SwissProt_TrEMBL.out.ENCHformat_top1.tsv"
## [8] "08_PW363.tr_SwissProt_TrEMBL.out.ENCHformat_top1.tsv"
# watch our for comment.chars # in tables
for (i in myfiles) {
  print(i)
  blast = read.table(file = paste0("../output/intermediate/PW363/",i), 
               header = TRUE, 
               sep = "\t", 
               quote = NULL,
               dec = ".", 
               stringsAsFactors = FALSE,
               comment.char = "@")
  
  blast = blast[,-1]
  # if (grepl(".tr_", i)) {
  #   remove = which(!grepl("evgclass", blast[,1])) # lost
  #   blast = blast[-remove,]
  # }
  tmp = strsplit(blast[,1], " ")
  tmp = data.table::transpose(tmp)[[1]]
  blast[,1] = tmp
  remove = which(duplicated(blast[,1])) # duplicated
  if (length(remove)) blast = blast[-remove,]
  
  colnames(blast)[1] = "ID"
  blast$Gaps.num = round(blast$Gaps.*blast$Aligned_seq_length/100)
  blast$Query_coverage = (blast$Aligned_seq_length-blast$Gaps.num)/(blast$Query_length/3)
  blast$Target_coverage = (blast$Aligned_seq_length-blast$Gaps.num)/blast$Target_length
  blast$Query_Target_ratio = (blast$Query_length/3)/blast$Target_length
  ind1 = which(blast$Query_coverage >= 0.50)
  ind2 = which(blast$Target_coverage >= 0.50)
  ind = intersect(ind1, ind2)
  blast = blast[ind, ]
  blast$Target_coverage[blast$Target_coverage > 1.0] = 1.0
  blast$Query_coverage[blast$Query_coverage > 1.0] = 1.0
  ind = which(colnames(blast) %in% c("ID", "Target_ID", "Target_description", "Aligned_seq_length", "Target_coverage", "Query_coverage", 
                         "E.value", "Score"))
  blast = blast[, ind]

  
  j = gsub(".out.ENCHformat_top1.tsv", "", i)
  j = gsub("_PW363", "", j)
  colnames(blast)[2:dim(blast)[2]] = paste0(j, "_", colnames(blast)[2:dim(blast)[2]])
  
  tr.3 = merge(tr.3, blast, by="ID", all.x =TRUE)
  no.gene_score.0 = merge(no.gene_score.0, blast, by="ID", all.x = TRUE)

  if ((grepl(".cds_", i))) {
    tr.3_utrorf = merge(tr.3_utrorf, blast, by="ID", all.x =TRUE)
    # cat (".cds:", dim(merge(tr.3_utrorf, blast, by="ID", all.x =TRUE)), "\n")
    # print(dim(tr.3_utrorf))
  } else {
    colnames(blast)[1] = "ID2"
    tr.3_utrorf = merge(tr.3_utrorf, blast, by="ID2", all.x =TRUE)
    # cat(".tr: ", dim(merge(tr.3_utrorf, blast, by="ID2", all.x =TRUE)), "\n")
    # print(dim(tr.3_utrorf))
  }
    
}
## [1] "01_PW363.cds_solanumTuberosum.out.ENCHformat_top1.tsv"
## [1] "02_PW363.tr_solanumTuberosum.out.ENCHformat_top1.tsv"
## [1] "03_PW363.cds_Solanaceae.out.ENCHformat_top1.tsv"
## [1] "04_PW363.tr_Solanaceae.out.ENCHformat_top1.tsv"
## [1] "05_PW363.cds_SP_TrEMBL_plants.out.ENCHformat_top1.tsv"
## [1] "06_PW363.tr_SP_TrEMBL_plants.out.ENCHformat_top1.tsv"
## [1] "07_PW363.cds_SwissProt_TrEMBL.out.ENCHformat_top1.tsv"
## [1] "08_PW363.tr_SwissProt_TrEMBL.out.ENCHformat_top1.tsv"
### save.image("../output/PW363/PW363_withFiltering.RData")

save(tr.3, tr.3_utrorf, no.gene_score.0, file = "../output/PW363/tr3s.RData")
rm(list=ls())
gc()
##           used (Mb) gc trigger  (Mb)  max used  (Mb)
## Ncells  939040 50.2    4351378 232.4   5439222 290.5
## Vcells 3729346 28.5   91994926 701.9 114993657 877.4
load("../output/PW363/tr3s.RData")

library(data.table)

library(tidyverse)

library(grr)

write

i <- sapply(tr.3, is.factor) 
tr.3[i] <- lapply(tr.3[i], as.character) 
tr.3[is.na(tr.3)] = NA
tr.3[tr.3 == ""] = NA

i <- sapply(no.gene_score.0, is.factor) 
no.gene_score.0[i] <- lapply(no.gene_score.0[i], as.character) 
no.gene_score.0[is.na(no.gene_score.0)] = NA
no.gene_score.0[no.gene_score.0 == ""] = NA
# write.table(no.gene_score.0, "../output/PW363/no.gene_score.0_all.tsv", 
#             sep = "\t", row.names = FALSE,
#             quote = FALSE)

hist(as.numeric(no.gene_score.0$cds.length), 
     breaks = seq(0, max(as.numeric(no.gene_score.0$cds.length), na.rm = TRUE) + 50, 50),
     main = "PW363 evigene MA score NA/0 CDS sequence length", xlab = "CDS length",
     xlim=c(0, max(300,max(as.numeric(no.gene_score.0$cds.length), na.rm = TRUE))),
     col = "grey45", las=2)
abline(v=180, col = "red")
abline(v=200, col = "blue")
abline(v=300, col = "green")
axis(1, at=c(180), labels=c(180), las=2, col = c("red"), cex.axis=0.75, col.axis = c("red"))
axis(1, at=c(200), labels=c(200), las=2, col = c("blue"), cex.axis=0.75, col.axis = c("blue"))
axis(1, at=c(300), labels=c(300), las=2, col = c("green"), cex.axis=0.75, col.axis = c("green"))

table(as.numeric(no.gene_score.0$cds.length) >= 300)
## 
## FALSE  TRUE 
## 53928 15204
table(as.numeric(no.gene_score.0$cds.length) >= 200)
## 
## FALSE  TRUE 
## 46261 22871
table(as.numeric(no.gene_score.0$cds.length) >= 180)
## 
## FALSE  TRUE 
## 41976 27156
i <- sapply(tr.3_utrorf, is.factor) 
tr.3_utrorf[i] <- lapply(tr.3_utrorf[i], as.character) 
tr.3_utrorf[is.na(tr.3_utrorf)] = NA
tr.3_utrorf[tr.3_utrorf == ""] = NA

keep or remove

# colnames(tr.3_utrorf)

okay = which(!is.na(tr.3_utrorf$DM.trID))
length(okay)
## [1] 108
print("check if motif is in contamination")
## [1] "check if motif is in contamination"
okay = c(okay, which(!is.na(tr.3_utrorf$IPS_Analysis)))
length(okay)
## [1] 2326
ind = grep("Target_ID", colnames(tr.3_utrorf))

okay = c(okay, 
         which(!is.na(tr.3_utrorf[,ind[1]])),
         which(!is.na(tr.3_utrorf[,ind[2]])),
         which(!is.na(tr.3_utrorf[,ind[3]])),
         which(!is.na(tr.3_utrorf[,ind[4]])),
         which(!is.na(tr.3_utrorf[,ind[5]])),
         which(!is.na(tr.3_utrorf[,ind[6]])))

# which(as.numeric(tr.3_utrorf$cds.length) >= 180)
# which(as.numeric(tr.3_utrorf$cds.length) < 180)

okay = sort(unique(okay))
length(okay)
## [1] 2820
tr.3_utrorf.annot = tr.3_utrorf[okay,]
tr.3_utrorf.unknown = tr.3_utrorf[-okay,]

hist(as.numeric(tr.3_utrorf.annot$tr.length), 
     breaks = seq(0, max(as.numeric(tr.3_utrorf.annot$tr.length), na.rm = TRUE) + 50, 50),
      main = "OK PW363 evigene transcript with secondary ORF length", xlab = "transcript length",
     xlim=c(0, max(300,max(as.numeric(tr.3_utrorf.annot$tr.length), na.rm = TRUE))))

hist(as.numeric(tr.3_utrorf.annot$cds.length), 
     breaks = seq(0, max(as.numeric(tr.3_utrorf.annot$cds.length), na.rm = TRUE) + 50, 50),
     main = "OK PW363 evigene secondary ORF coding sequence length", xlab = "CDS length",
     xlim=c(0, max(300,max(as.numeric(tr.3_utrorf.annot$cds.length), na.rm = TRUE))),
     col = "grey45", las=2)
abline(v=180, col = "red")
abline(v=200, col = "blue")
abline(v=300, col = "green")
axis(1, at=c(180), labels=c(180), las=2, col = c("red"), cex.axis=0.75, col.axis = c("red"))
axis(1, at=c(200), labels=c(200), las=2, col = c("blue"), cex.axis=0.75, col.axis = c("blue"))
axis(1, at=c(300), labels=c(300), las=2, col = c("green"), cex.axis=0.75, col.axis = c("green"))

table(as.numeric(tr.3_utrorf.annot$cds.length) >= 300)
## 
## FALSE  TRUE 
##     2  2818
table(as.numeric(tr.3_utrorf.annot$cds.length) >= 200)
## 
## TRUE 
## 2820
table(as.numeric(tr.3_utrorf.annot$cds.length) >= 180)
## 
## TRUE 
## 2820
hist(as.numeric(tr.3_utrorf.unknown$tr.length), 
     breaks = seq(0, max(as.numeric(tr.3_utrorf.unknown$tr.length), na.rm = TRUE) + 50, 50),
      main = "notOK PW363 evigene transcript with secondary ORF length", xlab = "transcript length",
     xlim=c(0, max(300,max(as.numeric(tr.3_utrorf.unknown$tr.length), na.rm = TRUE))))

hist(as.numeric(tr.3_utrorf.unknown$cds.length), 
     breaks = seq(0, max(as.numeric(tr.3_utrorf.unknown$cds.length), na.rm = TRUE) + 50, 50),
     main = "notOK PW363 evigene secondary ORF coding sequence length", xlab = "CDS length",
     xlim=c(0, max(300,max(as.numeric(tr.3_utrorf.unknown$cds.length), na.rm = TRUE))),
     col = "grey45", las=2)
abline(v=180, col = "red")
abline(v=200, col = "blue")
abline(v=300, col = "green")
axis(1, at=c(180), labels=c(180), las=2, col = c("red"), cex.axis=0.75, col.axis = c("red"))
axis(1, at=c(200), labels=c(200), las=2, col = c("blue"), cex.axis=0.75, col.axis = c("blue"))
axis(1, at=c(300), labels=c(300), las=2, col = c("green"), cex.axis=0.75, col.axis = c("green"))

table(as.numeric(tr.3_utrorf.unknown$cds.length) >= 300)
## 
## TRUE 
##  704
table(as.numeric(tr.3_utrorf.unknown$cds.length) >= 200)
## 
## TRUE 
##  704
table(as.numeric(tr.3_utrorf.unknown$cds.length) >= 180)
## 
## TRUE 
##  704
# shorty = tr.3_utrorf.annot[which(as.numeric(tr.3_utrorf.annot$cds.length) < 300),]

# weight =  㤼㸶2.1 + 0.02 搼㸷 Similarity% (or Positives%)) + 0.01 搼㸷 Coverage% > 0
# S ... Similarity% (BLAST-tsv file: Positives%) 
# C ... Coverage% = Aligned_seq_length / minimum(Query_length/1 OR 3 for prot, Target_length) 


tr.3_utrorf.annot[is.na(tr.3_utrorf.annot)] = "-"
tr.3_utrorf.unknown[is.na(tr.3_utrorf.unknown)] = "-"

# write.table(tr.3_utrorf.annot, file = "../output/PW363/PW363_tr.3_utrorf_annot_all.tsv", 
#             append = FALSE, quote = FALSE, sep = "\t",
#             eol = "\n", na = "NA", dec = ".", row.names = FALSE,
#             col.names = TRUE, qmethod = c("escape", "double"),
#             fileEncoding = "")

# write.table(tr.3_utrorf.unknown, file = "../output/PW363/remove/PW363_tr.3_utrorf_unknown.tsv", 
#             append = FALSE, quote = FALSE, sep = "\t",
#             eol = "\n", na = "NA", dec = ".", row.names = FALSE,
#             col.names = TRUE, qmethod = c("escape", "double"),
#             fileEncoding = "")

# write.table(tr.3_utrorf.annot$ID, file = "../output/PW363/PW363_tr.3_utrorf_annot_IDs_all.tsv", 
#             append = FALSE, quote = FALSE, sep = "\t",
#             eol = "\n", na = "NA", dec = ".", row.names = FALSE,
#             col.names = FALSE, qmethod = c("escape", "double"),
#             fileEncoding = "")

write.table(tr.3_utrorf.unknown$ID, file = "../output/PW363/PW363_tr.3_utrorf_drop_IDs.tsv",
            append = FALSE, quote = FALSE, sep = "\t",
            eol = "\n", na = "NA", dec = ".", row.names = FALSE,
            col.names = FALSE, qmethod = c("escape", "double"),
            fileEncoding = "")


tr.3_utrorf.annot$status="discard"
tr.3_utrorf.unknown$status="discard"

ind = which(tr.3_utrorf.annot$tr.pass == "okay")

# write.table(tr.3_utrorf.annot[ind,], file = "../output/PW363/keep/PW363_tr.3_utrorf_annot_keep.tsv", 
#             append = FALSE, quote = FALSE, sep = "\t",
#             eol = "\n", na = "NA", dec = ".", row.names = FALSE,
#             col.names = TRUE, qmethod = c("escape", "double"),
#             fileEncoding = "")
tr.3_utrorf.annot[ind,]$status = "keep"

write.table(tr.3_utrorf.annot$ID[ind], file = "../output/PW363/PW363_tr.3_utrorf_keep_IDs.tsv",
            append = FALSE, quote = FALSE, sep = "\t",
            eol = "\n", na = "NA", dec = ".", row.names = FALSE,
            col.names = FALSE, qmethod = c("escape", "double"),
            fileEncoding = "")

# write.table(tr.3_utrorf.annot[-ind,], file = "../output/PW363/remove/PW363_tr.3_utrorf_annot_drop.tsv", 
#             append = FALSE, quote = FALSE, sep = "\t",
#             eol = "\n", na = "NA", dec = ".", row.names = FALSE,
#             col.names = TRUE, qmethod = c("escape", "double"),
#             fileEncoding = "")

tr.3_utrorf.annot[-ind,]$status = "discard"

write.table(tr.3_utrorf.annot$ID[-ind], file = "../output/PW363/PW363_tr.3_utrorf_drop_IDs.tsv", 
            append = TRUE, quote = FALSE, sep = "\t",
            eol = "\n", na = "NA", dec = ".", row.names = FALSE,
            col.names = FALSE, qmethod = c("escape", "double"),
            fileEncoding = "")
# colnames(tr.3)

okay = which(!is.na(tr.3$DM.trID))
length(okay)
## [1] 108587
okay = c(okay, which(!is.na(tr.3$IPS_Analysis)))
length(okay)
## [1] 246558
ind = grep("Target_ID", colnames(tr.3))

okay = c(okay, 
         which(!is.na(tr.3[,ind[1]])),
         which(!is.na(tr.3[,ind[2]])),
         which(!is.na(tr.3[,ind[3]])),
         which(!is.na(tr.3[,ind[4]])),
         which(!is.na(tr.3[,ind[5]])),
         which(!is.na(tr.3[,ind[6]])))

# which(as.numeric(tr.3$cds.length) >= 180)
# which(as.numeric(tr.3$cds.length) < 180)

okay = sort(unique(okay))
length(okay)
## [1] 178151
tr.3.annot = tr.3[okay,]
tr.3.unknown = tr.3[-okay,]

tr.3.annot$status="keep"
tr.3.unknown$status="discard"

hist(as.numeric(tr.3.annot$tr.length), 
     breaks = seq(0, max(as.numeric(tr.3.annot$tr.length), na.rm = TRUE) + 50, 50),
      main = "OK PW363 evigene transcript length", xlab = "transcript length",
     xlim=c(0, max(300,max(as.numeric(tr.3.annot$tr.length), na.rm = TRUE))))

hist(as.numeric(tr.3.annot$cds.length), 
     breaks = seq(0, max(as.numeric(tr.3.annot$cds.length), na.rm = TRUE) + 50, 50),
     main = "OK PW363 evigene coding sequence length", xlab = "CDS length",
     xlim=c(0, max(300,max(as.numeric(tr.3.annot$cds.length), na.rm = TRUE))))
abline(v=180, col = "red")
abline(v=200, col = "blue")
abline(v=300, col = "green")

table(as.numeric(tr.3.annot$cds.length) >= 300)
## 
##  FALSE   TRUE 
##  34511 143640
table(as.numeric(tr.3.annot$cds.length) >= 200)
## 
##  FALSE   TRUE 
##  23147 155004
table(as.numeric(tr.3.annot$cds.length) >= 180)
## 
##  FALSE   TRUE 
##  19774 158377
hist(as.numeric(tr.3.unknown$tr.length), 
     breaks = seq(0, max(as.numeric(tr.3.unknown$tr.length), na.rm = TRUE) + 50, 50),
      main = "notOK PW363 evigene transcript with secondary ORF length", xlab = "transcript length",
     xlim=c(0, max(300,max(as.numeric(tr.3.unknown$tr.length), na.rm = TRUE))))

hist(as.numeric(tr.3.unknown$cds.length), 
     breaks = seq(0, max(as.numeric(tr.3.unknown$cds.length), na.rm = TRUE) + 50, 50),
     main = "notOK PW363 evigene secondary ORF coding sequence length", xlab = "CDS length",
     xlim=c(0, max(300,max(as.numeric(tr.3.unknown$cds.length), na.rm = TRUE))),     
     col = "grey45", las=2)
abline(v=180, col = "red")
abline(v=200, col = "blue")
abline(v=300, col = "green")
axis(1, at=c(180), labels=c(180), las=2, col = c("red"), cex.axis=0.75, col.axis = c("red"))
axis(1, at=c(200), labels=c(200), las=2, col = c("blue"), cex.axis=0.75, col.axis = c("blue"))
axis(1, at=c(300), labels=c(300), las=2, col = c("green"), cex.axis=0.75, col.axis = c("green"))

table(as.numeric(tr.3.unknown$cds.length) >= 300)
## 
## FALSE  TRUE 
## 77827 13714
table(as.numeric(tr.3.unknown$cds.length) >= 200)
## 
## FALSE  TRUE 
## 65083 26458
table(as.numeric(tr.3.unknown$cds.length) >= 180)
## 
## FALSE  TRUE 
## 58548 32993
tr.3.annot[is.na(tr.3.annot)] = "-"
tr.3.unknown[is.na(tr.3.unknown)] = "-"

# write.table(tr.3.annot, file = "../output/PW363/keep/PW363_tr.3_annot.tsv", 
#             append = FALSE, quote = FALSE, sep = "\t",
#             eol = "\n", na = "NA", dec = ".", row.names = FALSE,
#             col.names = TRUE, qmethod = c("escape", "double"),
#             fileEncoding = "")
# 
# write.table(tr.3.unknown, file = "../output/PW363/remove/PW363_tr.3_unknown.tsv", 
#             append = FALSE, quote = FALSE, sep = "\t",
#             eol = "\n", na = "NA", dec = ".", row.names = FALSE,
#             col.names = TRUE, qmethod = c("escape", "double"),
#             fileEncoding = "")

write.table(tr.3.annot$ID, file = "../output/PW363/PW363_tr.3_keep_IDs.tsv", 
            append = FALSE, quote = FALSE, sep = "\t",
            eol = "\n", na = "NA", dec = ".", row.names = FALSE,
            col.names = FALSE, qmethod = c("escape", "double"),
            fileEncoding = "")

write.table(tr.3.unknown$ID, file = "../output/PW363/PW363_tr.3_drop_IDs.tsv", 
            append = FALSE, quote = FALSE, sep = "\t",
            eol = "\n", na = "NA", dec = ".", row.names = FALSE,
            col.names = FALSE, qmethod = c("escape", "double"),
            fileEncoding = "")

combo table

tr.3.unknown$ID2 = tr.3.unknown$ID
tr.3.annot$ID2 = tr.3.annot$ID
orderC = match(colnames(tr.3_utrorf.annot), colnames(tr.3.annot))
combo = rbind(tr.3.unknown[,orderC], tr.3_utrorf.unknown, tr.3_utrorf.annot, tr.3.annot[,orderC])
nrow(combo)==nrow(tr.3_utrorf)+nrow(tr.3)
## [1] TRUE
colnames(combo)[1:2]=c("tr.ID", "cds.ID")

write.table(combo, file = "../output/PW363/PW363_tr.cds.tsv", 
            append = FALSE, quote = FALSE, sep = "\t",
            eol = "\n", na = "NA", dec = ".", row.names = FALSE,
            col.names = TRUE, qmethod = c("escape", "double"),
            fileEncoding = "")
# colnames(tr.3.annot)

take = c(which(is.na(no.gene_score.0$IPS_Analysis)))
length(take)
## [1] 58841
ind = grep("Target_ID", colnames(no.gene_score.0))

take = intersect(take, (
         which(is.na(no.gene_score.0[,ind[1]]) &
         (is.na(no.gene_score.0[,ind[2]])) &
         (is.na(no.gene_score.0[,ind[3]])) &
         (is.na(no.gene_score.0[,ind[4]])) &
         (is.na(no.gene_score.0[,ind[5]])) &
         (is.na(no.gene_score.0[,ind[6]])))))


take = sort(unique(take))
length(take)
## [1] 56502
no.gene_score.0 = no.gene_score.0[take,]
no.gene_score.0 = no.gene_score.0[, 1:which(colnames(no.gene_score.0) == "blastnVector")]
no.gene_score.0 = no.gene_score.0[, -c(grep("IPR", colnames(no.gene_score.0)), grep("IPS", colnames(no.gene_score.0)))]

write.table(no.gene_score.0, file = "../output/PW363/PW363_justDMmapped_score.0.tsv", 
            append = FALSE, quote = FALSE, sep = "\t",
            eol = "\n", na = "NA", dec = ".", row.names = FALSE,
            col.names = TRUE, qmethod = c("escape", "double"),
            fileEncoding = "")
sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17134)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Slovenian_Slovenia.1250  LC_CTYPE=Slovenian_Slovenia.1250   
## [3] LC_MONETARY=Slovenian_Slovenia.1250 LC_NUMERIC=C                       
## [5] LC_TIME=Slovenian_Slovenia.1250    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] grr_0.9.5         forcats_0.4.0     stringr_1.4.0    
##  [4] dplyr_0.8.3       purrr_0.3.2       readr_1.3.1      
##  [7] tidyr_0.8.3       tibble_2.1.3      ggplot2_3.2.1    
## [10] tidyverse_1.2.1   data.table_1.12.2
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.2       cellranger_1.1.0 pillar_1.4.2     compiler_3.6.1  
##  [5] tools_3.6.1      zeallot_0.1.0    digest_0.6.20    lubridate_1.7.4 
##  [9] jsonlite_1.6     evaluate_0.14    nlme_3.1-141     gtable_0.3.0    
## [13] lattice_0.20-38  pkgconfig_2.0.2  rlang_0.4.0      cli_1.1.0       
## [17] rstudioapi_0.10  yaml_2.2.0       haven_2.1.1      xfun_0.8        
## [21] withr_2.1.2      xml2_1.2.2       httr_1.4.1       knitr_1.24      
## [25] vctrs_0.2.0      hms_0.5.1        generics_0.0.2   grid_3.6.1      
## [29] tidyselect_0.2.5 glue_1.3.1       R6_2.4.0         readxl_1.3.1    
## [33] rmarkdown_1.14   modelr_0.1.5     magrittr_1.5     backports_1.1.4 
## [37] scales_1.0.0     htmltools_0.3.6  rvest_0.3.4      assertthat_0.2.1
## [41] colorspace_1.4-1 stringi_1.4.3    lazyeval_0.2.2   munsell_0.5.0   
## [45] broom_0.5.2      crayon_1.3.4
rm(list=ls())
gc()
##           used (Mb) gc trigger  (Mb)  max used  (Mb)
## Ncells  942840 50.4    3481103 186.0   5439222 290.5
## Vcells 3715039 28.4  106435075 812.1 114993657 877.4