Mercurial > repos > davidvanzessen > mutation_analysis

view mutation_analysis.r @ 21:c9f9623f1f76 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Thu, 02 Apr 2015 03:31:23 -0400
parents cb7c65e3e43f
children d84c9791d8c4
line wrap: on
line source

args <- commandArgs(trailingOnly = TRUE)

input = args[1]
genes = unlist(strsplit(args[2], ","))
outputdir = args[3]
setwd(outputdir)

dat = read.table(input, header=T, sep="\t", fill=T, stringsAsFactors=F)

if(length(dat$Sequence.ID) == 0){
  setwd(outputdir)
  result = data.frame(x = rep(0, 5), y = rep(0, 5), z = rep(NA, 5))
  row.names(result) = c("Number of Mutations (%)", "Transition (%)", "Transversions (%)", "Transitions at G C (%)", "Targeting of C G (%)")
  write.table(x=result, file="mutations.txt", sep=",",quote=F,row.names=T,col.names=F)
  transitionTable = data.frame(A=rep(0, 4),C=rep(0, 4),G=rep(0, 4),T=rep(0, 4))
  row.names(transitionTable) = c("A", "C", "G", "T")
  transitionTable["A","A"] = NA
  transitionTable["C","C"] = NA
  transitionTable["G","G"] = NA
  transitionTable["T","T"] = NA
  write.table(x=transitionTable, file="transitions.txt", sep=",",quote=F,row.names=T,col.names=NA)
  cat("0", file="n.txt")
  stop("No data")
}



cleanup_columns = c("FR1.IMGT.c.a",
                    "FR2.IMGT.g.t",
                    "CDR1.IMGT.Nb.of.nucleotides",
                    "CDR2.IMGT.t.a",
                    "FR1.IMGT.c.g",
                    "CDR1.IMGT.c.t",
                    "FR2.IMGT.a.c",
                    "FR2.IMGT.Nb.of.mutations",
                    "FR2.IMGT.g.c",
                    "FR2.IMGT.a.g",
                    "FR3.IMGT.t.a",
                    "FR3.IMGT.t.c",
                    "FR2.IMGT.g.a",
                    "FR3.IMGT.c.g",
                    "FR1.IMGT.Nb.of.mutations",
                    "CDR1.IMGT.g.a",
                    "CDR1.IMGT.t.g",
                    "CDR1.IMGT.g.c",
                    "CDR2.IMGT.Nb.of.nucleotides",
                    "FR2.IMGT.a.t",
                    "CDR1.IMGT.Nb.of.mutations",
                    "CDR1.IMGT.a.g",
                    "FR3.IMGT.a.c",
                    "FR1.IMGT.g.a",
                    "FR3.IMGT.a.g",
                    "FR1.IMGT.a.t",
                    "CDR2.IMGT.a.g",
                    "CDR2.IMGT.Nb.of.mutations",
                    "CDR2.IMGT.g.t",
                    "CDR2.IMGT.a.c",
                    "CDR1.IMGT.t.c",
                    "FR3.IMGT.g.c",
                    "FR1.IMGT.g.t",
                    "FR3.IMGT.g.t",
                    "CDR1.IMGT.a.t",
                    "FR1.IMGT.a.g",
                    "FR3.IMGT.a.t",
                    "FR3.IMGT.Nb.of.nucleotides",
                    "FR2.IMGT.t.c",
                    "CDR2.IMGT.g.a",
                    "FR2.IMGT.t.a",
                    "CDR1.IMGT.t.a",
                    "FR2.IMGT.t.g",
                    "FR3.IMGT.t.g",
                    "FR2.IMGT.Nb.of.nucleotides",
                    "FR1.IMGT.t.a",
                    "FR1.IMGT.t.g",
                    "FR3.IMGT.c.t",
                    "FR1.IMGT.t.c",
                    "CDR2.IMGT.a.t",
                    "FR2.IMGT.c.t",
                    "CDR1.IMGT.g.t",
                    "CDR2.IMGT.t.g",
                    "FR1.IMGT.Nb.of.nucleotides",
                    "CDR1.IMGT.c.g",
                    "CDR2.IMGT.t.c",
                    "FR3.IMGT.g.a",
                    "CDR1.IMGT.a.c",
                    "FR2.IMGT.c.a",
                    "FR3.IMGT.Nb.of.mutations",
                    "FR2.IMGT.c.g",
                    "CDR2.IMGT.g.c",
                    "FR1.IMGT.g.c",
                    "CDR2.IMGT.c.t",
                    "FR3.IMGT.c.a",
                    "CDR1.IMGT.c.a",
                    "CDR2.IMGT.c.g",
                    "CDR2.IMGT.c.a",
                    "FR1.IMGT.c.t")

