Mercurial > repos > davidvanzessen > test_plotting_merged
changeset 0:0e6f33e21c60 draft
Uploaded
| author | davidvanzessen |
|---|---|
| date | Wed, 09 Oct 2013 04:31:48 -0400 |
| parents | |
| children | c43c63f4ef80 |
| files | RScript.r plotting_merged.xml r_wrapper.sh |
| diffstat | 3 files changed, 224 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/RScript.r Wed Oct 09 04:31:48 2013 -0400 @@ -0,0 +1,174 @@ +options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } ) + +args <- commandArgs(trailingOnly = TRUE) + +inFile = args[1] +outFile = args[2] +outDir = args[3] + +if (!require("gridExtra")) { +install.packages("gridExtra", repos="http://cran.xl-mirror.nl/") +} +library (gridExtra) +if (!require("ggplot2")) { +install.packages("ggplot2", repos="http://cran.xl-mirror.nl/") +} +require(ggplot2) +if (!require("plyr")) { +install.packages("plyr", repos="http://cran.xl-mirror.nl/") +} +require(plyr) + +test = read.table(inFile, sep="\t", header=TRUE) + +test$Top.V.Gene = gsub("[*]([0-9]+)", "", test$Top.V.Gene) +test$Top.D.Gene = gsub("[*]([0-9]+)", "", test$Top.D.Gene) +test$Top.J.Gene = gsub("[*]([0-9]+)", "", test$Top.J.Gene) + +test$VDJCDR3 = do.call(paste, c(test[c("Top.V.Gene", "Top.D.Gene", "Top.J.Gene","CDR3.Seq.DNA")], sep = ":")) + +PROD = test[test$VDJ.Frame != "In-frame with stop codon" & test$VDJ.Frame != "Out-of-frame" & test$CDR3.Found.How != "NOT_FOUND" , ] + +NONPROD = test[test$VDJ.Frame == "In-frame with stop codon" | test$VDJ.Frame == "Out-of-frame" | test$CDR3.Found.How == "NOT_FOUND" , ] + +PRODF = PROD[ -1] + +#unique(PRODF[duplicated(PRODF),]) +#length(row.names(PRODF[duplicated(PRODF),])) + +#length(row.names(PRODF)) +PRODF = unique(PRODF) +#length(row.names(PRODF)) + +PRODFV = ddply(PRODF, c("Sample", "Top.V.Gene"), function(x) summary(x$VDJCDR3)) +PRODFV$Length = as.numeric(PRODFV$Length) +Total = 0 +Total = ddply(PRODFV, .(Sample), function(x) data.frame(Total = sum(x$Length))) +PRODFV = merge(PRODFV, Total, by.x='Sample', by.y='Sample', all.x=TRUE) +PRODFV = ddply(PRODFV, c("Sample", "Top.V.Gene"), summarise, relFreq= (Length*100 / Total)) + +PRODFD = ddply(PRODF, c("Sample", "Top.D.Gene"), function(x) summary(x$VDJCDR3)) +PRODFD$Length = as.numeric(PRODFD$Length) +Total = 0 +Total = ddply(PRODFD, .(Sample), function(x) data.frame(Total = sum(x$Length))) +PRODFD = merge(PRODFD, Total, by.x='Sample', by.y='Sample', all.x=TRUE) +PRODFD = ddply(PRODFD, c("Sample", "Top.D.Gene"), summarise, relFreq= (Length*100 / Total)) + +PRODFJ = ddply(PRODF, c("Sample", "Top.J.Gene"), function(x) summary(x$VDJCDR3)) +PRODFJ$Length = as.numeric(PRODFJ$Length) +Total = 0 +Total = ddply(PRODFJ, .(Sample), function(x) data.frame(Total = sum(x$Length))) +PRODFJ = merge(PRODFJ, Total, by.x='Sample', by.y='Sample', all.x=TRUE) +PRODFJ = ddply(PRODFJ, c("Sample", "Top.J.Gene"), summarise, relFreq= (Length*100 / Total)) + +V = c("v.name\tchr.orderV\nIGHV7-81\t1\nIGHV3-74\t2\nIGHV3-73\t3\nIGHV3-72\t4\nIGHV3-71\t5\nIGHV2-70\t6\nIGHV1-69\t7\nIGHV3-66\t8\nIGHV3-64\t9\nIGHV4-61\t10\nIGHV4-59\t11\nIGHV1-58\t12\nIGHV3-53\t13\nIGHV3-52\t14\nIGHV5-a\t15\nIGHV5-51\t16\nIGHV3-49\t17\nIGHV3-48\t18\nIGHV3-47\t19\nIGHV1-46\t20\nIGHV1-45\t21\nIGHV3-43\t22\nIGHV4-39\t23\nIGHV3-35\t24\nIGHV4-34\t25\nIGHV3-33\t26\nIGHV4-31\t27\nIGHV4-30-4\t28\nIGHV4-30-2\t29\nIGHV3-30-3\t30\nIGHV3-30\t31\nIGHV4-28\t32\nIGHV2-26\t33\nIGHV1-24\t34\nIGHV3-23\t35\nIGHV3-22\t36\nIGHV3-21\t37\nIGHV3-20\t38\nIGHV3-19\t39\nIGHV1-18\t40\nIGHV3-15\t41\nIGHV3-13\t42\nIGHV3-11\t43\nIGHV3-9\t44\nIGHV1-8\t45\nIGHV3-7\t46\nIGHV2-5\t47\nIGHV7-4-1\t48\nIGHV4-4\t49\nIGHV4-b\t50\nIGHV1-3\t51\nIGHV1-2\t52\nIGHV6-1\t53") +tcV = textConnection(V) +Vchain = read.table(tcV, sep="\t", header=TRUE) +PRODFV = merge(PRODFV, Vchain, by.x='Top.V.Gene', by.y='v.name', all.x=TRUE) +close(tcV) + +D = c("v.name\tchr.orderD\nIGHD1-1\t1\nIGHD2-2\t2\nIGHD3-3\t3\nIGHD6-6\t4\nIGHD1-7\t5\nIGHD2-8\t6\nIGHD3-9\t7\nIGHD3-10\t8\nIGHD4-11\t9\nIGHD5-12\t10\nIGHD6-13\t11\nIGHD1-14\t12\nIGHD2-15\t13\nIGHD3-16\t14\nIGHD4-17\t15\nIGHD5-18\t16\nIGHD6-19\t17\nIGHD1-20\t18\nIGHD2-21\t19\nIGHD3-22\t20\nIGHD4-23\t21\nIGHD5-24\t22\nIGHD6-25\t23\nIGHD1-26\t24\nIGHD7-27\t25") +tcD = textConnection(D) +Dchain = read.table(tcD, sep="\t", header=TRUE) +PRODFD = merge(PRODFD, Dchain, by.x='Top.D.Gene', by.y='v.name', all.x=TRUE) +close(tcD) + + +J = c("v.name\tchr.orderJ\nIGHJ1\t1\nIGHJ2\t2\nIGHJ3\t3\nIGHJ4\t4\nIGHJ5\t5\nIGHJ6\t6") +tcJ = textConnection(J) +Jchain = read.table(tcJ, sep="\t", header=TRUE) +PRODFJ = merge(PRODFJ, Jchain, by.x='Top.J.Gene', by.y='v.name', all.x=TRUE) +close(tcJ) + +setwd(outDir) + +pV = ggplot(PRODFV) +pV = pV + geom_bar( aes( x=factor(reorder(Top.V.Gene, chr.orderV)), y=relFreq, fill=Sample), stat='identity', position="dodge") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + +png("VPlot.png",width = 1280, height = 720) +pV +dev.off(); + +pD = ggplot(PRODFD) +pD = pD + geom_bar( aes( x=factor(reorder(Top.D.Gene, chr.orderD)), y=relFreq, fill=Sample), stat='identity', position="dodge") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + +png("DPlot.png",width = 800, height = 600) +pD +dev.off(); + +pJ = ggplot(PRODFJ) +pJ = pJ + geom_bar( aes( x=factor(reorder(Top.J.Gene, chr.orderJ)), y=relFreq, fill=Sample), stat='identity', position="dodge") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + +png("JPlot.png",width = 800, height = 600) +pJ +dev.off(); + + +plotVD <- function(dat){ + ggplot() + + geom_tile(data=dat, aes(x=factor(Top.V.Gene), y=factor(Top.D.Gene), fill=log)) + + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + + scale_fill_gradient(low="white", high="red") + + ggtitle(unique(dat$Sample)) +} + + +VandDCount = ddply(PRODF, c("Top.V.Gene", "Top.D.Gene", "Sample"), function(x) summary(x$VDJCDR3)) +cartegianProductVD = expand.grid(Top.V.Gene = Vchain$v.name, Top.D.Gene = Dchain$v.name, Sample = unique(test$Sample)) + +completeVD = merge(VandDCount, cartegianProductVD, all.y=TRUE) +completeVD$Length = as.numeric(completeVD$Length) +completeVD$log = log(completeVD$Length) +completeVD$log[is.na(completeVD$log)] = 0 +l = split(completeVD, f=completeVD[,"Sample"]) +png("HeatmapVD%d.png") +lapply(l, FUN=plotVD) +dev.off() + + +plotVJ <- function(dat){ + ggplot() + + geom_tile(data=dat, aes(x=factor(Top.V.Gene), y=factor(Top.J.Gene), fill=log)) + + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + + scale_fill_gradient(low="white", high="red") + + ggtitle(unique(dat$Sample)) +} + +VandJCount = ddply(PRODF, c("Top.V.Gene", "Top.J.Gene", "Sample"), function(x) summary(x$VDJCDR3)) +cartegianProductVJ = expand.grid(Top.V.Gene = Vchain$v.name, Top.J.Gene = Jchain$v.name, Sample = unique(test$Sample)) + +completeVJ = merge(VandJCount, cartegianProductVJ, all.y=TRUE) +completeVJ$Length = as.numeric(completeVJ$Length) +completeVJ$log = log(completeVJ$Length) +completeVJ$log[is.na(completeVJ$log)] = 0 +l = split(completeVJ, f=completeVJ[,"Sample"]) +png("HeatmapVJ%d.png") +lapply(l, FUN=plotVJ) +dev.off() + +plotDJ <- function(dat){ + ggplot() + + geom_tile(data=dat, aes(x=factor(Top.D.Gene), y=factor(Top.J.Gene), fill=log)) + + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + + scale_fill_gradient(low="white", high="red") + + ggtitle(unique(dat$Sample)) +} + +DandJCount = ddply(PRODF, c("Top.D.Gene", "Top.J.Gene", "Sample"), function(x) summary(x$VDJCDR3)) +cartegianProductDJ = expand.grid(Top.D.Gene = Dchain$v.name, Top.J.Gene = Jchain$v.name, Sample = unique(test$Sample)) + +completeDJ = merge(DandJCount, cartegianProductDJ, all.y=TRUE) +completeDJ$Length = as.numeric(completeDJ$Length) +completeDJ$log = log(completeDJ$Length) +completeDJ$log[is.na(completeDJ$log)] = 0 +l = split(completeDJ, f=completeDJ[,"Sample"]) +png("HeatmapDJ%d.png") +lapply(l, FUN=plotDJ) +dev.off() + + +sampleFile <- file("samples.txt") +un = unique(test$Sample) +un = paste(un, sep="\n") +writeLines(un, sampleFile) +close(sampleFile)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/plotting_merged.xml Wed Oct 09 04:31:48 2013 -0400 @@ -0,0 +1,15 @@ +<tool id="test_plotting_merged" name="TEST-Plotting the merged data-TEST" version="1.0"> + <description>Creating a graph of the merged reports.</description> + <command interpreter="bash"> + r_wrapper.sh $in_file $out_file $out_file.files_path + </command> + <inputs> + <param name="in_file" format="tabular" type="data" label="Data to Process" /> + </inputs> + <outputs> + <data format="html" name="out_file" /> + </outputs> + <help> + Creating a graph of the merged and parsed igblast reports. + </help> +</tool>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/r_wrapper.sh Wed Oct 09 04:31:48 2013 -0400 @@ -0,0 +1,35 @@ +#!/bin/bash +echo $1 +echo $2 +echo $3 + +inputFile=$1 +outputFile=$2 +outputDir=$3 +dir="$(cd "$(dirname "$0")" && pwd)" +mkdir $3 +Rscript --verbose $dir/RScript.r $inputFile $outputDir $outputDir 2>&1 +echo "<html>" >> $2 +echo "<img src='VPlot.png'/>" >> $2 +echo "<img src='DPlot.png'/>" >> $2 +echo "<img src='JPlot.png'/>" >> $2 + +samples=`cat $outputDir/samples.txt` +count=1 +echo "<table border='1'>" >> $outputFile +for sample in $samples; do + echo "<tr><td colspan='3' height='100'></td>" >> $outputFile + echo "<tr><td colspan='3'><h1>$sample</h1></td></tr>" >> $outputFile + echo "<tr><td><h2>V-D Heatmap:</h2></td><td><h2>V-J Heatmap:</h2></td><td><h2>D-J Heatmap:</h2></td></tr><tr>" >> $outputFile + mv "$outputDir/HeatmapVD$count.png" "$outputDir/VD_$sample.png" + echo "<td><img src='VD_$sample.png'/></td>" >> $outputFile + mv "$outputDir/HeatmapVJ$count.png" "$outputDir/VJ_$sample.png" + echo "<td><img src='VJ_$sample.png'/></td>" >> $outputFile + mv "$outputDir/HeatmapDJ$count.png" "$outputDir/DJ_$sample.png" + echo "<td><img src='DJ_$sample.png'/></td></tr>" >> $outputFile + count=$((count+1)) +done +echo "</table>" >> $outputFile + +echo "</html>" >> $2 +