knitr::opts_chunk$set(echo = TRUE)
oldWC = read.table(file = "../intermediate_2_prepare-network/3cv_evigene-header-network_unfiltered_weak-components_perWC-filtered_with-rep.txt",
header = FALSE,
sep = "\t", quote = NULL,
dec = ".",
stringsAsFactors = FALSE,
comment.char = "#")
head(oldWC)
colnames(oldWC) = c("stCuSTr.wcID", "cv", "geneCnt", "geneID", "status")
newWC = cbind(oldWC, rep(NA, nrow(oldWC)), rep(NA, nrow(oldWC)))
dim(oldWC)
dim(newWC)
colnames(newWC) = c("stCuSTr.wcID", "cv", "geneCnt", "geneID", "status", "cdhit2d.wcID", "cdhit2d.status")
unique(newWC$status)
newWC$cdhit2d.status[newWC$status == "utrorf"] = "utrorf"
newWC$cdhit2d.status[newWC$status == "representative"] = "representative"
interWC = read.table(file = "../output/fasta_1B_alt-per-cv/stCuSTr-D_aa_altKept.3cols",
header = FALSE,
sep = "\t", quote = NULL,
dec = ".",
stringsAsFactors = FALSE,
comment.char = "#")
colnames(interWC) = c("cdhit2d.wcID", "aaLen", "geneID")
interWC$cdhit2d.wcID = paste0("stCuSTr-D_",formatC(interWC$cdhit2d.wcID, width = 5, format = "d", flag = "0"))
head(interWC)
ind = match(interWC$geneID, newWC$geneID)
any(is.na(ind))
newWC$cdhit2d.wcID[ind] = interWC$cdhit2d.wcID
newWC$cdhit2d.status[ind] = "alternative"
interWC = read.table(file = "../output/fasta_1A_rep-per-cv/stCuSTr-D_aa_rep.3cols",
header = FALSE,
sep = "\t", quote = NULL,
dec = ".",
stringsAsFactors = FALSE,
comment.char = "#")
colnames(interWC) = c("cdhit2d.wcID", "aaLen", "geneID")
interWC$cdhit2d.wcID = paste0("stCuSTr-D_",formatC(interWC$cdhit2d.wcID, width = 5, format = "d", flag = "0"))
head(interWC)
ind = match(interWC$geneID, newWC$geneID)
any(is.na(ind))
newWC$cdhit2d.wcID[ind] = interWC$cdhit2d.wcID
PW363
interWC = read.table(file = "../output/fasta_1B_alt-per-cv/stCuSTr-P_aa_altKept.3cols",
header = FALSE,
sep = "\t", quote = NULL,
dec = ".",
stringsAsFactors = FALSE,
comment.char = "#")
colnames(interWC) = c("cdhit2d.wcID", "aaLen", "geneID")
interWC$cdhit2d.wcID = paste0("stCuSTr-P_",formatC(interWC$cdhit2d.wcID, width = 5, format = "d", flag = "0"))
head(interWC)
ind = match(interWC$geneID, newWC$geneID)
any(is.na(ind))
newWC$cdhit2d.wcID[ind] = interWC$cdhit2d.wcID
newWC$cdhit2d.status[ind] = "alternative"
interWC = read.table(file = "../output/fasta_1A_rep-per-cv/stCuSTr-P_aa_rep.3cols",
header = FALSE,
sep = "\t", quote = NULL,
dec = ".",
stringsAsFactors = FALSE,
comment.char = "#")
colnames(interWC) = c("cdhit2d.wcID", "aaLen", "geneID")
interWC$cdhit2d.wcID = paste0("stCuSTr-P_",formatC(interWC$cdhit2d.wcID, width = 5, format = "d", flag = "0"))
head(interWC)
ind = match(interWC$geneID, newWC$geneID)
any(is.na(ind))
newWC$cdhit2d.wcID[ind] = interWC$cdhit2d.wcID
interWC = read.table(file = "../output/fasta_1B_alt-per-cv/stCuSTr-R_aa_altKept.3cols",
header = FALSE,
sep = "\t", quote = NULL,
dec = ".",
stringsAsFactors = FALSE,
comment.char = "#")
colnames(interWC) = c("cdhit2d.wcID", "aaLen", "geneID")
interWC$cdhit2d.wcID = paste0("stCuSTr-R_",formatC(interWC$cdhit2d.wcID, width = 5, format = "d", flag = "0"))
head(interWC)
ind = match(interWC$geneID, newWC$geneID)
any(is.na(ind))
newWC$cdhit2d.wcID[ind] = interWC$cdhit2d.wcID
newWC$cdhit2d.status[ind] = "alternative"
interWC = read.table(file = "../output/fasta_1A_rep-per-cv/stCuSTr-R_aa_rep.3cols",
header = FALSE,
sep = "\t", quote = NULL,
dec = ".",
stringsAsFactors = FALSE,
comment.char = "#")
colnames(interWC) = c("cdhit2d.wcID", "aaLen", "geneID")
interWC$cdhit2d.wcID = paste0("stCuSTr-R_",formatC(interWC$cdhit2d.wcID, width = 5, format = "d", flag = "0"))
head(interWC)
ind = match(interWC$geneID, newWC$geneID)
any(is.na(ind))
newWC$cdhit2d.wcID[ind] = interWC$cdhit2d.wcID
newWC$cdhit2d.status[is.na(newWC$cdhit2d.status)] = "eliminated"
newWC$cdhit2d.geneCnt = 1
tmp = aggregate(newWC$cdhit2d.geneCnt, by = list(y = newWC$cdhit2d.wcID), FUN = sum)
newWC$cdhit2d.geneCnt = tmp$x[match(newWC$cdhit2d.wcID, tmp$y)]
newWC$cdhit2d.geneCnt[is.na(newWC$cdhit2d.geneCnt)] = "NAp"
kept = newWC[(newWC$cdhit2d.status != "utrorf") & (newWC$cdhit2d.status != "eliminated"),]
translation = cbind(kept$stCuSTr.wcID[kept$status == "representative"],
kept$cdhit2d.wcID[kept$status == "representative"])
kept$cdhit2d.wcID = translation[,1][match(kept$cdhit2d.wcID, translation[,2])]
kept = rbind(kept, newWC[(newWC$cdhit2d.status == "utrorf"),])
kept = rbind(kept, newWC[(newWC$cdhit2d.status == "eliminated"),])
kept$match.wc = (kept$stCuSTr.wcID == kept$cdhit2d.wcID)
kept$match.wc[is.na(kept$match.wc)] = "NAp"
unique(kept$match.wc)
sum(kept$match.wc == TRUE)
sum(kept$match.wc == FALSE)
sum(kept$match.wc == "NAp")
kept = kept[with(kept, order(rev(match.wc), rev(cdhit2d.status), cdhit2d.wcID)), ]
write.table(kept,
file = "../output/3cv_cdhit-2d_weak-components.txt",
append = FALSE, quote = FALSE, sep = "\t",
eol = "\n", na = "NA", dec = ".", row.names = FALSE,
col.names = TRUE)
kept$cdhit2d.wcID[is.na(kept$cdhit2d.wcID)] = kept$stCuSTr.wcID[is.na(kept$cdhit2d.wcID)]
i1 = which(kept$cdhit2d.geneCnt == "NAp")
i2 = which(kept$cdhit2d.geneCnt != "NAp")
kept$cdhit2d.geneCnt[i1] = ""
alias = cbind(kept$geneID[i2],
paste0(kept$geneID[i2], " [", kept$cdhit2d.wcID[i2], ", ",
kept$cdhit2d.status[i2], ", #", kept$cdhit2d.geneCnt[i2], "]"))
alias = rbind(alias,
cbind(kept$geneID[i1],
paste0(kept$geneID[i1], " [", kept$cdhit2d.wcID[i1], ", ",
kept$cdhit2d.status[i1], "]")))
write.table(alias,
file = "../output/3cv_cdhit-2d_weak-components_alias.txt",
append = FALSE, quote = FALSE, sep = "\t",
eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE)
sessionInfo()