Mercurial > repos > davidvanzessen > mutation_analysis
view sequence_overview.r @ 81:a778156dad3d draft
Uploaded
author | davidvanzessen |
---|---|
date | Tue, 17 May 2016 04:17:38 -0400 |
parents | 0513b46178c4 |
children | 564c4f6da203 |
line wrap: on
line source
library(reshape2) args <- commandArgs(trailingOnly = TRUE) gene.matches = args[1] sequence.file = args[2] outputdir = args[3] NToverview.file = paste(outputdir, "ntoverview.txt", sep="/") NTsum.file = paste(outputdir, "ntsum.txt", sep="/") main.html = "index.html" setwd(outputdir) genes = read.table(gene.matches, header=T, sep="\t", fill=T) sequences = read.table(sequence.file, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="") dat = merge(sequences, genes, by="Sequence.ID") dat$seq_conc = paste(dat$CDR1.IMGT, dat$CDR2.IMGT, dat$CDR3.IMGT, dat$FR2.IMGT, dat$FR3.IMGT) IDs = dat[,c("Sequence.ID", "seq_conc", "best_match", "Functionality")] IDs$best_match = as.character(IDs$best_match) #dat = data.frame(data.table(dat)[, list(freq=.N), by=c("best_match", "seq_conc")]) dat = data.frame(table(dat$seq_conc)) dat = dat[dat$Freq > 1,] names(dat) = c("seq_conc", "Freq") dat$seq_conc = factor(dat$seq_conc) dat = dat[order(nchar(as.character(dat$seq_conc))),] #writing html from R... td = function(val) { paste("<td>", val, "</td>", sep="") } tr = function(val) { capture.output(cat("<tr>", td(val), "</tr>", sep="")) } 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'>", file=main.html, append=F) cat("<caption>CDR1+CDR2+CDR3+FR2+FR3 sequences that show up more than once</caption>", file=main.html, append=T) cat("<tr><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></tr>", file=main.html, append=T) for(i in 1:nrow(dat)){ ca1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "ca1",] ca2 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "ca2",] cg1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "cg1",] cg2 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "cg2",] cg3 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "cg3",] cg4 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "cg4",] cm = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "cm",] id = as.numeric(dat[i,"seq_conc"]) functionality = dat[i,"Functionality"] if(nrow(ca1) > 0){ cat(tbl(ca1), file=paste("ca1_", id, ".html", sep="")) } if(nrow(ca2) > 0){ cat(tbl(ca2), file=paste("ca2_", id, ".html", sep="")) } if(nrow(cg1) > 0){ cat(tbl(cg1), file=paste("cg1_", id, ".html", sep="")) } if(nrow(cg2) > 0){ cat(tbl(cg2), file=paste("cg2_", id, ".html", sep="")) } if(nrow(cg3) > 0){ cat(tbl(cg3), file=paste("cg3_", id, ".html", sep="")) } if(nrow(cg4) > 0){ cat(tbl(cg4), file=paste("cg4_", id, ".html", sep="")) } if(nrow(cm) > 0){ cat(tbl(cm), file=paste("cm_", id, ".html", sep="")) } ca1.html = make.link(id, "ca1", nrow(ca1)) ca2.html = make.link(id, "ca2", nrow(ca2)) cg1.html = make.link(id, "cg1", nrow(cg1)) cg2.html = make.link(id, "cg2", nrow(cg2)) cg3.html = make.link(id, "cg3", nrow(cg3)) cg4.html = make.link(id, "cg4", nrow(cg4)) cm.html = make.link(id, "cm", nrow(cm)) rw = c(as.character(dat[i,"seq_conc"]), as.character(functionality), ca1.html, ca2.html, cg1.html, cg2.html, cg3.html, cg4.html, cm.html) cat(tr(rw), file=main.html, append=T) } cat("</table>", file=main.html, append=T) #ACGT overview NToverview = genes[,c("Sequence.ID", "best_match")] sequences$seq = paste(sequences$CDR2.IMGT, sequences$CDR2.IMGT, sequences$FR2.IMGT, sequences$FR3.IMGT, sep="_") NToverview = merge(NToverview, sequences[,c("Sequence.ID", "seq")], by="Sequence.ID") 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)) NTsum = data.frame(Sequence.ID="-", best_match="Sum", seq="-", A = sum(NToverview$A), C = sum(NToverview$C), G = sum(NToverview$G), T = sum(NToverview$T)) print(names(NToverview)) print(names(NTsum)) NToverview = rbind(NToverview, NTsum) write.table(NToverview, NToverview.file, quote=F, sep="\t", row.names=F, col.names=T) #write.table(NTsum, NTsum.file, quote=F, sep="\t", row.names=F, col.names=T)