for(col in cleanup_columns){
  dat[,col] = gsub("\\(.*\\)", "", dat[,col])
  #dat[dat[,col] == "",] = "0"
  dat[,col] = as.numeric(dat[,col])
  dat[is.na(dat[,col]),] = 0
}

dat$VRegionMutations =  dat$CDR1.IMGT.Nb.of.mutations + 
                        dat$FR2.IMGT.Nb.of.mutations + 
                        dat$CDR2.IMGT.Nb.of.mutations +
                        dat$FR3.IMGT.Nb.of.mutations

dat$VRegionNucleotides =  dat$CDR1.IMGT.Nb.of.nucleotides +
                          dat$FR2.IMGT.Nb.of.nucleotides +
                          dat$CDR2.IMGT.Nb.of.nucleotides +
                          dat$FR3.IMGT.Nb.of.nucleotides

dat$transitionMutations = dat$CDR1.IMGT.a.g +
                          dat$CDR1.IMGT.g.a +
                          dat$CDR1.IMGT.c.t +
                          dat$CDR1.IMGT.t.c +
                          dat$FR2.IMGT.a.g +
                          dat$FR2.IMGT.g.a +
                          dat$FR2.IMGT.c.t +
                          dat$FR2.IMGT.t.c +
                          dat$CDR2.IMGT.a.g +
                          dat$CDR2.IMGT.g.a +
                          dat$CDR2.IMGT.c.t +
                          dat$CDR2.IMGT.t.c +
                          dat$FR3.IMGT.a.g +
                          dat$FR3.IMGT.g.a +
                          dat$FR3.IMGT.c.t +
                          dat$FR3.IMGT.t.c

dat$transversionMutations = dat$CDR1.IMGT.a.c +
                            dat$CDR1.IMGT.c.a +
                            dat$CDR1.IMGT.a.t +
                            dat$CDR1.IMGT.t.a +
                            dat$CDR1.IMGT.g.c +
                            dat$CDR1.IMGT.c.g +
                            dat$CDR1.IMGT.g.t +
                            dat$CDR1.IMGT.t.g +
                            dat$FR2.IMGT.a.c +
                            dat$FR2.IMGT.c.a +
                            dat$FR2.IMGT.a.t +
                            dat$FR2.IMGT.t.a +
                            dat$FR2.IMGT.g.c +
                            dat$FR2.IMGT.c.g +
                            dat$FR2.IMGT.g.t +
                            dat$FR2.IMGT.t.g +
                            dat$CDR2.IMGT.a.c +
                            dat$CDR2.IMGT.c.a +
                            dat$CDR2.IMGT.a.t +
                            dat$CDR2.IMGT.t.a +
                            dat$CDR2.IMGT.g.c +
                            dat$CDR2.IMGT.c.g +
                            dat$CDR2.IMGT.g.t +
                            dat$CDR2.IMGT.t.g +
                            dat$FR3.IMGT.a.c +
                            dat$FR3.IMGT.c.a +
                            dat$FR3.IMGT.a.t +
                            dat$FR3.IMGT.t.a +
                            dat$FR3.IMGT.g.c +
                            dat$FR3.IMGT.c.g +
                            dat$FR3.IMGT.g.t +
                            dat$FR3.IMGT.t.g


dat$transitionMutationsAtGC = dat$CDR1.IMGT.g.a +
                              dat$CDR1.IMGT.c.t +
                              dat$FR2.IMGT.g.a +
                              dat$FR2.IMGT.c.t +
                              dat$CDR2.IMGT.g.a +
                              dat$CDR2.IMGT.c.t +
                              dat$FR3.IMGT.g.a +
                              dat$FR3.IMGT.c.t

