[Patchwork-commits] r136 - pkg/patchwork/R

[noreply at r-forge.r-project.org]

Author: mayrhofer
Date: 2012-06-25 14:51:00 +0200 (Mon, 25 Jun 2012)
New Revision: 136

Modified:
   pkg/patchwork/R/patchwork.copynumbers.r
Log:
Updated CN calculation @ X and Y.

Modified: pkg/patchwork/R/patchwork.copynumbers.r
===================================================================
--- pkg/patchwork/R/patchwork.copynumbers.r	2012-06-25 12:49:40 UTC (rev 135)
+++ pkg/patchwork/R/patchwork.copynumbers.r	2012-06-25 12:51:00 UTC (rev 136)
@@ -1,4 +1,4 @@
-patchwork.copynumbers = function(name = "copynumbers_",cn2,delta,het,hom,maxCn=8,ceiling=1,forcedelta=F)
+patchwork.copynumbers = function(name = "copynumbers_",cn2,delta,het,hom,maxCn=8,ceiling=1,forcedelta=F,male.sample=F,male2femref=F)
 	{
 	
 	#packagepath = system.file(package="patchwork")
@@ -264,19 +264,50 @@
     			fullCN[i] <- paste('cn',Cn[i],'m',mCn[i],sep='') # Full description
   		}
  
+
     regions$Cn <- Cn
     regions$mCn <- mCn
     regions$fullCN <- fullCN
 	
-	meanCn <- weightedMean(regions$Cn,regions$np)
-	ix1 <- regions$Cn==trunc(meanCn)
-	ix2 <- regions$Cn==ceiling(meanCn)
-	xdelta <- weightedMean(regions$median[ix2],regions$np[ix2]) - weightedMean(regions$median[ix1],regions$np[ix1])
+    temp <- regions[is.autosome(regions$chr),]
+	meanCn <- weightedMean(temp$Cn,temp$np)
+	ix1 <- temp$Cn==trunc(meanCn)
+	ix2 <- temp$Cn==ceiling(meanCn)
+	xdelta <- weightedMean(temp$median[ix2],temp$np[ix2]) - weightedMean(temp$median[ix1],temp$np[ix1])
 	expected_delta <- 1/meanCn
 	
 	tumor_percentDNA <- xdelta / expected_delta
     tumor_percent <- tumor_percentDNA/(meanCn/2) / ( tumor_percentDNA/(meanCn/2) + (1-tumor_percentDNA) )
 
+    ## Fix sex chromosomes in case of male sample.
+    if (male.sample) {
+        ix <- which(!is.autosome(regions$chr))
+        regions$mCn[ix] <- 0
+        ## Treats differently depending on sex of the reference
+        if (!male2femref) {
+            # X and Y (originally cn1) were compared to cn1 (matched normal or male references).
+            # Unchanged copy number (=1) is expected at relative coverage 1. 
+            # Gains/loss are expected at xdelta (the delta of the autosomes)*2.
+            rawcopychange <- (regions$median[ix]-1) / (xdelta*2)
+            regions$Cn[ix] <- round(1+rawcopychange)
+            regions$fullCN[ix] <- paste('cn',round(1+rawcopychange),'m0')
+        } else 
+        {
+            # X and Y (originally cn1) were compared to cn2 (female references).
+            # Unchanged copy number (=1) is expected at relative coverage 1/2. 
+            # Gains/loss are expected at xdelta from there.
+            rawcopychange <- (regions$median[ix]-1/2) / (xdelta)
+            regions$Cn[ix] <- round(1+rawcopychange)
+            regions$fullCN[ix] <- paste('cn',round(1+rawcopychange),'m0')
+        }
+    }
+    ## And a small fix on female samples
+    if (!male.sample) {
+        ix <- which(regions$chr=='chrY')
+        regions$mCn[ix] <- regions$Cn[ix] <- 0
+    }
+
+
     regions$meanCn <- meanCn
     regions$tumor_percent <- tumor_percent
     write.csv(regions,file='Copynumbers.csv')