extensions/net.sf.basedb.reggie/trunk/src/net/sf/basedb/reggie/baf/bafPlot.R

Code
Comments
Other
Rev Date Author Line
5049 23 Oct 18 nicklas 1 # 
5049 23 Oct 18 nicklas 2 #  R script for generating a plot with mBAF data
5049 23 Oct 18 nicklas 3 #  The plot is divided into two parts:
5049 23 Oct 18 nicklas 4 #  1) Average mBAF per region for a sample. Regions that differ too much (pval < 0.01) 
5049 23 Oct 18 nicklas 5 #     compared to the normal reference samples are highlighted
5049 23 Oct 18 nicklas 6 #  2) Distribution of average mBAF values compared to a set of reference normal samples
5049 23 Oct 18 nicklas 7 #
5049 23 Oct 18 nicklas 8 #  Parameters
5049 23 Oct 18 nicklas 9 #  sampleName: Name of the sample (alignment)
5049 23 Oct 18 nicklas 10 #  mBafFile: Path to tab-separated file with mBAF data. The file should have the following
5049 23 Oct 18 nicklas 11 #            columns: Region, Count, AvgMBaf, PVal
5049 23 Oct 18 nicklas 12 #  datadir: Directory with reference data for the normal samples
5049 23 Oct 18 nicklas 13 #  plotFile: Name of file to create
5049 23 Oct 18 nicklas 14 #  width, height: Size of the plot
5049 23 Oct 18 nicklas 15 #  pval: p-value limit for regions that should be highlighted
5051 24 Oct 18 nicklas 16 #  minCount: If the number of SNPs in a region is below this value that bar is not colored or highlighted
5049 23 Oct 18 nicklas 17 #
5051 24 Oct 18 nicklas 18 mBafPlot = function(sampleName, mBafFile="mbafdata.txt", datadir=".", plotFile="mbafplot.png", width=1600, height=1000, pval=0.01, minCount=20)
5049 23 Oct 18 nicklas 19 {
5049 23 Oct 18 nicklas 20   # Load mBAF data
5049 23 Oct 18 nicklas 21   # Tab-separated with headers: Region, Count, AvgMBaf, PVal
5049 23 Oct 18 nicklas 22   mbafdata <- read.table(mBafFile, header=T, sep="\t", check.names=F, stringsAsFactors=F);
5049 23 Oct 18 nicklas 23   
5049 23 Oct 18 nicklas 24   # Load mBAF data for normal references
5049 23 Oct 18 nicklas 25   # Tab-separated with headers: Sample, 1p, 1q, ...
5049 23 Oct 18 nicklas 26   normalrefs <- read.table(paste(datadir, "normal-refs.txt", sep="/"), header=T, sep="\t", check.names=F, stringsAsFactors=F);
5049 23 Oct 18 nicklas 27   # We only need the mBAF values as one long list of values
5049 23 Oct 18 nicklas 28   normalrefs <- na.omit(unlist(normalrefs[,2:length(names(normalrefs))], use.names=F));
5049 23 Oct 18 nicklas 29   
5049 23 Oct 18 nicklas 30   png(plotFile, width=width, height=height);
5049 23 Oct 18 nicklas 31   
5049 23 Oct 18 nicklas 32   # Prepare for plottting 2 plots in a single window above each other
5049 23 Oct 18 nicklas 33   # and make room for adding title
5049 23 Oct 18 nicklas 34   par(mfrow=c(2, 1), oma=c(0,0,2,0), mar=c(4.5, 4.5, 2.0, 2.5));
5049 23 Oct 18 nicklas 35   
5049 23 Oct 18 nicklas 36   # *** Plot 1: Bar plot with Region on X axis and average mBAF on Y axis
5051 24 Oct 18 nicklas 37   # Regions with low p-value are highlighted (but only if we have enough SNPs)
5051 24 Oct 18 nicklas 38   barcolors <- barcolors <- ifelse(mbafdata$Count < minCount, "white", ifelse(mbafdata$PVal < pval, "darkblue", "gray"));
5051 24 Oct 18 nicklas 39   barborders <- ifelse(mbafdata$Count < minCount, "gray", "black");
5051 24 Oct 18 nicklas 40   bp <- barplot(mbafdata$AvgMBaf, main="Average mBAF per region", cex.main=1.25, font.main=1, ylab="mBAF", ylim=c(0.5,0.9), names.arg=mbafdata$Region, col=barcolors, border=barborders, xpd=F, xaxs = "i", xlim=c(0, 1.2*length(mbafdata$AvgMBaf)));
5049 23 Oct 18 nicklas 41   # Add Count (of SNPS making up the average) above each bar
5051 24 Oct 18 nicklas 42   text(x=bp, y=ifelse(is.na(mbafdata$AvgMBaf), 0.5, mbafdata$AvgMBaf), label=mbafdata$Count, cex=0.85, pos=3);
5049 23 Oct 18 nicklas 43   # Draw x-axis at y=0.5
5049 23 Oct 18 nicklas 44   axis(1, pos=0.5, at=c(0,1.2*length(mbafdata$AvgMBaf)), tick=T, labels=F);
5049 23 Oct 18 nicklas 45   
5049 23 Oct 18 nicklas 46   # *** Plot 2: Density plot of the average mBAF values
5051 24 Oct 18 nicklas 47   densitySample = density(na.omit(mbafdata$AvgMBaf));
5049 23 Oct 18 nicklas 48   densityNormal = density(normalrefs);
5049 23 Oct 18 nicklas 49   maxY <- max(densitySample$y, densityNormal$y);
5049 23 Oct 18 nicklas 50   plot(densitySample, main="Distribution of mBAF in sample vs. normals", cex.main=1.25, font.main=1, xlab="mBAF", ylab="", yaxt="n", bty="n", xaxs = "i", xlim=c(0.5, 0.9), ylim=c(0, maxY), lwd=2);
5049 23 Oct 18 nicklas 51   # Fill sample distribution with semi-transparant gray
5049 23 Oct 18 nicklas 52   polygon(densitySample, col=adjustcolor("lightgray", alpha.f=0.5));
5049 23 Oct 18 nicklas 53   # Draw normal reference distribution with dashed line
5049 23 Oct 18 nicklas 54   lines(densityNormal, col="red", lty=2, lwd=2);
5049 23 Oct 18 nicklas 55   # Mark individual mBAF values with small ticks on the x-axis
5049 23 Oct 18 nicklas 56   rug(mbafdata$AvgMBaf)
5049 23 Oct 18 nicklas 57   
5049 23 Oct 18 nicklas 58   # Add title
5049 23 Oct 18 nicklas 59   title(main=sampleName, outer=T, cex.main=1.75);
5049 23 Oct 18 nicklas 60   off <- dev.off();
5049 23 Oct 18 nicklas 61 }

Generated on Tue Nov 21 03:32:49 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)