dat$totalMutationsAtGC = dat$CDR1.IMGT.g.a +
                         dat$CDR1.IMGT.c.t +
                         dat$CDR1.IMGT.c.a +
                         dat$CDR1.IMGT.g.c +
                         dat$CDR1.IMGT.c.g +
                         dat$CDR1.IMGT.g.t +
                         dat$FR2.IMGT.g.a +
                         dat$FR2.IMGT.c.t +
                         dat$FR2.IMGT.c.a +
                         dat$FR2.IMGT.g.c +
                         dat$FR2.IMGT.c.g +
                         dat$FR2.IMGT.g.t +
                         dat$CDR2.IMGT.g.a +
                         dat$CDR2.IMGT.c.t +
                         dat$CDR2.IMGT.c.a +
                         dat$CDR2.IMGT.g.c +
                         dat$CDR2.IMGT.c.g +
                         dat$CDR2.IMGT.g.t +
                         dat$FR3.IMGT.g.a +
                         dat$FR3.IMGT.c.t +
                         dat$FR3.IMGT.c.a +
                         dat$FR3.IMGT.g.c +
                         dat$FR3.IMGT.c.g +
                         dat$FR3.IMGT.g.t



setwd(outputdir)

matrx = matrix(data = 0, ncol=((length(genes) + 1) * 3),nrow=5)
for(i in 1:length(genes)){
  gene = genes[i]
  tmp = dat[grepl(paste(".*", gene, ".*", sep=""), dat$best_match),]
  if(gene == "."){
    tmp = dat
  }
  if(length(tmp) == 0){
		cat("0", file=paste(gene, "_value.txt" ,sep=""))
    next
  }
  j = i - 1
  x = (j * 3) + 1
  y = (j * 3) + 2
  z = (j * 3) + 3
  matrx[1,x] = sum(tmp$VRegionMutations)
  matrx[1,y] = sum(tmp$VRegionNucleotides)
  matrx[1,z] = round(matrx[1,x] / matrx[1,y] * 100, digits=1)
  matrx[2,x] = sum(tmp$transitionMutations)
  matrx[2,y] = sum(tmp$VRegionMutations)
  matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1)
  matrx[3,x] = sum(tmp$transversionMutations)
  matrx[3,y] = sum(tmp$VRegionMutations)
  matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1)
  matrx[4,x] = sum(tmp$transitionMutationsAtGC)
  matrx[4,y] = sum(tmp$totalMutationsAtGC)
  matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1)
  matrx[5,x] = sum(tmp$totalMutationsAtGC)
  matrx[5,y] = sum(tmp$VRegionMutations)
  matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1)
  
  transitionTable = data.frame(A=1:4,C=1:4,G=1:4,T=1:4)
  row.names(transitionTable) = c("A", "C", "G", "T")
  transitionTable["A","A"] = NA
  transitionTable["C","C"] = NA
  transitionTable["G","G"] = NA
  transitionTable["T","T"] = NA
  nts = c("a", "c", "g", "t")
  
  
  for(nt1 in nts){
    for(nt2 in nts){
      if(nt1 == nt2){
        next
      }
      NT1 = LETTERS[letters == nt1]
      NT2 = LETTERS[letters == nt2]
      FR1 = 0 #paste("FR1.IMGT.", nt1, ".", nt2, sep="")
      CDR1 = paste("CDR1.IMGT.", nt1, ".", nt2, sep="")
      FR2 = paste("FR2.IMGT.", nt1, ".", nt2, sep="")
      CDR2 = paste("CDR2.IMGT.", nt1, ".", nt2, sep="")
      FR3 = paste("FR3.IMGT.", nt1, ".", nt2, sep="")
      transitionTable[NT1,NT2] = sum( tmp[,CDR1] +
                                      tmp[,FR2] +
                                      tmp[,CDR2] +
                                      tmp[,FR3])
    }
  }
  write.table(x=transitionTable, file=paste("transitions_", gene ,".txt", sep=""), sep=",",quote=F,row.names=T,col.names=NA)
  write.table(x=tmp[,c("Sequence.ID", "best_match", "chunk_hit_percentage", "nt_hit_percentage", "start_locations")], file=paste("matched_", gene ,".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T)
  
  cat(matrx[1,x], file=paste(gene, "_value.txt" ,sep=""))
  cat(length(tmp$Sequence.ID), file=paste(gene, "_n.txt" ,sep=""))
}

#again for all of the data
tmp = dat
j = i
x = (j * 3) + 1
y = (j * 3) + 2
z = (j * 3) + 3
matrx[1,x] = sum(tmp$VRegionMutations)
matrx[1,y] = sum(tmp$VRegionNucleotides)
matrx[1,z] = round(matrx[1,x] / matrx[1,y] * 100, digits=1)
matrx[2,x] = sum(tmp$transitionMutations)
matrx[2,y] = sum(tmp$VRegionMutations)
matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1)
matrx[3,x] = sum(tmp$transversionMutations)
matrx[3,y] = sum(tmp$VRegionMutations)
matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1)
matrx[4,x] = sum(tmp$transitionMutationsAtGC)
matrx[4,y] = sum(tmp$totalMutationsAtGC)
matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1)
matrx[5,x] = sum(tmp$totalMutationsAtGC)
matrx[5,y] = sum(tmp$VRegionMutations)
matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1)

