Mercurial > repos > davidvanzessen > shm_csr
changeset 1:faae21ba5c63 draft
Uploaded
| author | davidvanzessen |
|---|---|
| date | Tue, 25 Oct 2016 07:28:43 -0400 |
| parents | c33d93683a09 |
| children | e85fec274cde |
| files | aa_histogram.r merge_and_filter.r pattern_plots.r shm_csr.py shm_csr.r shm_csr.xml wrapper.sh |
| diffstat | 7 files changed, 591 insertions(+), 567 deletions(-) [+] |
line wrap: on
line diff
--- a/aa_histogram.r Thu Oct 13 10:52:24 2016 -0400 +++ b/aa_histogram.r Tue Oct 25 07:28:43 2016 -0400 @@ -37,7 +37,7 @@ print("---------------- plot ----------------") - m = ggplot(dat_dt, aes(x=i, y=freq)) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + m = ggplot(dat_dt, aes(x=i, y=freq)) + theme(axis.text.x = element_text(angle = 90, hjust = 1), text = element_text(size=13, colour="black")) m = m + geom_bar(stat="identity", colour = "black", fill = "darkgrey", alpha=0.8) + scale_x_continuous(breaks=dat_dt$i, labels=dat_dt$i) m = m + annotate("segment", x = 0.5, y = -0.05, xend=26.5, yend=-0.05, colour="darkgreen", size=1) + annotate("text", x = 13, y = -0.1, label="FR1") m = m + annotate("segment", x = 26.5, y = -0.07, xend=38.5, yend=-0.07, colour="darkblue", size=1) + annotate("text", x = 32.5, y = -0.15, label="CDR1")
--- a/merge_and_filter.r Thu Oct 13 10:52:24 2016 -0400 +++ b/merge_and_filter.r Tue Oct 25 07:28:43 2016 -0400 @@ -32,7 +32,6 @@ gene_identification$chunk_hit_percentage = (gene_identification$length / gene_identification$ref.length) * 100 gene_identification = gene_identification[,c("qseqid", "chunk_hit_percentage", "pident", "qstart", "sseqid")] colnames(gene_identification) = c("Sequence.ID", "chunk_hit_percentage", "nt_hit_percentage", "start_locations", "best_match") - } input.sequence.count = nrow(summ) @@ -62,32 +61,9 @@ summ = summ[summ$Functionality != "No results" & summ$Functionality != "unknown (see comment)" & summ$Functionality != "unknown",] } -print(paste("Number of sequences after productive filter:", nrow(summ))) - -filtering.steps = rbind(filtering.steps, c("After productive filter", nrow(summ))) - -splt = strsplit(class.filter, "_")[[1]] -chunk_hit_threshold = as.numeric(splt[1]) -nt_hit_threshold = as.numeric(splt[2]) - -higher_than=(summ$chunk_hit_percentage >= chunk_hit_threshold & summ$nt_hit_percentage >= nt_hit_threshold) - -unmatched=summ[NULL,c("Sequence.ID", "chunk_hit_percentage", "nt_hit_percentage", "start_locations", "best_match")] +print(paste("Number of sequences after functionality filter:", nrow(summ))) -if(!all(higher_than, na.rm=T)){ #check for 'not all' because that would mean the unmatched set is empty - unmatched = summ[!higher_than,] - unmatched = unmatched[,c("Sequence.ID", "chunk_hit_percentage", "nt_hit_percentage", "start_locations", "best_match")] - unmatched$best_match = paste("unmatched,", unmatched$best_match) - summ[!higher_than,"best_match"] = paste("unmatched,", summ[!higher_than,"best_match"]) -} - -if(any(higher_than, na.rm=T)){ - #summ = summ[higher_than,] -} - -if(nrow(summ) == 0){ - stop("No data remaining after filter") -} +filtering.steps = rbind(filtering.steps, c("After functionality filter", nrow(summ))) result = merge(summ, mutationanalysis[,!(names(mutationanalysis) %in% names(summ)[-1])], by="Sequence.ID") @@ -114,22 +90,16 @@ result$JGene = gsub("^Homsap ", "", result$J.GENE.and.allele) result$JGene = gsub("[*].*", "", result$JGene) -result$past = do.call(paste, c(result[unlist(strsplit(unique.type, ","))], sep = ":")) - -result = result[!(duplicated(result$past)), ] - -result = result[,!(names(result) %in% c("past"))] - -print(paste("Number of sequences in result after", unique.type, "filtering:", nrow(result))) - -filtering.steps = rbind(filtering.steps, c("After duplicate filter", nrow(result))) - print(paste("Number of empty CDR1 sequences:", sum(result$CDR1.IMGT.seq == ""))) print(paste("Number of empty FR2 sequences:", sum(result$FR2.IMGT.seq == ""))) print(paste("Number of empty CDR2 sequences:", sum(result$CDR2.IMGT.seq == ""))) print(paste("Number of empty FR3 sequences:", sum(result$FR3.IMGT.seq == ""))) -if(empty.region.filter == "FR1"){ +if(empty.region.filter == "leader"){ + result = result[result$FR1.IMGT.seq != "" & result$CDR1.IMGT.seq != "" & result$FR2.IMGT.seq != "" & result$CDR2.IMGT.seq != "" & result$FR3.IMGT.seq != "", ] + print(paste("Number of sequences after empty FR1, CDR1, FR2, CDR2 and FR3 column filter:", nrow(result))) + filtering.steps = rbind(filtering.steps, c("After empty FR1, CDR1, FR2, CDR2, FR3 filter", nrow(result))) +} else if(empty.region.filter == "FR1"){ result = result[result$CDR1.IMGT.seq != "" & result$FR2.IMGT.seq != "" & result$CDR2.IMGT.seq != "" & result$FR3.IMGT.seq != "", ] print(paste("Number of sequences after empty CDR1, FR2, CDR2 and FR3 column filter:", nrow(result))) filtering.steps = rbind(filtering.steps, c("After empty CDR1, FR2, CDR2, FR3 filter", nrow(result))) @@ -143,7 +113,12 @@ filtering.steps = rbind(filtering.steps, c("After empty CDR2, FR3 filter", nrow(result))) } -if(empty.region.filter == "FR1"){ +print(paste("Number of sequences in result after CDR/FR filtering:", nrow(result))) +print(paste("Number of matched sequences in result after CDR/FR filtering:", nrow(result[!grepl("unmatched", result$best_match),]))) + +if(empty.region.filter == "leader"){ + result = result[!(grepl("n|N", result$FR1.IMGT.seq) | grepl("n|N", result$FR2.IMGT.seq) | grepl("n|N", result$FR3.IMGT.seq) | grepl("n|N", result$CDR1.IMGT.seq) | grepl("n|N", result$CDR2.IMGT.seq) | grepl("n|N", result$CDR3.IMGT.seq)),] +} else if(empty.region.filter == "FR1"){ result = result[!(grepl("n|N", result$FR2.IMGT.seq) | grepl("n|N", result$FR3.IMGT.seq) | grepl("n|N", result$CDR1.IMGT.seq) | grepl("n|N", result$CDR2.IMGT.seq) | grepl("n|N", result$CDR3.IMGT.seq)),] } else if(empty.region.filter == "CDR1"){ result = result[!(grepl("n|N", result$FR2.IMGT.seq) | grepl("n|N", result$FR3.IMGT.seq) | grepl("n|N", result$CDR2.IMGT.seq) | grepl("n|N", result$CDR3.IMGT.seq)),] @@ -171,7 +146,9 @@ if(filter.unique != "no"){ clmns = names(result) - if(empty.region.filter == "FR1"){ + if(empty.region.filter == "leader"){ + result$unique.def = paste(result$FR1.IMGT.seq, result$CDR1.IMGT.seq, result$FR2.IMGT.seq, result$CDR2.IMGT.seq, result$FR3.IMGT.seq, result$CDR3.IMGT.seq) + } else if(empty.region.filter == "FR1"){ result$unique.def = paste(result$CDR1.IMGT.seq, result$FR2.IMGT.seq, result$CDR2.IMGT.seq, result$FR3.IMGT.seq, result$CDR3.IMGT.seq) } else if(empty.region.filter == "CDR1"){ rresult$unique.def = paste(result$FR2.IMGT.seq, result$CDR2.IMGT.seq, result$FR3.IMGT.seq, result$CDR3.IMGT.seq) @@ -199,10 +176,41 @@ #write.table(inputdata.removed, "unique_removed.csv", sep=",",quote=F,row.names=F,col.names=T) } -print(paste("Number of sequences in result after CDR/FR filtering:", nrow(result))) -print(paste("Number of matched sequences in result after CDR/FR filtering:", nrow(result[!grepl("unmatched", result$best_match),]))) +filtering.steps = rbind(filtering.steps, c("After filter unique sequences", nrow(result))) + + +splt = strsplit(class.filter, "_")[[1]] +chunk_hit_threshold = as.numeric(splt[1]) +nt_hit_threshold = as.numeric(splt[2]) + +higher_than=(summ$chunk_hit_percentage >= chunk_hit_threshold & summ$nt_hit_percentage >= nt_hit_threshold) + +unmatched=summ[NULL,c("Sequence.ID", "chunk_hit_percentage", "nt_hit_percentage", "start_locations", "best_match")] + +if(!all(higher_than, na.rm=T)){ #check for 'not all' because that would mean the unmatched set is empty + unmatched = summ[!higher_than,] + unmatched = unmatched[,c("Sequence.ID", "chunk_hit_percentage", "nt_hit_percentage", "start_locations", "best_match")] + unmatched$best_match = paste("unmatched,", unmatched$best_match) + summ[!higher_than,"best_match"] = paste("unmatched,", summ[!higher_than,"best_match"]) +} -filtering.steps = rbind(filtering.steps, c("After unique filter", nrow(result))) +if(any(higher_than, na.rm=T)){ + #summ = summ[higher_than,] +} + +if(nrow(summ) == 0){ + stop("No data remaining after filter") +} + +result$past = do.call(paste, c(result[unlist(strsplit(unique.type, ","))], sep = ":")) + +result = result[!(duplicated(result$past)), ] + +result = result[,!(names(result) %in% c("past"))] + +print(paste("Number of sequences in result after", unique.type, "filtering:", nrow(result))) + +filtering.steps = rbind(filtering.steps, c("After remove duplicates based on filter", nrow(result))) print(paste("Number of rows in result:", nrow(result))) print(paste("Number of rows in unmatched:", nrow(unmatched)))
--- a/pattern_plots.r Thu Oct 13 10:52:24 2016 -0400 +++ b/pattern_plots.r Tue Oct 25 07:28:43 2016 -0400 @@ -45,7 +45,7 @@ write.table(data1, plot1.txt, quote=F, sep="\t", na="", row.names=F, col.names=T) p = ggplot(data1, aes(Class, value)) + geom_bar(aes(fill=Type), stat="identity", position="dodge", colour = "black") + ylab("% of mutations") + guides(fill=guide_legend(title=NULL)) -p = p + theme(panel.background = element_rect(fill = "white", colour="black")) + scale_fill_manual(values=c("RGYW.WRCY" = "white", "TW.WA" = "blue4")) +p = p + theme(panel.background = element_rect(fill = "white", colour="black"),text = element_text(size=13, colour="black")) + scale_fill_manual(values=c("RGYW.WRCY" = "white", "TW.WA" = "blue4")) #p = p + scale_colour_manual(values=c("RGYW.WRCY" = "black", "TW.WA" = "blue4")) png(filename=plot1.png, width=480, height=300) print(p) @@ -79,7 +79,7 @@ write.table(data2, plot2.txt, quote=F, sep="\t", na="", row.names=F, col.names=T) p = ggplot(data2, aes(x=Class, y=value, fill=Type)) + geom_bar(position="fill", stat="identity", colour = "black") + scale_y_continuous(labels=percent_format()) + guides(fill=guide_legend(title=NULL)) + ylab("% of mutations") -p = p + theme(panel.background = element_rect(fill = "white", colour="black")) + scale_fill_manual(values=c("A/T" = "blue4", "G/C transversions" = "gray74", "G/C transitions" = "white")) +p = p + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=13, colour="black")) + scale_fill_manual(values=c("A/T" = "blue4", "G/C transversions" = "gray74", "G/C transitions" = "white")) #p = p + scale_colour_manual(values=c("A/T" = "blue4", "G/C transversions" = "gray74", "G/C transitions" = "black")) png(filename=plot2.png, width=480, height=300) print(p) @@ -115,7 +115,7 @@ write.table(data3, plot3.txt, quote=F, sep="\t", na="", row.names=F, col.names=T) p = ggplot(data3, aes(Class, value)) + geom_bar(aes(fill=Type), stat="identity", position="dodge", colour = "black") + ylab("% of nucleotides") + guides(fill=guide_legend(title=NULL)) -p = p + theme(panel.background = element_rect(fill = "white", colour="black")) + scale_fill_manual(values=c("A/T" = "blue4", "G/C transversions" = "gray74", "G/C transitions" = "white")) +p = p + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=13, colour="black")) + scale_fill_manual(values=c("A/T" = "blue4", "G/C transversions" = "gray74", "G/C transitions" = "white")) #p = p + scale_colour_manual(values=c("A/T" = "blue4", "G/C transversions" = "gray74", "G/C transitions" = "black")) png(filename=plot3.png, width=480, height=300) print(p)
--- a/shm_csr.py Thu Oct 13 10:52:24 2016 -0400 +++ b/shm_csr.py Tue Oct 25 07:28:43 2016 -0400 @@ -7,15 +7,14 @@ parser.add_argument("--input", help="The '7_V-REGION-mutation-and-AA-change-table' and '10_V-REGION-mutation-hotspots' merged together, with an added 'best_match' annotation") parser.add_argument("--genes", help="The genes available in the 'best_match' column") -parser.add_argument("--includefr1", help="Should the mutation/nucleotides in the FR1 region be included?") +parser.add_argument("--empty_region_filter", help="Where does the sequence start?", choices=['leader', 'FR1', 'CDR1', 'FR2']) parser.add_argument("--output", help="Output file") args = parser.parse_args() infile = args.input genes = str(args.genes).split(",") -print "includefr1 =", args.includefr1 -include_fr1 = True if args.includefr1 == "yes" else False +empty_region_filter = args.empty_region_filter outfile = args.output genedic = dict() @@ -59,16 +58,14 @@ ID = linesplt[IDIndex] genedic[ID] = linesplt[best_matchIndex] try: - if linesplt[fr1Index] != "NA": - mutationdic[ID + "_FR1"] = [mutationMatcher.match(x).groups() for x in linesplt[fr1Index].split("|") if x] if include_fr1 else [] - else: - mutationdic[ID + "_FR1"] = [] - mutationdic[ID + "_CDR1"] = [mutationMatcher.match(x).groups() for x in linesplt[cdr1Index].split("|") if x] if linesplt[cdr1Index] != "NA" else [] - mutationdic[ID + "_FR2"] = [mutationMatcher.match(x).groups() for x in linesplt[fr2Index].split("|") if x] if linesplt[fr2Index] != "NA" else [] - mutationdic[ID + "_CDR2"] = [mutationMatcher.match(x).groups() for x in linesplt[cdr2Index].split("|") if x] if linesplt[cdr2Index] != "NA" else [] + mutationdic[ID + "_FR1"] = [mutationMatcher.match(x).groups() for x in linesplt[fr1Index].split("|") if x] if (linesplt[fr1Index] != "NA" and empty_region_filter == "leader") else [] + mutationdic[ID + "_CDR1"] = [mutationMatcher.match(x).groups() for x in linesplt[cdr1Index].split("|") if x] if (linesplt[cdr1Index] != "NA" and empty_region_filter in ["leader", "FR1"]) else [] + mutationdic[ID + "_FR2"] = [mutationMatcher.match(x).groups() for x in linesplt[fr2Index].split("|") if x] if (linesplt[fr2Index] != "NA" and empty_region_filter in ["leader", "FR1", "CDR1"]) else [] + mutationdic[ID + "_CDR2"] = [mutationMatcher.match(x).groups() for x in linesplt[cdr2Index].split("|") if x] if (linesplt[cdr2Index] != "NA") else [] mutationdic[ID + "_FR2-CDR2"] = mutationdic[ID + "_FR2"] + mutationdic[ID + "_CDR2"] mutationdic[ID + "_FR3"] = [mutationMatcher.match(x).groups() for x in linesplt[fr3Index].split("|") if x] if linesplt[fr3Index] != "NA" else [] - except e: + except Exception as e: + print "Something went wrong while processing this line:" print linesplt print linecount print e @@ -189,8 +186,6 @@ linesplt = line.split("\t") gene = linesplt[best_matchIndex] ID = linesplt[IDIndex] - if ID == "ca2": - print linesplt RGYW = [(int(x), int(y), z) for (x, y, z) in [hotspotMatcher.match(x).groups() for x in linesplt[aggctatIndex].split("|") if x]] WRCY = [(int(x), int(y), z) for (x, y, z) in @@ -201,8 +196,7 @@ [hotspotMatcher.match(x).groups() for x in linesplt[tatIndex].split("|") if x]] RGYWCount[ID], WRCYCount[ID], WACount[ID], TWCount[ID] = 0, 0, 0, 0 - mutationList = (mutationdic[ID + "_FR1"] if include_fr1 else []) + mutationdic[ID + "_CDR1"] + mutationdic[ - ID + "_FR2"] + mutationdic[ID + "_CDR2"] + mutationdic[ID + "_FR3"] + mutationList = mutationdic[ID + "_FR1"] + mutationdic[ID + "_CDR1"] + mutationdic[ID + "_FR2"] + mutationdic[ID + "_CDR2"] + mutationdic[ID + "_FR3"] for mutation in mutationList: frm, where, to, AAfrm, AAwhere, AAto, junk = mutation mutation_in_RGYW = any([(start <= int(where) <= end) for (start, end, region) in RGYW]) @@ -271,7 +265,7 @@ o.write(typ + " (%)") curr = dic[typ] for gene in genes: - geneMatcher = geneMatchers[gene] #re.compile("^" + gene + ".*") #recompile every loop.... + geneMatcher = geneMatchers[gene] if valuedic[gene + "_" + fname] is 0: o.write(",0,0,0") else:
--- a/shm_csr.r Thu Oct 13 10:52:24 2016 -0400 +++ b/shm_csr.r Tue Oct 25 07:28:43 2016 -0400 @@ -1,493 +1,516 @@ -library(data.table) -library(ggplot2) -library(reshape2) - -args <- commandArgs(trailingOnly = TRUE) - -input = args[1] -genes = unlist(strsplit(args[2], ",")) -outputdir = args[3] -include_fr1 = ifelse(args[4] == "yes", T, F) -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", - "CDR3.IMGT.Nb.of.nucleotides", - "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", - "FR1.IMGT.Nb.of.silent.mutations", - "FR2.IMGT.Nb.of.silent.mutations", - "FR3.IMGT.Nb.of.silent.mutations", - "FR1.IMGT.Nb.of.nonsilent.mutations", - "FR2.IMGT.Nb.of.nonsilent.mutations", - "FR3.IMGT.Nb.of.nonsilent.mutations") - - -print("Cleaning up columns") -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]),col] = 0 -} - -regions = c("FR1", "CDR1", "FR2", "CDR2", "FR3") -if(!include_fr1){ - regions = c("CDR1", "FR2", "CDR2", "FR3") -} - -sum_by_row = function(x, columns) { sum(as.numeric(x[columns]), na.rm=T) } - -print("aggregating data into new columns") - -VRegionMutations_columns = paste(regions, ".IMGT.Nb.of.mutations", sep="") -dat$VRegionMutations = apply(dat, FUN=sum_by_row, 1, columns=VRegionMutations_columns) - -VRegionNucleotides_columns = paste(regions, ".IMGT.Nb.of.nucleotides", sep="") -dat$FR3.IMGT.Nb.of.nucleotides = nchar(dat$FR3.IMGT.seq) -dat$VRegionNucleotides = apply(dat, FUN=sum_by_row, 1, columns=VRegionNucleotides_columns) - -transitionMutations_columns = paste(rep(regions, each=4), c(".IMGT.a.g", ".IMGT.g.a", ".IMGT.c.t", ".IMGT.t.c"), sep="") -dat$transitionMutations = apply(dat, FUN=sum_by_row, 1, columns=transitionMutations_columns) - -transversionMutations_columns = paste(rep(regions, each=8), c(".IMGT.a.c",".IMGT.c.a",".IMGT.a.t",".IMGT.t.a",".IMGT.g.c",".IMGT.c.g",".IMGT.g.t",".IMGT.t.g"), sep="") -dat$transversionMutations = apply(dat, FUN=sum_by_row, 1, columns=transversionMutations_columns) - - -transitionMutationsAtGC_columns = paste(rep(regions, each=2), c(".IMGT.g.a",".IMGT.c.t"), sep="") -dat$transitionMutationsAtGC = apply(dat, FUN=sum_by_row, 1, columns=transitionMutationsAtGC_columns) - - -totalMutationsAtGC_columns = paste(rep(regions, each=6), c(".IMGT.c.g",".IMGT.c.t",".IMGT.c.a",".IMGT.g.c",".IMGT.g.a",".IMGT.g.t"), sep="") -#totalMutationsAtGC_columns = paste(rep(regions, each=6), c(".IMGT.g.a",".IMGT.c.t",".IMGT.c.a",".IMGT.c.g",".IMGT.g.t"), sep="") -dat$totalMutationsAtGC = apply(dat, FUN=sum_by_row, 1, columns=totalMutationsAtGC_columns) - -transitionMutationsAtAT_columns = paste(rep(regions, each=2), c(".IMGT.a.g",".IMGT.t.c"), sep="") -dat$transitionMutationsAtAT = apply(dat, FUN=sum_by_row, 1, columns=transitionMutationsAtAT_columns) - -totalMutationsAtAT_columns = paste(rep(regions, each=6), c(".IMGT.a.g",".IMGT.a.c",".IMGT.a.t",".IMGT.t.g",".IMGT.t.c",".IMGT.t.a"), sep="") -#totalMutationsAtAT_columns = paste(rep(regions, each=5), c(".IMGT.a.g",".IMGT.t.c",".IMGT.a.c",".IMGT.g.c",".IMGT.t.g"), sep="") -dat$totalMutationsAtAT = apply(dat, FUN=sum_by_row, 1, columns=totalMutationsAtAT_columns) - - -FRRegions = regions[grepl("FR", regions)] -CDRRegions = regions[grepl("CDR", regions)] - -FR_silentMutations_columns = paste(FRRegions, ".IMGT.Nb.of.silent.mutations", sep="") -dat$silentMutationsFR = apply(dat, FUN=sum_by_row, 1, columns=FR_silentMutations_columns) - -CDR_silentMutations_columns = paste(CDRRegions, ".IMGT.Nb.of.silent.mutations", sep="") -dat$silentMutationsCDR = apply(dat, FUN=sum_by_row, 1, columns=CDR_silentMutations_columns) - -FR_nonSilentMutations_columns = paste(FRRegions, ".IMGT.Nb.of.nonsilent.mutations", sep="") -dat$nonSilentMutationsFR = apply(dat, FUN=sum_by_row, 1, columns=FR_nonSilentMutations_columns) - -CDR_nonSilentMutations_columns = paste(CDRRegions, ".IMGT.Nb.of.nonsilent.mutations", sep="") -dat$nonSilentMutationsCDR = apply(dat, FUN=sum_by_row, 1, columns=CDR_nonSilentMutations_columns) - -mutation.sum.columns = c("Sequence.ID", "VRegionMutations", "VRegionNucleotides", "transitionMutations", "transversionMutations", "transitionMutationsAtGC", "transitionMutationsAtAT", "silentMutationsFR", "nonSilentMutationsFR", "silentMutationsCDR", "nonSilentMutationsCDR") - -write.table(dat[,mutation.sum.columns], "mutation_by_id.txt", sep="\t",quote=F,row.names=F,col.names=T) - -setwd(outputdir) - -base.order = data.frame(base=c("A", "T", "C", "G"), order=1:4) - -calculate_result = function(i, gene, dat, matrx, f, fname, name){ - tmp = dat[grepl(paste("^", gene, ".*", sep=""), dat$best_match),] - - j = i - 1 - x = (j * 3) + 1 - y = (j * 3) + 2 - z = (j * 3) + 3 - - if(nrow(tmp) > 0){ - - if(fname == "sum"){ - matrx[1,x] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) - matrx[1,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1) - matrx[1,z] = round(f(matrx[1,x] / matrx[1,y]) * 100, digits=1) - } else { - matrx[1,x] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) - matrx[1,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1) - matrx[1,z] = round(f(tmp$VRegionMutations / tmp$VRegionNucleotides) * 100, digits=1) - } - - matrx[2,x] = round(f(tmp$transitionMutations, na.rm=T), digits=1) - matrx[2,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) - matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1) - - matrx[3,x] = round(f(tmp$transversionMutations, na.rm=T), digits=1) - matrx[3,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) - matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1) - - matrx[4,x] = round(f(tmp$transitionMutationsAtGC, na.rm=T), digits=1) - matrx[4,y] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1) - matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1) - - matrx[5,x] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1) - matrx[5,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) - matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1) - - matrx[6,x] = round(f(tmp$transitionMutationsAtAT, na.rm=T), digits=1) - matrx[6,y] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1) - matrx[6,z] = round(matrx[6,x] / matrx[6,y] * 100, digits=1) - - matrx[7,x] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1) - matrx[7,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) - matrx[7,z] = round(matrx[7,x] / matrx[7,y] * 100, digits=1) - - matrx[8,x] = round(f(tmp$nonSilentMutationsFR, na.rm=T), digits=1) - matrx[8,y] = round(f(tmp$silentMutationsFR, na.rm=T), digits=1) - matrx[8,z] = round(matrx[8,x] / matrx[8,y], digits=1) - - matrx[9,x] = round(f(tmp$nonSilentMutationsCDR, na.rm=T), digits=1) - matrx[9,y] = round(f(tmp$silentMutationsCDR, na.rm=T), digits=1) - matrx[9,z] = round(matrx[9,x] / matrx[9,y], digits=1) - - if(fname == "sum"){ - matrx[10,x] = round(f(rowSums(tmp[,c("FR2.IMGT.Nb.of.nucleotides", "FR3.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1) - matrx[10,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1) - matrx[10,z] = round(matrx[10,x] / matrx[10,y] * 100, digits=1) - - matrx[11,x] = round(f(rowSums(tmp[,c("CDR1.IMGT.Nb.of.nucleotides", "CDR2.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1) - matrx[11,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1) - matrx[11,z] = round(matrx[11,x] / matrx[11,y] * 100, digits=1) - } - } - - transitionTable = data.frame(A=zeros,C=zeros,G=zeros,T=zeros) - row.names(transitionTable) = c("A", "C", "G", "T") - transitionTable["A","A"] = NA - transitionTable["C","C"] = NA - transitionTable["G","G"] = NA - transitionTable["T","T"] = NA - - if(nrow(tmp) > 0){ - 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="") - if(include_fr1){ - transitionTable[NT1,NT2] = sum(tmp[,c(FR1, CDR1, FR2, CDR2, FR3)]) - } else { - transitionTable[NT1,NT2] = sum(tmp[,c(CDR1, FR2, CDR2, FR3)]) - } - } - } - transition = transitionTable - transition$id = names(transition) - - transition2 = melt(transition, id.vars="id") - - transition2 = merge(transition2, base.order, by.x="id", by.y="base") - transition2 = merge(transition2, base.order, by.x="variable", by.y="base") - - transition2[is.na(transition2$value),]$value = 0 - - if(!all(transition2$value == 0)){ #having rows of data but a transition table filled with 0 is bad - - print("Plotting stacked transition") - - png(filename=paste("transitions_stacked_", name, ".png", sep="")) - p = ggplot(transition2, aes(factor(reorder(id, order.x)), y=value, fill=factor(reorder(variable, order.y)))) + geom_bar(position="fill", stat="identity", colour="black") #stacked bar - p = p + xlab("From base") + ylab("To base") + ggtitle("Mutations frequency from base to base") + guides(fill=guide_legend(title=NULL)) - p = p + theme(panel.background = element_rect(fill = "white", colour="black")) + scale_fill_manual(values=c("A" = "blue4", "G" = "lightblue1", "C" = "olivedrab3", "T" = "olivedrab4")) - #p = p + scale_colour_manual(values=c("A" = "black", "G" = "black", "C" = "black", "T" = "black")) - print(p) - dev.off() - - print("Plotting heatmap transition") - - png(filename=paste("transitions_heatmap_", name, ".png", sep="")) - p = ggplot(transition2, aes(factor(reorder(id, order.x)), factor(reorder(variable, order.y)))) + geom_tile(aes(fill = value)) + scale_fill_gradient(low="white", high="steelblue") #heatmap - p = p + xlab("From base") + ylab("To base") + ggtitle("Mutations frequency from base to base") + theme(panel.background = element_rect(fill = "white", colour="black")) - print(p) - dev.off() - } else { - print("No data to plot") - } - } - - #print(paste("writing value file: ", name, "_", fname, "_value.txt" ,sep="")) - - write.table(x=transitionTable, file=paste("transitions_", name ,"_", fname, ".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_", name , "_", fname, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T) - - cat(matrx[1,x], file=paste(name, "_", fname, "_value.txt" ,sep="")) - cat(nrow(tmp), file=paste(name, "_", fname, "_n.txt" ,sep="")) - - #print(paste(fname, name, nrow(tmp))) - - matrx -} - -nts = c("a", "c", "g", "t") -zeros=rep(0, 4) - -funcs = c(median, sum, mean) -fnames = c("median", "sum", "mean") - -print("Creating result tables") - -for(i in 1:length(funcs)){ - func = funcs[[i]] - fname = fnames[[i]] - - rows = 9 - if(fname == "sum"){ - rows = 11 - } - matrx = matrix(data = 0, ncol=((length(genes) + 1) * 3),nrow=rows) - - for(i in 1:length(genes)){ - print(paste("Creating table for", fname, genes[i])) - matrx = calculate_result(i, genes[i], dat, matrx, func, fname, genes[i]) - } - - matrx = calculate_result(i + 1, ".*", dat[!grepl("unmatched", dat$best_match),], matrx, func, fname, name="all") - - result = data.frame(matrx) - if(fname == "sum"){ - row.names(result) = c("Number of Mutations (%)", "Transitions (%)", "Transversions (%)", "Transitions at G C (%)", "Targeting of C G (%)", "Transitions at A T (%)", "Targeting of A T (%)", "FR R/S (ratio)", "CDR R/S (ratio)", "nt in FR", "nt in CDR") - } else { - row.names(result) = c("Number of Mutations (%)", "Transitions (%)", "Transversions (%)", "Transitions at G C (%)", "Targeting of C G (%)", "Transitions at A T (%)", "Targeting of A T (%)", "FR R/S (ratio)", "CDR R/S (ratio)") - } - - write.table(x=result, file=paste("mutations_", fname, ".txt", sep=""), sep=",",quote=F,row.names=T,col.names=F) -} - -print("Adding median number of mutations to sum table") - -sum.table = read.table("mutations_sum.txt", sep=",", header=F) -median.table = read.table("mutations_median.txt", sep=",", header=F) - -new.table = sum.table[1,] -new.table[2,] = median.table[1,] -new.table[3:12,] = sum.table[2:11,] -new.table[,1] = as.character(new.table[,1]) -new.table[2,1] = "Median of Number of Mutations (%)" - -#sum.table = sum.table[c("Number of Mutations (%)", "Median of Number of Mutations (%)", "Transition (%)", "Transversions (%)", "Transitions at G C (%)", "Targeting of C G (%)", "Transitions at A T (%)", "Targeting of A T (%)", "FR R/S (ratio)", "CDR R/S (ratio)", "nt in FR", "nt in CDR"),] - -write.table(x=new.table, file="mutations_sum.txt", sep=",",quote=F,row.names=F,col.names=F) - - -print("Plotting IGA piechart") - -dat = dat[!grepl("^unmatched", dat$best_match),] - -#blegh -genesForPlot = dat[grepl("IGA", 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=Gene)) - pc = pc + geom_bar(width = 1, stat = "identity") + scale_fill_manual(labels=genesForPlot$label, values=c("IGA1" = "lightblue1", "IGA2" = "blue4")) - pc = pc + coord_polar(theta="y") - pc = pc + theme(panel.background = element_rect(fill = "white", colour="black")) - pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IGA subclasses", "( n =", sum(genesForPlot$Freq), ")")) - write.table(genesForPlot, "IGA.txt", sep="\t",quote=F,row.names=F,col.names=T) - - png(filename="IGA.png") - print(pc) - dev.off() -} - -print("Plotting IGG piechart") - -genesForPlot = dat[grepl("IGG", 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=Gene)) - pc = pc + geom_bar(width = 1, stat = "identity") + scale_fill_manual(labels=genesForPlot$label, values=c("IGG1" = "olivedrab3", "IGG2" = "red", "IGG3" = "gold", "IGG4" = "darkred")) - pc = pc + coord_polar(theta="y") - pc = pc + theme(panel.background = element_rect(fill = "white", colour="black")) - pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IGG subclasses", "( n =", sum(genesForPlot$Freq), ")")) - write.table(genesForPlot, "IGG.txt", sep="\t",quote=F,row.names=F,col.names=T) - - png(filename="IGG.png") - print(pc) - dev.off() -} - - -print("Plotting scatterplot") - -dat$percentage_mutations = round(dat$VRegionMutations / dat$VRegionNucleotides * 100, 2) - -p = ggplot(dat, aes(best_match, percentage_mutations)) -p = p + geom_point(aes(colour=best_match), position="jitter") + geom_boxplot(aes(middle=mean(percentage_mutations)), alpha=0.1, outlier.shape = NA) -p = p + xlab("Subclass") + ylab("Frequency") + ggtitle("Frequency scatter plot") + theme(panel.background = element_rect(fill = "white", colour="black")) -p = p + scale_fill_manual(values=c("IGA1" = "lightblue1", "IGA2" = "blue4", "IGG1" = "olivedrab3", "IGG2" = "red", "IGG3" = "gold", "IGG4" = "darkred", "IGM" = "black")) -p = p + scale_colour_manual(values=c("IGA1" = "lightblue1", "IGA2" = "blue4", "IGG1" = "olivedrab3", "IGG2" = "red", "IGG3" = "gold", "IGG4" = "darkred", "IGM" = "black")) - -png(filename="scatter.png") -print(p) -dev.off() - -write.table(dat[,c("Sequence.ID", "best_match", "VRegionMutations", "VRegionNucleotides", "percentage_mutations")], "scatter.txt", sep="\t",quote=F,row.names=F,col.names=T) - -write.table(dat, input, sep="\t",quote=F,row.names=F,col.names=T) - - -print("Plotting frequency ranges plot") - -dat$best_match_class = substr(dat$best_match, 0, 3) -freq_labels = c("0", "0-2", "2-5", "5-10", "10-15", "15-20", "20") -dat$frequency_bins = cut(dat$percentage_mutations, breaks=c(-Inf, 0, 2,5,10,15,20, Inf), labels=freq_labels) - -frequency_bins_sum = data.frame(data.table(dat)[, list(class_sum=sum(.N)), by=c("best_match_class")]) - -frequency_bins_data = data.frame(data.table(dat)[, list(frequency_count=.N), by=c("best_match_class", "frequency_bins")]) - -frequency_bins_data = merge(frequency_bins_data, frequency_bins_sum, by="best_match_class") - -frequency_bins_data$frequency = round(frequency_bins_data$frequency_count / frequency_bins_data$class_sum * 100, 2) - -p = ggplot(frequency_bins_data, aes(frequency_bins, frequency)) -p = p + geom_bar(aes(fill=best_match_class), stat="identity", position="dodge") + theme(panel.background = element_rect(fill = "white", colour="black")) -p = p + xlab("Frequency ranges") + ylab("Frequency") + ggtitle("Mutation Frequencies by class") + scale_fill_manual(values=c("IGA" = "blue4", "IGG" = "olivedrab3", "IGM" = "black")) - -png(filename="frequency_ranges.png") -print(p) -dev.off() - -frequency_bins_data_by_class = frequency_bins_data - -write.table(frequency_bins_data_by_class, "frequency_ranges_classes.txt", sep="\t",quote=F,row.names=F,col.names=T) - -frequency_bins_data = data.frame(data.table(dat)[, list(frequency_count=.N), by=c("best_match", "best_match_class", "frequency_bins")]) - -frequency_bins_data = merge(frequency_bins_data, frequency_bins_sum, by="best_match_class") - -frequency_bins_data$frequency = round(frequency_bins_data$frequency_count / frequency_bins_data$class_sum * 100, 2) - -write.table(frequency_bins_data, "frequency_ranges_subclasses.txt", sep="\t",quote=F,row.names=F,col.names=T) - - -#frequency_bins_data_by_class -#frequency_ranges_subclasses.txt - - - - - - - - - - - - - - - - - - - - - - - - - - - +library(data.table) +library(ggplot2) +library(reshape2) + +args <- commandArgs(trailingOnly = TRUE) + +input = args[1] +genes = unlist(strsplit(args[2], ",")) +outputdir = args[3] +empty.region.filter = args[4] +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", + "CDR3.IMGT.Nb.of.nucleotides", + "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", + "FR1.IMGT.Nb.of.silent.mutations", + "FR2.IMGT.Nb.of.silent.mutations", + "FR3.IMGT.Nb.of.silent.mutations", + "FR1.IMGT.Nb.of.nonsilent.mutations", + "FR2.IMGT.Nb.of.nonsilent.mutations", + "FR3.IMGT.Nb.of.nonsilent.mutations") + +print("Cleaning up columns") + +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]),col] = 0 +} + +regions = c("FR1", "CDR1", "FR2", "CDR2", "FR3") +if(empty.region.filter == "FR1") { + regions = c("CDR1", "FR2", "CDR2", "FR3") +} else if (empty.region.filter == "CDR1") { + regions = c("FR2", "CDR2", "FR3", "CDR3") +} else if (empty.region.filter == "FR2") { + regions = c("CDR2", "FR3", "CDR3") +} + +sum_by_row = function(x, columns) { sum(as.numeric(x[columns]), na.rm=T) } + +print("aggregating data into new columns") + +VRegionMutations_columns = paste(regions, ".IMGT.Nb.of.mutations", sep="") +dat$VRegionMutations = apply(dat, FUN=sum_by_row, 1, columns=VRegionMutations_columns) + +VRegionNucleotides_columns = paste(regions, ".IMGT.Nb.of.nucleotides", sep="") +dat$FR3.IMGT.Nb.of.nucleotides = nchar(dat$FR3.IMGT.seq) +dat$VRegionNucleotides = apply(dat, FUN=sum_by_row, 1, columns=VRegionNucleotides_columns) + +transitionMutations_columns = paste(rep(regions, each=4), c(".IMGT.a.g", ".IMGT.g.a", ".IMGT.c.t", ".IMGT.t.c"), sep="") +dat$transitionMutations = apply(dat, FUN=sum_by_row, 1, columns=transitionMutations_columns) + +transversionMutations_columns = paste(rep(regions, each=8), c(".IMGT.a.c",".IMGT.c.a",".IMGT.a.t",".IMGT.t.a",".IMGT.g.c",".IMGT.c.g",".IMGT.g.t",".IMGT.t.g"), sep="") +dat$transversionMutations = apply(dat, FUN=sum_by_row, 1, columns=transversionMutations_columns) + +transitionMutationsAtGC_columns = paste(rep(regions, each=2), c(".IMGT.g.a",".IMGT.c.t"), sep="") +dat$transitionMutationsAtGC = apply(dat, FUN=sum_by_row, 1, columns=transitionMutationsAtGC_columns) + +totalMutationsAtGC_columns = paste(rep(regions, each=6), c(".IMGT.c.g",".IMGT.c.t",".IMGT.c.a",".IMGT.g.c",".IMGT.g.a",".IMGT.g.t"), sep="") +#totalMutationsAtGC_columns = paste(rep(regions, each=6), c(".IMGT.g.a",".IMGT.c.t",".IMGT.c.a",".IMGT.c.g",".IMGT.g.t"), sep="") +dat$totalMutationsAtGC = apply(dat, FUN=sum_by_row, 1, columns=totalMutationsAtGC_columns) + +transitionMutationsAtAT_columns = paste(rep(regions, each=2), c(".IMGT.a.g",".IMGT.t.c"), sep="") +dat$transitionMutationsAtAT = apply(dat, FUN=sum_by_row, 1, columns=transitionMutationsAtAT_columns) + +totalMutationsAtAT_columns = paste(rep(regions, each=6), c(".IMGT.a.g",".IMGT.a.c",".IMGT.a.t",".IMGT.t.g",".IMGT.t.c",".IMGT.t.a"), sep="") +#totalMutationsAtAT_columns = paste(rep(regions, each=5), c(".IMGT.a.g",".IMGT.t.c",".IMGT.a.c",".IMGT.g.c",".IMGT.t.g"), sep="") +dat$totalMutationsAtAT = apply(dat, FUN=sum_by_row, 1, columns=totalMutationsAtAT_columns) + +FRRegions = regions[grepl("FR", regions)] +CDRRegions = regions[grepl("CDR", regions)] + +FR_silentMutations_columns = paste(FRRegions, ".IMGT.Nb.of.silent.mutations", sep="") +dat$silentMutationsFR = apply(dat, FUN=sum_by_row, 1, columns=FR_silentMutations_columns) + +CDR_silentMutations_columns = paste(CDRRegions, ".IMGT.Nb.of.silent.mutations", sep="") +dat$silentMutationsCDR = apply(dat, FUN=sum_by_row, 1, columns=CDR_silentMutations_columns) + +FR_nonSilentMutations_columns = paste(FRRegions, ".IMGT.Nb.of.nonsilent.mutations", sep="") +dat$nonSilentMutationsFR = apply(dat, FUN=sum_by_row, 1, columns=FR_nonSilentMutations_columns) + +CDR_nonSilentMutations_columns = paste(CDRRegions, ".IMGT.Nb.of.nonsilent.mutations", sep="") +dat$nonSilentMutationsCDR = apply(dat, FUN=sum_by_row, 1, columns=CDR_nonSilentMutations_columns) + +mutation.sum.columns = c("Sequence.ID", "VRegionMutations", "VRegionNucleotides", "transitionMutations", "transversionMutations", "transitionMutationsAtGC", "transitionMutationsAtAT", "silentMutationsFR", "nonSilentMutationsFR", "silentMutationsCDR", "nonSilentMutationsCDR") +write.table(dat[,mutation.sum.columns], "mutation_by_id.txt", sep="\t",quote=F,row.names=F,col.names=T) + +setwd(outputdir) + +base.order = data.frame(base=c("A", "T", "C", "G"), order=1:4) + +calculate_result = function(i, gene, dat, matrx, f, fname, name){ + tmp = dat[grepl(paste("^", gene, ".*", sep=""), dat$best_match),] + + j = i - 1 + x = (j * 3) + 1 + y = (j * 3) + 2 + z = (j * 3) + 3 + + if(nrow(tmp) > 0){ + if(fname == "sum"){ + matrx[1,x] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) + matrx[1,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1) + matrx[1,z] = round(f(matrx[1,x] / matrx[1,y]) * 100, digits=1) + } else { + matrx[1,x] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) + matrx[1,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1) + matrx[1,z] = round(f(tmp$VRegionMutations / tmp$VRegionNucleotides) * 100, digits=1) + } + + matrx[2,x] = round(f(tmp$transitionMutations, na.rm=T), digits=1) + matrx[2,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) + matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1) + + matrx[3,x] = round(f(tmp$transversionMutations, na.rm=T), digits=1) + matrx[3,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) + matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1) + + matrx[4,x] = round(f(tmp$transitionMutationsAtGC, na.rm=T), digits=1) + matrx[4,y] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1) + matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1) + + matrx[5,x] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1) + matrx[5,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) + matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1) + + matrx[6,x] = round(f(tmp$transitionMutationsAtAT, na.rm=T), digits=1) + matrx[6,y] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1) + matrx[6,z] = round(matrx[6,x] / matrx[6,y] * 100, digits=1) + + matrx[7,x] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1) + matrx[7,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) + matrx[7,z] = round(matrx[7,x] / matrx[7,y] * 100, digits=1) + + matrx[8,x] = round(f(tmp$nonSilentMutationsFR, na.rm=T), digits=1) + matrx[8,y] = round(f(tmp$silentMutationsFR, na.rm=T), digits=1) + matrx[8,z] = round(matrx[8,x] / matrx[8,y], digits=1) + + matrx[9,x] = round(f(tmp$nonSilentMutationsCDR, na.rm=T), digits=1) + matrx[9,y] = round(f(tmp$silentMutationsCDR, na.rm=T), digits=1) + matrx[9,z] = round(matrx[9,x] / matrx[9,y], digits=1) + + if(fname == "sum"){ + matrx[10,x] = round(f(rowSums(tmp[,c("FR2.IMGT.Nb.of.nucleotides", "FR3.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1) + matrx[10,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1) + matrx[10,z] = round(matrx[10,x] / matrx[10,y] * 100, digits=1) + + matrx[11,x] = round(f(rowSums(tmp[,c("CDR1.IMGT.Nb.of.nucleotides", "CDR2.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1) + matrx[11,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1) + matrx[11,z] = round(matrx[11,x] / matrx[11,y] * 100, digits=1) + } + } + + transitionTable = data.frame(A=zeros,C=zeros,G=zeros,T=zeros) + row.names(transitionTable) = c("A", "C", "G", "T") + transitionTable["A","A"] = NA + transitionTable["C","C"] = NA + transitionTable["G","G"] = NA + transitionTable["T","T"] = NA + + if(nrow(tmp) > 0){ + 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="") + if (empty.region.filter == "leader"){ + transitionTable[NT1,NT2] = sum(tmp[,c(FR1, CDR1, FR2, CDR2, FR3)]) + } else if (empty.region.filter == "FR1") { + transitionTable[NT1,NT2] = sum(tmp[,c(CDR1, FR2, CDR2, FR3)]) + } else if (empty.region.filter == "CDR1") { + transitionTable[NT1,NT2] = sum(tmp[,c(FR2, CDR2, FR3)]) + } else if (empty.region.filter == "FR2") { + transitionTable[NT1,NT2] = sum(tmp[,c(CDR2, FR3)]) + } + } + } + transition = transitionTable + transition$id = names(transition) + + transition2 = melt(transition, id.vars="id") + + transition2 = merge(transition2, base.order, by.x="id", by.y="base") + + transition2 = merge(transition2, base.order, by.x="variable", by.y="base") + + transition2[is.na(transition2$value),]$value = 0 + + if(any(transition2$value == 0)){ #having rows of data but a transition table filled with 0 is bad + print("Plotting stacked transition") + png(filename=paste("transitions_stacked_", name, ".png", sep="")) + p = ggplot(transition2, aes(factor(reorder(id, order.x)), y=value, fill=factor(reorder(variable, order.y)))) + geom_bar(position="fill", stat="identity", colour="black") #stacked bar + p = p + xlab("From base") + ylab("To base") + ggtitle("Mutations frequency from base to base") + guides(fill=guide_legend(title=NULL)) + p = p + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=13, colour="black")) + scale_fill_manual(values=c("A" = "blue4", "G" = "lightblue1", "C" = "olivedrab3", "T" = "olivedrab4")) + #p = p + scale_colour_manual(values=c("A" = "black", "G" = "black", "C" = "black", "T" = "black")) + print(p) + dev.off() + + print("Plotting heatmap transition") + + png(filename=paste("transitions_heatmap_", name, ".png", sep="")) + p = ggplot(transition2, aes(factor(reorder(id, order.x)), factor(reorder(variable, order.y)))) + geom_tile(aes(fill = value)) + scale_fill_gradient(low="white", high="steelblue") #heatmap + p = p + xlab("From base") + ylab("To base") + ggtitle("Mutations frequency from base to base") + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=13, colour="black")) + print(p) + dev.off() + } else { + print("No data to plot") + } + } + + #print(paste("writing value file: ", name, "_", fname, "_value.txt" ,sep="")) + write.table(x=transitionTable, file=paste("transitions_", name ,"_", fname, ".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_", name , "_", fname, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T) + cat(matrx[1,x], file=paste(name, "_", fname, "_value.txt" ,sep="")) + cat(nrow(tmp), file=paste(name, "_", fname, "_n.txt" ,sep="")) + #print(paste(fname, name, nrow(tmp))) + matrx +} +nts = c("a", "c", "g", "t") +zeros=rep(0, 4) +funcs = c(median, sum, mean) +fnames = c("median", "sum", "mean") + +print("Creating result tables") + +for(i in 1:length(funcs)){ + func = funcs[[i]] + fname = fnames[[i]] + + rows = 9 + if(fname == "sum"){ + rows = 11 + } + matrx = matrix(data = 0, ncol=((length(genes) + 1) * 3),nrow=rows) + for(i in 1:length(genes)){ + print(paste("Creating table for", fname, genes[i])) + matrx = calculate_result(i, genes[i], dat, matrx, func, fname, genes[i]) + } + matrx = calculate_result(i + 1, ".*", dat[!grepl("unmatched", dat$best_match),], matrx, func, fname, name="all") + + result = data.frame(matrx) + if(fname == "sum"){ + row.names(result) = c("Number of Mutations (%)", "Transitions (%)", "Transversions (%)", "Transitions at G C (%)", "Targeting of C G (%)", "Transitions at A T (%)", "Targeting of A T (%)", "FR R/S (ratio)", "CDR R/S (ratio)", "nt in FR", "nt in CDR") + } else { + row.names(result) = c("Number of Mutations (%)", "Transitions (%)", "Transversions (%)", "Transitions at G C (%)", "Targeting of C G (%)", "Transitions at A T (%)", "Targeting of A T (%)", "FR R/S (ratio)", "CDR R/S (ratio)") + } + write.table(x=result, file=paste("mutations_", fname, ".txt", sep=""), sep=",",quote=F,row.names=T,col.names=F) +} + +print("Adding median number of mutations to sum table") +sum.table = read.table("mutations_sum.txt", sep=",", header=F) +median.table = read.table("mutations_median.txt", sep=",", header=F) + +new.table = sum.table[1,] +new.table[2,] = median.table[1,] +new.table[3:12,] = sum.table[2:11,] +new.table[,1] = as.character(new.table[,1]) +new.table[2,1] = "Median of Number of Mutations (%)" + +#sum.table = sum.table[c("Number of Mutations (%)", "Median of Number of Mutations (%)", "Transition (%)", "Transversions (%)", "Transitions at G C (%)", "Targeting of C G (%)", "Transitions at A T (%)", "Targeting of A T (%)", "FR R/S (ratio)", "CDR R/S (ratio)", "nt in FR", "nt in CDR"),] + +write.table(x=new.table, file="mutations_sum.txt", sep=",",quote=F,row.names=F,col.names=F) + +print("Plotting IGA piechart") + +dat = dat[!grepl("^unmatched", dat$best_match),] + +#blegh + +genesForPlot = dat[grepl("IGA", 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=Gene)) + pc = pc + geom_bar(width = 1, stat = "identity") + scale_fill_manual(labels=genesForPlot$label, values=c("IGA1" = "lightblue1", "IGA2" = "blue4")) + pc = pc + coord_polar(theta="y") + scale_y_continuous(breaks=NULL) + pc = pc + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=13, colour="black")) + pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IGA subclasses", "( n =", sum(genesForPlot$Freq), ")")) + write.table(genesForPlot, "IGA.txt", sep="\t",quote=F,row.names=F,col.names=T) + + png(filename="IGA.png") + print(pc) + dev.off() +} + +print("Plotting IGG piechart") + +genesForPlot = dat[grepl("IGG", 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=Gene)) + pc = pc + geom_bar(width = 1, stat = "identity") + scale_fill_manual(labels=genesForPlot$label, values=c("IGG1" = "olivedrab3", "IGG2" = "red", "IGG3" = "gold", "IGG4" = "darkred")) + pc = pc + coord_polar(theta="y") + scale_y_continuous(breaks=NULL) + pc = pc + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=13, colour="black")) + pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IGG subclasses", "( n =", sum(genesForPlot$Freq), ")")) + write.table(genesForPlot, "IGG.txt", sep="\t",quote=F,row.names=F,col.names=T) + + png(filename="IGG.png") + print(pc) + dev.off() +} + +print("Plotting scatterplot") + +dat$percentage_mutations = round(dat$VRegionMutations / dat$VRegionNucleotides * 100, 2) +dat.clss = dat + +dat.clss$best_match = substr(dat.clss$best_match, 0, 3) + +dat.clss = rbind(dat, dat.clss) + +p = ggplot(dat.clss, aes(best_match, percentage_mutations)) +p = p + geom_point(aes(colour=best_match), position="jitter") + geom_boxplot(aes(middle=mean(percentage_mutations)), alpha=0.1, outlier.shape = NA) +p = p + xlab("Subclass") + ylab("Frequency") + ggtitle("Frequency scatter plot") + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=13, colour="black")) +p = p + scale_fill_manual(values=c("IGA" = "blue4", "IGA1" = "lightblue1", "IGA2" = "blue4", "IGG" = "olivedrab3", "IGG1" = "olivedrab3", "IGG2" = "red", "IGG3" = "gold", "IGG4" = "darkred", "IGM" = "darkviolet")) +p = p + scale_colour_manual(values=c("IGA" = "blue4", "IGA1" = "lightblue1", "IGA2" = "blue4", "IGG" = "olivedrab3", "IGG1" = "olivedrab3", "IGG2" = "red", "IGG3" = "gold", "IGG4" = "darkred", "IGM" = "darkviolet")) + +png(filename="scatter.png") +print(p) +dev.off() + +write.table(dat[,c("Sequence.ID", "best_match", "VRegionMutations", "VRegionNucleotides", "percentage_mutations")], "scatter.txt", sep="\t",quote=F,row.names=F,col.names=T) + +write.table(dat, input, sep="\t",quote=F,row.names=F,col.names=T) + +print("Plotting frequency ranges plot") + +dat$best_match_class = substr(dat$best_match, 0, 3) +freq_labels = c("0", "0-2", "2-5", "5-10", "10-15", "15-20", "20") +dat$frequency_bins = cut(dat$percentage_mutations, breaks=c(-Inf, 0, 2,5,10,15,20, Inf), labels=freq_labels) + +frequency_bins_sum = data.frame(data.table(dat)[, list(class_sum=sum(.N)), by=c("best_match_class")]) + +frequency_bins_data = data.frame(data.table(dat)[, list(frequency_count=.N), by=c("best_match_class", "frequency_bins")]) + +frequency_bins_data = merge(frequency_bins_data, frequency_bins_sum, by="best_match_class") + +frequency_bins_data$frequency = round(frequency_bins_data$frequency_count / frequency_bins_data$class_sum * 100, 2) + +p = ggplot(frequency_bins_data, aes(frequency_bins, frequency)) +p = p + geom_bar(aes(fill=best_match_class), stat="identity", position="dodge") + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=13, colour="black")) +p = p + xlab("Frequency ranges") + ylab("Frequency") + ggtitle("Mutation Frequencies by class") + scale_fill_manual(values=c("IGA" = "blue4", "IGG" = "olivedrab3", "IGM" = "black")) + +png(filename="frequency_ranges.png") +print(p) +dev.off() + +frequency_bins_data_by_class = frequency_bins_data + +write.table(frequency_bins_data_by_class, "frequency_ranges_classes.txt", sep="\t",quote=F,row.names=F,col.names=T) + +frequency_bins_data = data.frame(data.table(dat)[, list(frequency_count=.N), by=c("best_match", "best_match_class", "frequency_bins")]) + +frequency_bins_data = merge(frequency_bins_data, frequency_bins_sum, by="best_match_class") + +frequency_bins_data$frequency = round(frequency_bins_data$frequency_count / frequency_bins_data$class_sum * 100, 2) + +write.table(frequency_bins_data, "frequency_ranges_subclasses.txt", sep="\t",quote=F,row.names=F,col.names=T) + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
--- a/shm_csr.xml Thu Oct 13 10:52:24 2016 -0400 +++ b/shm_csr.xml Tue Oct 25 07:28:43 2016 -0400 @@ -1,24 +1,20 @@ <tool id="shm_csr" name="SHM & CSR pipeline" version="1.0"> <description></description> <command interpreter="bash"> - wrapper.sh $in_file custom $out_file $out_file.files_path ${in_file.name} ${include_fr1} $functionality $unique $naive_output_ca $naive_output_cg $naive_output_cm $filter_uniques $class_filter $empty_region_filter + wrapper.sh $in_file custom $out_file $out_file.files_path ${in_file.name} "-" $functionality $unique $naive_output_ca $naive_output_cg $naive_output_cm $filter_uniques $class_filter $empty_region_filter </command> <inputs> <param name="in_file" type="data" label="IMGT zip file to be analysed" /> <param name="empty_region_filter" type="select" label="Sequence starts at" help="" > + <option value="leader" selected="true">Leader: include FR1, CDR1, FR2, CDR2, FR3 in filters</option> <option value="FR1" selected="true">FR1: include CDR1,FR2,CDR2,FR3 in filters</option> <option value="CDR1">CDR1: include FR2,CDR2,FR3 in filters</option> <option value="FR2">FR2: include CDR2,FR3 in filters</option> </param> - <param name="include_fr1" type="select" label="Include mutations in FR1 region" help="" > - <option value="yes">yes</option> - <option value="no" selected="true">no</option> - </param> <param name="functionality" type="select" label="Functionality filter" help="" > - <option value="productive" selected="true">Productive: Keep "productive" and "productive (see comment)"</option> - <option value="unproductive">Unproductive: Keep "unproductive" and "unproductive (see comment)"</option> - <option value="remove_unknown">Remove "unknown" and "unknown (see comment)"</option> - <option value="dont_filter">Don't filter</option> + <option value="productive" selected="true">Productive (Productive and Productive see comment)</option> + <option value="unproductive">Unproductive (Unproductive and Unproductive see comment)</option> + <option value="remove_unknown">Productive and Unproductive (Productive, Productive see comment, Unproductive, Unproductive and Unproductive see comment)</option> </param> <param name="filter_uniques" type="select" label="Filter unique sequences" help="See below for an example."> <option value="remove">Remove uniques (Based on nucleotide sequence + C)</option> @@ -42,7 +38,7 @@ <option value="60_55">>60% class and >55% subclass</option> <option value="70_0">>70% class</option> <option value="60_0">>60% class</option> - <option value="101_101">No class</option> + <option value="101_101">Do not assign (sub)class</option> </param> <conditional name="naive_output_cond"> <param name="naive_output" type="select" label="Output new IMGT archives per class into your history?">
--- a/wrapper.sh Thu Oct 13 10:52:24 2016 -0400 +++ b/wrapper.sh Tue Oct 25 07:28:43 2016 -0400 @@ -158,13 +158,13 @@ classes="IGA,IGA1,IGA2,IGG,IGG1,IGG2,IGG3,IGG4,IGM,unmatched" echo "R mutation analysis" -Rscript $dir/shm_csr.r $outdir/merged.txt $classes $outdir ${include_fr1} 2>&1 +Rscript $dir/shm_csr.r $outdir/merged.txt $classes $outdir ${empty_region_filter} 2>&1 echo "---------------- shm_csr.py ----------------" echo "---------------- shm_csr.py ----------------<br />" >> $log -python $dir/shm_csr.py --input $outdir/merged.txt --genes $classes --includefr1 "${include_fr1}" --output $outdir/hotspot_analysis.txt +python $dir/shm_csr.py --input $outdir/merged.txt --genes $classes --empty_region_filter "${empty_region_filter}" --output $outdir/hotspot_analysis.txt echo "---------------- aa_histogram.r ----------------" echo "---------------- aa_histogram.r ----------------<br />" >> $log @@ -575,6 +575,9 @@ cd $tmp +echo "Cleaning up *.RData files" +find $outdir/baseline -name "*.RData" -type f -delete + echo "---------------- naive_output.r ----------------" echo "---------------- naive_output.r ----------------<br />" >> $log