Mercurial > repos > davidvanzessen > shm_csr
view sequence_overview.r @ 17:b95fa7e426c3 draft
Uploaded
author | davidvanzessen |
---|---|
date | Thu, 24 Nov 2016 10:24:19 -0500 |
parents | 3968d04b5724 |
children | fe44a905aee9 |
line wrap: on
line source
library(reshape2) args <- commandArgs(trailingOnly = TRUE) before.unique.file = args[1] merged.file = args[2] outputdir = args[3] gene.classes = unlist(strsplit(args[4], ",")) hotspot.analysis.sum.file = args[5] NToverview.file = paste(outputdir, "ntoverview.txt", sep="/") NTsum.file = paste(outputdir, "ntsum.txt", sep="/") main.html = "index.html" empty.region.filter = args[6] setwd(outputdir) before.unique = read.table(before.unique.file, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="") merged = read.table(merged.file, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="") hotspot.analysis.sum = read.table(hotspot.analysis.sum.file, header=F, sep=",", fill=T, stringsAsFactors=F, quote="") #before.unique = before.unique[!grepl("unmatched", before.unique$best_match),] if(empty.region.filter == "leader"){ before.unique$seq_conc = paste(before.unique$FR1.IMGT.seq, before.unique$CDR1.IMGT.seq, before.unique$FR2.IMGT.seq, before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq) } else if(empty.region.filter == "FR1"){ before.unique$seq_conc = paste(before.unique$CDR1.IMGT.seq, before.unique$FR2.IMGT.seq, before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq) } else if(empty.region.filter == "CDR1"){ before.unique$seq_conc = paste(before.unique$FR2.IMGT.seq, before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq) } else if(empty.region.filter == "FR2"){ before.unique$seq_conc = paste(before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq) } IDs = before.unique[,c("Sequence.ID", "seq_conc", "best_match", "Functionality")] IDs$best_match = as.character(IDs$best_match) dat = data.frame(table(before.unique$seq_conc)) names(dat) = c("seq_conc", "Freq") dat$seq_conc = factor(dat$seq_conc) dat = dat[order(as.character(dat$seq_conc)),] #writing html from R... get.bg.color = function(val){ if(val %in% c("TRUE", "FALSE", "T", "F")){ #if its a logical value, give the background a green/red color return(ifelse(val,"#eafaf1","#f9ebea")) } else if (!is.na(as.numeric(val))) { #if its a numerical value, give it a grey tint if its >0 return(ifelse(val > 0,"#eaecee","white")) } else { return("white") } } td = function(val) { return(paste("<td bgcolor='", get.bg.color(val), "'>", val, "</td>", sep="")) } tr = function(val) { return(paste(c("<tr>", sapply(val, td), "</tr>"), collapse="")) } make.link = function(id, clss, val) { paste("<a href='", clss, "_", id, ".html'>", val, "</a>", sep="") } tbl = function(df) { res = "<table border='1'>" for(i in 1:nrow(df)){ res = paste(res, tr(df[i,]), sep="") } res = paste(res, "</table>") } cat("<table border='1' class='pure-table pure-table-striped'>", file=main.html, append=F) if(empty.region.filter == "leader"){ cat("<caption>FR1+CDR1+FR2+CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T) } else if(empty.region.filter == "FR1"){ cat("<caption>CDR1+FR2+CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T) } else if(empty.region.filter == "CDR1"){ cat("<caption>FR2+CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T) } else if(empty.region.filter == "FR2"){ cat("<caption>CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T) } cat("<tr>", file=main.html, append=T) cat("<th>Sequence</th><th>Functionality</th><th>ca1</th><th>ca2</th><th>cg1</th><th>cg2</th><th>cg3</th><th>cg4</th><th>cm</th><th>un</th>", file=main.html, append=T) cat("<th>total CA</th><th>total CG</th><th>number of subclasses</th><th>present in both Ca and Cg</th><th>Ca1+Ca2</th>", file=main.html, append=T) cat("<th>Cg1+Cg2</th><th>Cg1+Cg3</th><th>Cg1+Cg4</th><th>Cg2+Cg3</th><th>Cg2+Cg4</th><th>Cg3+Cg4</th>", file=main.html, append=T) cat("<th>Cg1+Cg2+Cg3</th><th>Cg2+Cg3+Cg4</th><th>Cg1+Cg2+Cg4</th><th>Cg1+Cg3+Cg4</th><th>Cg1+Cg2+Cg3+Cg4</th>", file=main.html, append=T) cat("</tr>", file=main.html, append=T) single.sequences=0 #sequence only found once, skipped in.multiple=0 #same sequence across multiple subclasses multiple.in.one=0 #same sequence multiple times in one subclass unmatched=0 #all of the sequences are unmatched some.unmatched=0 #one or more sequences in a clone are unmatched matched=0 #should be the same als matched sequences sequence.id.page="by_id.html" for(i in 1:nrow(dat)){ ca1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGA1", IDs$best_match),] ca2 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGA2", IDs$best_match),] cg1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGG1", IDs$best_match),] cg2 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGG2", IDs$best_match),] cg3 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGG3", IDs$best_match),] cg4 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGG4", IDs$best_match),] cm = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGM", IDs$best_match),] un = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^unmatched", IDs$best_match),] allc = rbind(ca1, ca2, cg1, cg2, cg3, cg4, cm, un) ca1.n = nrow(ca1) ca2.n = nrow(ca2) cg1.n = nrow(cg1) cg2.n = nrow(cg2) cg3.n = nrow(cg3) cg4.n = nrow(cg4) cm.n = nrow(cm) un.n = nrow(un) classes = c(ca1.n, ca2.n, cg1.n, cg2.n, cg3.n, cg4.n, cm.n, un.n) classes.sum = sum(classes) if(classes.sum == 1){ single.sequences = single.sequences + 1 next } if(un.n == classes.sum){ unmatched = unmatched + 1 next } in.classes = sum(classes > 0) matched = matched + in.classes #count in how many subclasses the sequence occurs. if(any(classes == classes.sum)){ multiple.in.one = multiple.in.one + 1 } else if (un.n > 0) { some.unmatched = some.unmatched + 1 } else { in.multiple = in.multiple + 1 } id = as.numeric(dat[i,"seq_conc"]) functionality = paste(unique(allc[,"Functionality"]), collapse=",") by.id.row = c() if(ca1.n > 0){ cat(tbl(ca1), file=paste("IGA1_", id, ".html", sep="")) } if(ca2.n > 0){ cat(tbl(ca2), file=paste("IGA2_", id, ".html", sep="")) } if(cg1.n > 0){ cat(tbl(cg1), file=paste("IGG1_", id, ".html", sep="")) } if(cg2.n > 0){ cat(tbl(cg2), file=paste("IGG2_", id, ".html", sep="")) } if(cg3.n > 0){ cat(tbl(cg3), file=paste("IGG3_", id, ".html", sep="")) } if(cg4.n > 0){ cat(tbl(cg4), file=paste("IGG4_", id, ".html", sep="")) } if(cm.n > 0){ cat(tbl(cm), file=paste("IGM_", id, ".html", sep="")) } if(un.n > 0){ cat(tbl(un), file=paste("un_", id, ".html", sep="")) } ca1.html = make.link(id, "IGA1", ca1.n) ca2.html = make.link(id, "IGA2", ca2.n) cg1.html = make.link(id, "IGG1", cg1.n) cg2.html = make.link(id, "IGG2", cg2.n) cg3.html = make.link(id, "IGG3", cg3.n) cg4.html = make.link(id, "IGG4", cg4.n) cm.html = make.link(id, "IGM", cm.n) un.html = make.link(id, "un", un.n) #extra columns ca.n = ca1.n + ca2.n cg.n = cg1.n + cg2.n + cg3.n + cg4.n #in.classes in.ca.cg = (ca.n > 0 & cg.n > 0) in.ca1.ca2 = (ca1.n > 0 & ca2.n > 0) in.cg1.cg2 = (cg1.n > 0 & cg2.n > 0) in.cg1.cg3 = (cg1.n > 0 & cg3.n > 0) in.cg1.cg4 = (cg1.n > 0 & cg4.n > 0) in.cg2.cg3 = (cg2.n > 0 & cg3.n > 0) in.cg2.cg4 = (cg2.n > 0 & cg4.n > 0) in.cg3.cg4 = (cg3.n > 0 & cg4.n > 0) in.cg1.cg2.cg3 = (cg1.n > 0 & cg2.n > 0 & cg3.n > 0) in.cg2.cg3.cg4 = (cg2.n > 0 & cg3.n > 0 & cg4.n > 0) in.cg1.cg2.cg4 = (cg1.n > 0 & cg2.n > 0 & cg4.n > 0) in.cg1.cg3.cg4 = (cg1.n > 0 & cg3.n > 0 & cg4.n > 0) in.cg.all = (cg1.n > 0 & cg2.n > 0 & cg3.n > 0 & cg4.n > 0) #rw = c(as.character(dat[i,"seq_conc"]), functionality, ca1.html, ca2.html, cg1.html, cg2.html, cg3.html, cg4.html, cm.html, un.html) rw = c(as.character(dat[i,"seq_conc"]), functionality, ca1.html, ca2.html, cg1.html, cg2.html, cg3.html, cg4.html, cm.html, un.html) rw = c(rw, ca.n, cg.n, in.classes, in.ca.cg, in.ca1.ca2, in.cg1.cg2, in.cg1.cg3, in.cg1.cg4, in.cg2.cg3, in.cg2.cg4, in.cg3.cg4, in.cg1.cg2.cg3, in.cg2.cg3.cg4, in.cg1.cg2.cg4, in.cg1.cg3.cg4, in.cg.all) cat(tr(rw), file=main.html, append=T) for(i in 1:nrow(allc)){ #generate html by id html = make.link(id, allc[i,"best_match"], allc[i,"Sequence.ID"]) cat(paste(html, "<br />"), file=sequence.id.page, append=T) } } cat("</table>", file=main.html, append=T) print(paste("Single sequences:", single.sequences)) print(paste("Sequences in multiple subclasses:", in.multiple)) print(paste("Multiple sequences in one subclass:", multiple.in.one)) print(paste("Matched with unmatched:", some.unmatched)) print(paste("Count that should match 'matched' sequences:", matched)) #ACGT overview #NToverview = merged[!grepl("^unmatched", merged$best_match),] NToverview = merged if(empty.region.filter == "leader"){ NToverview$seq = paste(NToverview$FR1.IMGT.seq, NToverview$CDR1.IMGT.seq, NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq) } else if(empty.region.filter == "FR1"){ NToverview$seq = paste(NToverview$CDR1.IMGT.seq, NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq) } else if(empty.region.filter == "CDR1"){ NToverview$seq = paste(NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq) } else if(empty.region.filter == "FR2"){ NToverview$seq = paste(NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq) } NToverview$A = nchar(gsub("[^Aa]", "", NToverview$seq)) NToverview$C = nchar(gsub("[^Cc]", "", NToverview$seq)) NToverview$G = nchar(gsub("[^Gg]", "", NToverview$seq)) NToverview$T = nchar(gsub("[^Tt]", "", NToverview$seq)) #Nsum = data.frame(Sequence.ID="-", best_match="Sum", seq="-", A = sum(NToverview$A), C = sum(NToverview$C), G = sum(NToverview$G), T = sum(NToverview$T)) #NToverview = rbind(NToverview, NTsum) NTresult = data.frame(nt=c("A", "C", "T", "G")) for(clazz in gene.classes){ print(paste("class:", clazz)) NToverview.sub = NToverview[grepl(paste("^", clazz, sep=""), NToverview$best_match),] print(paste("nrow:", nrow(NToverview.sub))) new.col.x = c(sum(NToverview.sub$A), sum(NToverview.sub$C), sum(NToverview.sub$T), sum(NToverview.sub$G)) new.col.y = sum(new.col.x) new.col.z = round(new.col.x / new.col.y * 100, 2) tmp = names(NTresult) NTresult = cbind(NTresult, data.frame(new.col.x, new.col.y, new.col.z)) names(NTresult) = c(tmp, paste(clazz, c("x", "y", "z"), sep="")) } write.table(NToverview[,c("Sequence.ID", "best_match", "seq", "A", "C", "G", "T")], NToverview.file, quote=F, sep="\t", row.names=F, col.names=T) NToverview = NToverview[!grepl("unmatched", NToverview$best_match),] new.col.x = c(sum(NToverview$A), sum(NToverview$C), sum(NToverview$T), sum(NToverview$G)) new.col.y = sum(new.col.x) new.col.z = round(new.col.x / new.col.y * 100, 2) tmp = names(NTresult) NTresult = cbind(NTresult, data.frame(new.col.x, new.col.y, new.col.z)) names(NTresult) = c(tmp, paste("all", c("x", "y", "z"), sep="")) names(hotspot.analysis.sum) = names(NTresult) hotspot.analysis.sum = rbind(hotspot.analysis.sum, NTresult) write.table(hotspot.analysis.sum, hotspot.analysis.sum.file, quote=F, sep=",", row.names=F, col.names=F, na="0")