Mercurial > repos > davidvanzessen > mutation_analysis
changeset 82:564c4f6da203 draft
Uploaded
| author | davidvanzessen |
|---|---|
| date | Tue, 24 May 2016 03:52:14 -0400 |
| parents | a778156dad3d |
| children | e1114c791eb7 |
| files | aa_histogram.r merge_and_filter.r sequence_overview.r wrapper.sh |
| diffstat | 4 files changed, 71 insertions(+), 28 deletions(-) [+] |
line wrap: on
line diff
--- a/aa_histogram.r Tue May 17 04:17:38 2016 -0400 +++ b/aa_histogram.r Tue May 24 03:52:14 2016 -0400 @@ -34,7 +34,7 @@ print(paste("dat_dt", nrow(dat_dt))) m = ggplot(dat_dt, aes(x=i, y=freq)) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) -m = m + geom_histogram(stat="identity", colour = "black", fill = "darkgrey", alpha=0.8) + scale_x_continuous(breaks=1:length(dat_freq), labels=1:length(dat_freq)) +m = m + geom_bar(stat="identity", colour = "black", fill = "darkgrey", alpha=0.8) + scale_x_continuous(breaks=1:length(dat_freq), labels=1:length(dat_freq)) 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") m = m + annotate("segment", x = 38.5, y = -0.05, xend=55.5, yend=-0.05, colour="darkgreen", size=1) + annotate("text", x = 47, y = -0.1, label="FR2")
--- a/merge_and_filter.r Tue May 17 04:17:38 2016 -0400 +++ b/merge_and_filter.r Tue May 24 03:52:14 2016 -0400 @@ -59,7 +59,7 @@ unmatched=summ[NULL,c("Sequence.ID", "chunk_hit_percentage", "nt_hit_percentage", "start_locations", "best_match")] -if(!all(higher_than)){ #check for 'not all' because that would mean the unmatched set is empty +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)
--- a/sequence_overview.r Tue May 17 04:17:38 2016 -0400 +++ b/sequence_overview.r Tue May 24 03:52:14 2016 -0400 @@ -4,7 +4,10 @@ gene.matches = args[1] sequence.file = args[2] -outputdir = args[3] +merged.file = args[3] +outputdir = args[4] +gene.classes = unlist(strsplit(args[5], ",")) +hotspot.analysis.sum.file = args[6] NToverview.file = paste(outputdir, "ntoverview.txt", sep="/") NTsum.file = paste(outputdir, "ntsum.txt", sep="/") main.html = "index.html" @@ -13,9 +16,13 @@ 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="") +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="") dat = merge(sequences, genes, by="Sequence.ID") +dat = dat[dat$Sequence.ID %in% merged$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")] @@ -23,15 +30,15 @@ #dat = data.frame(data.table(dat)[, list(freq=.N), by=c("best_match", "seq_conc")]) -dat = data.frame(table(dat$seq_conc)) +dat = data.frame(table(dat$seq_conc, dat$Functionality)) dat = dat[dat$Freq > 1,] -names(dat) = c("seq_conc", "Freq") +names(dat) = c("seq_conc", "Functionality", "Freq") dat$seq_conc = factor(dat$seq_conc) -dat = dat[order(nchar(as.character(dat$seq_conc))),] +dat = dat[order(as.character(dat$seq_conc)),] #writing html from R... td = function(val) { paste("<td>", val, "</td>", sep="") } @@ -54,6 +61,14 @@ cm = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "cm",] + classes = c(nrow(ca1), nrow(ca2), nrow(cg1), nrow(cg2), nrow(cg3), nrow(cg4), nrow(cm)) + + classes.sum = sum(classes) + + if(!any(classes != 0 & classes != classes.sum)){ + next + } + id = as.numeric(dat[i,"seq_conc"]) functionality = dat[i,"Functionality"] @@ -113,20 +128,48 @@ NToverview = merge(NToverview, sequences[,c("Sequence.ID", "seq")], by="Sequence.ID") +NToverview = NToverview[NToverview$Sequence.ID %in% merged$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)) +#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){ + NToverview.sub = NToverview[grepl(clazz, NToverview$best_match),] + 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="")) +} -print(names(NToverview)) -print(names(NTsum)) +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="")) -NToverview = rbind(NToverview, NTsum) +names(hotspot.analysis.sum) = names(NTresult) + +hotspot.analysis.sum = rbind(hotspot.analysis.sum, NTresult) + +print(hotspot.analysis.sum) + +write.table(hotspot.analysis.sum, hotspot.analysis.sum.file, quote=F, sep=",", row.names=F, col.names=F, na="0") + 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)
--- a/wrapper.sh Tue May 17 04:17:38 2016 -0400 +++ b/wrapper.sh Tue May 24 03:52:14 2016 -0400 @@ -86,9 +86,9 @@ echo "---------------- mutation_analysis.r ----------------" echo "---------------- mutation_analysis.r ----------------<br />" >> $output -genes="ca,ca1,ca2,cg,cg1,cg2,cg3,cg4,cm,unmatched" +classes="ca,ca1,ca2,cg,cg1,cg2,cg3,cg4,cm,unmatched" echo "R mutation analysis" -Rscript $dir/mutation_analysis.r $outdir/merged.txt $genes $outdir ${include_fr1} 2>&1 +Rscript $dir/mutation_analysis.r $outdir/merged.txt $classes $outdir ${include_fr1} 2>&1 #echo "." > $output #exit 0 @@ -98,7 +98,7 @@ echo "---------------- mutation_analysis.py ----------------" echo "---------------- mutation_analysis.py ----------------<br />" >> $output -python $dir/mutation_analysis.py --input $outdir/merged.txt --genes $genes --includefr1 "${include_fr1}" --output $outdir/hotspot_analysis.txt +python $dir/mutation_analysis.py --input $outdir/merged.txt --genes $classes --includefr1 "${include_fr1}" --output $outdir/hotspot_analysis.txt echo "---------------- aa_histogram.r ----------------" echo "---------------- aa_histogram.r ----------------<br />" >> $output @@ -109,6 +109,19 @@ funcs=(sum mean median) +echo "---------------- sequence_overview.r ----------------" + +mkdir $outdir/sequence_overview + +Rscript $dir/sequence_overview.r $outdir/identified_genes.txt $PWD/sequences.txt $outdir/merged.txt $outdir/sequence_overview $classes $outdir/hotspot_analysis_sum.txt 2>&1 + +echo "<table border='1'>" > $outdir/base_overview.html + +while read ID class seq A C G T +do + echo "<tr><td>$ID</td><td>$class</td><td>$A</td><td>$C</td><td>$G</td><td>$T</td></tr>" >> $outdir/base_overview.html +done < $outdir/sequence_overview/ntoverview.txt + echo "<html><center><h1>$title</h1></center>" > $output @@ -222,7 +235,7 @@ #optional output for naive -echo "---------------- aa_histogram.r ----------------" +echo "---------------- naive_output.r ----------------" if [[ "$naive_output" != "None" ]] then @@ -233,19 +246,6 @@ Rscript $dir/naive_output.r $outdir/loader_output.txt $outdir/merged.txt ${naive_output_ca} ${naive_output_cg} ${naive_output_cm} $outdir/ntoverview.txt $outdir/ntsum.txt 2>&1 fi -echo "---------------- sequence_overview.r ----------------" - -mkdir $outdir/sequence_overview - -Rscript $dir/sequence_overview.r $outdir/identified_genes.txt $PWD/sequences.txt $outdir/sequence_overview 2>&1 - -echo "<table border='1'>" > $outdir/base_overview.html - -while read ID class seq A C G T -do - echo "<tr><td>$ID</td><td>$class</td><td>$A</td><td>$C</td><td>$G</td><td>$T</td></tr>" >> $outdir/base_overview.html -done < $outdir/sequence_overview/ntoverview.txt - echo "</table>" >> $outdir/base_overview.html echo "---------------- Done! ----------------"