transitionTable = data.frame(A=1:4,C=1:4,G=1:4,T=1:4)
row.names(transitionTable) = c("A", "C", "G", "T")
transitionTable["A","A"] = NA
transitionTable["C","C"] = NA
transitionTable["G","G"] = NA
transitionTable["T","T"] = NA
nts = c("a", "c", "g", "t")


for(nt1 in nts){
	for(nt2 in nts){
		if(nt1 == nt2){
			next
		}
		NT1 = LETTERS[letters == nt1]
		NT2 = LETTERS[letters == nt2]
		FR1 = paste("FR1.IMGT.", nt1, ".", nt2, sep="")
		CDR1 = paste("CDR1.IMGT.", nt1, ".", nt2, sep="")
		FR2 = paste("FR2.IMGT.", nt1, ".", nt2, sep="")
		CDR2 = paste("CDR2.IMGT.", nt1, ".", nt2, sep="")
		FR3 = paste("FR3.IMGT.", nt1, ".", nt2, sep="")
		transitionTable[NT1,NT2] = sum( tmp[,CDR1] +
																		tmp[,FR2] +
																		tmp[,CDR2] +
																		tmp[,FR3])
	}
}
write.table(x=transitionTable, file="transitions.txt", sep=",",quote=F,row.names=T,col.names=NA)
write.table(x=tmp[,c("Sequence.ID", "best_match", "chunk_hit_percentage", "nt_hit_percentage", "start_locations")], file="matched_all.txt", sep="\t",quote=F,row.names=F,col.names=T)
cat(matrx[1,x], file="total_value.txt")
cat(length(tmp$Sequence.ID), file="total_n.txt")



result = data.frame(matrx)
row.names(result) = c("Number of Mutations (%)", "Transition (%)", "Transversions (%)", "Transitions at G C (%)", "Targeting of C.G (%)")

write.table(x=result, file="mutations.txt", sep=",",quote=F,row.names=T,col.names=F)


if (!("ggplot2" %in% rownames(installed.packages()))) {
	install.packages("ggplot2", repos="http://cran.xl-mirror.nl/") 
}
library(ggplot2)

genesForPlot = gsub("[0-9]", "", dat$best_match)
genesForPlot = data.frame(table(genesForPlot))
colnames(genesForPlot) = c("Gene","Freq")
genesForPlot$label = paste(genesForPlot$Gene, "-", genesForPlot$Freq)

pc = ggplot(genesForPlot, aes(x = factor(1), y=Freq, fill=label))
pc = pc + geom_bar(width = 1, stat = "identity")
pc = pc + coord_polar(theta="y")
pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IgA", "( n =", sum(genesForPlot$Freq), ")"))

png(filename="all.png")
pc
dev.off()


#blegh
genesForPlot = dat[grepl("ca", dat$best_match),]$best_match
if(length(genesForPlot) > 0){
	genesForPlot = data.frame(table(genesForPlot))
	colnames(genesForPlot) = c("Gene","Freq")
	genesForPlot$label = paste(genesForPlot$Gene, "-", genesForPlot$Freq)

	pc = ggplot(genesForPlot, aes(x = factor(1), y=Freq, fill=label))
	pc = pc + geom_bar(width = 1, stat = "identity")
	pc = pc + coord_polar(theta="y")
	pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IgA", "( n =", sum(genesForPlot$Freq), ")"))


	png(filename="ca.png")
	print(pc)
	dev.off()
}

genesForPlot = dat[grepl("cg", dat$best_match),]$best_match
if(length(genesForPlot) > 0){
	genesForPlot = data.frame(table(genesForPlot))
	colnames(genesForPlot) = c("Gene","Freq")
	genesForPlot$label = paste(genesForPlot$Gene, "-", genesForPlot$Freq)

	pc = ggplot(genesForPlot, aes(x = factor(1), y=Freq, fill=label))
	pc = pc + geom_bar(width = 1, stat = "identity")
	pc = pc + coord_polar(theta="y")
	pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IgG", "( n =", sum(genesForPlot$Freq), ")"))


	png(filename="cg.png")
	print(pc)
	dev.off()
}