@@ -1,7 +1,7 @@

 Package: SICtools

 Type: Package

 Title: Find SNV/Indel differences between two bam files with near relationship

-Version: 1.1.0

+Version: 1.1.1

 Date: 2014-12-11

 Author: Xiaobin Xing, Wu Wei

 Maintainer: Xiaobin Xing <xiaobinxing0316@gmail.com>

@@ -15,7 +15,7 @@ Description: This package is to find SNV/Indel differences between two

 License: GPL (>=2)

 LazyLoad: Yes

 Depends: R (>= 3.0.0), methods, Rsamtools (>= 1.18.1), doParallel (>= 1.0.8), Biostrings (>= 2.32.1), stringr (>= 0.6.2), matrixStats

-               (>= 0.10.0)

+               (>= 0.10.0), plyr (>= 1.8.3)

 Imports: foreach, IRanges, GenomicRanges

 Suggests: knitr, RUnit, BiocGenerics

 biocViews: Alignment, Sequencing, Coverage, SequenceMatching, QualityControl, DataImport, Software, SNP, VariantDetection

@@ -6,5 +6,6 @@ import(Biostrings)

 import(matrixStats)

 import(foreach)

 import(IRanges)

+import(plyr)

 importFrom(GenomicRanges, GRanges)

 export(snpDiff,indelDiff)

@@ -73,7 +73,7 @@ function(ppIndel,rowIndex,refFsa,pValueCutOff,gtDistCutOff,verbose){

 		## select two most alts

 		if(ncol(gtMtx) >= 3){  

-			gtCountMtx <- gtMtx[,c('*',names(sort(colSums(gtMtx[,which(colnames(gtMtx)  != '*'),drop=FALSE]),

+			gtCountMtx <- gtMtx[,c('*',names(sort(colSums(gtMtx[,which(colnames(gtMtx)  != '*')]),

 											decreasing=TRUE)[1:2]))]     

 		}else if(ncol(gtMtx) == 2){

@@ -89,7 +89,7 @@ function(ppIndel,rowIndex,refFsa,pValueCutOff,gtDistCutOff,verbose){

 		gtPvalue <- ifelse(is.na(gtPvalue),1,gtPvalue)

 		gtDist <- ifelse(is.na(gtDist),0,gtDist)

-		if(gtPvalue < pValueCutOff & gtDist > gtDistCutOff){

+		if(gtPvalue <= pValueCutOff & gtDist >= gtDistCutOff){

 			## format genotype to output

 			delSeqLen <- as.numeric(str_extract(colnames(gtCountMtx),'-\\d+'))

@@ -113,21 +113,23 @@ function(ppIndel,rowIndex,refFsa,pValueCutOff,gtDistCutOff,verbose){

 			if(ncol(gtCountMtx) == 3){

 				## output table

-				tmpOut <- as.character(c(rowDat[,2:3],colnames(gtCountMtx),

+				tmpOut <- c(rowDat[,2:3],colnames(gtCountMtx),

 								gtCountMtx[1,1], gtCountMtx[1,2],gtCountMtx[1,3],

 								gtCountMtx[2,1], gtCountMtx[2,2],gtCountMtx[2,3],

-								pValue = format(gtPvalue,digits=3,scientific=TRUE),

-								gtDist = format(gtDist,digits=3,scientific=TRUE)))

-				

+								pValue = gtPvalue,

+								gtDist = gtDist)

+						

 			}else if(ncol(gtCountMtx) < 3){

 				## output table

-				tmpOut <- as.character(c(rowDat[,2:3],colnames(gtCountMtx),NA,

+				tmpOut <- c(rowDat[,2:3],colnames(gtCountMtx),NA,

 								gtCountMtx[1,1], gtCountMtx[1,2],NA,

 								gtCountMtx[2,1], gtCountMtx[2,2],NA,

-								format(gtPvalue,digits=3,scientific=TRUE),

-								format(gtDist,digits=3,scientific=TRUE)))

+								pValue = gtPvalue,

+								gtDist = gtDist)

 			}

+			

+			return(tmpOut)

 		}

 	}

 }

@@ -2,6 +2,29 @@ indelDiff <-

 function(bam1,bam2,refFsa,regChr,regStart,regEnd,minBaseQuality = 13,

 				minMapQuality = 0,nCores = 1,pValueCutOff = 0.05,gtDistCutOff = 0.1,verbose=TRUE){

+	## for test

+#	library(doParallel)

+#	library(Rsamtools)

+#	library(matrixStats)

+#	library(plyr)

+#	

+#	bam1 <- '/g/steinmetz/xing/SICtoolsDebug/20160319/debug/ER.sorted.bam'

+#	bam2 <- '/g/steinmetz/xing/SICtoolsDebug/20160319/debug/d1c.0.01.bam'

+#	refFsa <- '/g/steinmetz/xing/SICtoolsDebug/20160319/debug/S288c.fa'

+#	regChr <- 'chr2'

+#	regStart <- 1

+#	regEnd <- 230218

+#	minBaseQuality <- 13

+#	minMapQuality <- 0

+#	nCores <- 10

+#	pValueCutOff <- 0.05

+#	gtDistCutOff <- 0.1

+#	verbose <- TRUE

+#	

+	

+	## disable dimention dropping for matrix

+	`[` <- function(...) base::`[`(...,drop=FALSE)

+	

 	## set NULL

 	rowIndex <- NULL

@@ -11,12 +34,12 @@ function(bam1,bam2,refFsa,regChr,regStart,regEnd,minBaseQuality = 13,

 	if(!file.exists(refFsa)){stop("ERROR: input refFsa file doesn't exist!")}

 	if(regChr < 0){stop("ERROR: please input valid regChr!")}

-	if(as.numeric(regStart) < 0){stop("ERROR: regStart should be more than 0!")}

-	if(as.numeric(regEnd) < 0){stop("ERROR: regEnd should be more than 0!")}

-	if(as.numeric(regStart) > as.numeric(regEnd)){stop("ERROR: regStart should be smaller than regEnd!")}

+	regStart <- as.numeric(regStart);if(regStart < 0 ){stop("ERROR: regStart should be more than 0!")}

+	regEnd <- as.numeric(regEnd);if(regEnd < 0 ){stop("ERROR: regEnd should be more than 0!")}

+	if(regStart > regEnd){stop("ERROR: regStart should be smaller than regEnd!")}

-	regStart <- format(as.numeric(regStart),scientific=FALSE)

-	regEnd <- format(as.numeric(regEnd),scientific=FALSE)

+	regStart <- format(regStart,scientific=FALSE)

+	regEnd <- format(regEnd,scientific=FALSE)

 	minBaseQuality <- as.numeric(minBaseQuality)

 	minMapQuality <- as.numeric(minMapQuality)

 	pValueCutOff <- as.numeric(pValueCutOff)

@@ -31,46 +54,35 @@ function(bam1,bam2,refFsa,regChr,regStart,regEnd,minBaseQuality = 13,

 	ppIndel <- read.delim(pipe(paste(pathSICtools,' mpileup -Q ',minBaseQuality,' -q ',minMapQuality, ' -Ogf ',refFsa,' ',bam1,' ',bam2,' -r ',regChr,':',regStart,'-',regEnd,sep='')),

 			header=FALSE,colClasses = 'character',col.names=paste('V',1:12,sep=''))

-	## remove items with too low indel hits

-	ppIndel <- ppIndel[which(!(unlist(lapply(ppIndel$V6,

-												function(x){sum(gregexpr('[+-]',x)[[1]] > 0)}))/as.numeric(ppIndel$V5) < 0.05 

-								& unlist(lapply(ppIndel$V10,

-												function(x){sum(gregexpr('[+-]',x)[[1]] > 0)}))/as.numeric(ppIndel$V9) < 0.05) 

-							& as.numeric(ppIndel$V5) > 0 

-							& as.numeric(ppIndel$V9) > 0),,drop=FALSE]

+	## remove items with too low indel hits	

+	ppIndel <- ppIndel[which(!(unlist(lapply(gregexpr('[+-]',ppIndel$V6),function(x){sum(x > 0)}),use.names=FALSE)/as.numeric(ppIndel$V5) < 0.05 & 

+								

+	                           unlist(lapply(gregexpr('[+-]',ppIndel$V10),function(x){sum(x > 0)}),use.names=FALSE)/as.numeric(ppIndel$V9) < 0.05) & 

+						       as.numeric(ppIndel$V5) > 0 & as.numeric(ppIndel$V9) > 0),]

-	if(nrow(ppIndel) > 0){ 

+	if(nrow(ppIndel) > 0){

 		## test out

-		testOutIndel <- foreach(rowIndex=1:nrow(ppIndel)) %dopar% 

-				.indelDiffFunc(ppIndel,rowIndex,refFsa,pValueCutOff,gtDistCutOff,verbose)

-		

-		## remove NULL

-		testOutIndel <- testOutIndel[!sapply(testOutIndel,is.null)]

+		testOutIndel <- ldply(1:nrow(ppIndel),function(rowIndex){.indelDiffFunc(ppIndel,rowIndex,refFsa,pValueCutOff,gtDistCutOff,verbose)},.parallel=TRUE)

-		## combine the result

-		testOutIndel <- do.call(rbind,testOutIndel)

-		

-		if(!is.null(testOutIndel)){

+		if(nrow(testOutIndel) > 0){

-			testIndelMtx <- testOutIndel[,-c(1,3:5),drop=FALSE]

-			class(testIndelMtx) <- 'numeric'

-			testIndelDf <- data.frame(chr=testOutIndel[,1],ref=testOutIndel[,3],alt1=testOutIndel[,4],alt2 = testOutIndel[,5],testIndelMtx,stringsAsFactors=FALSE)

+			colnames(testOutIndel) <- c('chr','pos','ref','altGt1','altGt2','refBam1Count','altGt1Bam1Count',

+					'altGt2Bam1Count','refBam2Count','altGt1Bam2Count','altGt2Bam2Count','p.value','d.value')

+			class(testOutIndel$pos) <- 'numeric'

+			class(testOutIndel$refBam1Count) <- 'numeric'

+			class(testOutIndel$altGt1Bam1Count) <- 'numeric'

+			class(testOutIndel$altGt2Bam1Count) <- 'numeric'

+			class(testOutIndel$refBam2Count) <- 'numeric'

+			class(testOutIndel$altGt1Bam2Count) <- 'numeric'

+			class(testOutIndel$altGt2Bam2Count) <- 'numeric'

+			class(testOutIndel$p.value) <- 'numeric'

+			class(testOutIndel$d.value) <- 'numeric'

-			## at least one result

-			if(nrow(testIndelDf) > 0){

-				

-				rownames(testIndelDf) <- 1:nrow(testIndelDf)

-				testIndelDf <- testIndelDf[,c(1,5,2:4,6:13)]

-				colnames(testIndelDf) <- c('chr','pos','ref','altGt1','altGt2','refBam1Count','altGt1Bam1Count',

-						'altGt1Bam2Count','refBam2Count','altGt2Bam1Count','altGt2Bam2Count','p.value','d.value')

-				

-				testIndelDf <- testIndelDf[,c('chr','pos','ref','altGt1','altGt2','refBam1Count','altGt1Bam1Count','altGt2Bam1Count',

-								'refBam2Count','altGt1Bam2Count','altGt2Bam2Count','p.value','d.value')]

-				

-				return(testIndelDf[,,drop=FALSE])

-				

-			}

+			return(testOutIndel)

 		}

+		

+	}else{

+		NULL

 	}

 }

@@ -1,6 +1,29 @@

 snpDiff <-

-function(bam1,bam2,refFsa,regChr,regStart,regEnd,minBaseQuality = 13,

-		minMapQuality = 0,nCores = 1,pValueCutOff= 0.05,baseDistCutOff = 0.1,verbose=TRUE){

+		function(bam1,bam2,refFsa,regChr,regStart,regEnd,minBaseQuality = 13,

+				minMapQuality = 0,nCores = 1,pValueCutOff= 0.05,baseDistCutOff = 0.1,verbose=TRUE){

+	

+	

+	## for test

+#	library(doParallel)

+#	library(Rsamtools)

+#	library(matrixStats)

+#	library(plyr)

+#	

+#	bam1 <- '/g/steinmetz/xing/SICtoolsDebug/20160319/debug/ER.sorted.bam'

+#	bam2 <- '/g/steinmetz/xing/SICtoolsDebug/20160319/debug/d1c.0.01.bam'

+#	refFsa <- '/g/steinmetz/xing/SICtoolsDebug/20160319/debug/S288c.fa'

+#	regChr <- 'chr1'

+#	regStart <- 1

+#	regEnd <- 230218

+#	minBaseQuality <- 13

+#	minMapQuality <- 0

+#	nCores <- 10

+#	pValueCutOff <- 0.05

+#	baseDistCutOff <- 0.1

+#	verbose <- TRUE

+	

+	## disable dimention dropping for matrix

+	`[` <- function(...) base::`[`(...,drop=FALSE)

 	## validate the parameters

 	if(!file.exists(bam1)){stop("ERROR: input bam1 file doesn't exist!")}

@@ -8,9 +31,9 @@ function(bam1,bam2,refFsa,regChr,regStart,regEnd,minBaseQuality = 13,

 	if(!file.exists(refFsa)){stop("ERROR: input refFsa file doesn't exist!")}

 	if(regChr < 0 ){stop("ERROR: please input valid regChr!")}

-	if(as.numeric(regStart) < 0 ){stop("ERROR: regStart should be more than 0!")}

-	if(as.numeric(regEnd) < 0 ){stop("ERROR: regEnd should be more than 0!")}

-	if(as.numeric(regStart) > as.numeric(regEnd)){stop("ERROR: regStart should be smaller than regEnd!")}

+	regStart <- as.numeric(regStart);if(regStart < 0 ){stop("ERROR: regStart should be more than 0!")}

+	regEnd <- as.numeric(regEnd);if(regEnd < 0 ){stop("ERROR: regEnd should be more than 0!")}

+	if(regStart > regEnd){stop("ERROR: regStart should be smaller than regEnd!")}

 	pValueCutOff <- as.numeric(pValueCutOff)

 	baseDistCutOff <- as.numeric(baseDistCutOff) 

@@ -27,30 +50,24 @@ function(bam1,bam2,refFsa,regChr,regStart,regEnd,minBaseQuality = 13,

 		countMtx <- cbind(countMtx,x[['pos']])

 		## duplicated index

-		dupIndex <- duplicated(countMtx[,-c(5,10,11),drop=FALSE])

+		dupIndex <- duplicated(countMtx[,-c(5,10,11)])

 		## max bam1 == max bam2 index

-		maxBaseIndex <- rowMaxs(countMtx[,1:4,drop=FALSE])/rowSums(countMtx[,1:4,drop=FALSE]) > 0.95 & 

-				rowMaxs(countMtx[,6:9,drop=FALSE])/rowSums(countMtx[,6:9,drop=FALSE]) > 0.95 & 

-				max.col(countMtx[,1:4,drop=FALSE]) == max.col(countMtx[,6:9,drop=FALSE])

-		

-#		maxBaseIndex <- rowMaxs(countMtx[,,drop=FALSE],cols=1:4)/rowSums(countMtx[,,drop=FALSE],cols=1:4) > 0.95 & 

-#				rowMaxs(countMtx[,,drop=FALSE],cols=6:9)/rowSums(countMtx[,,drop=FALSE],cols=6:9) > 0.95 & 

-#				colMaxs(countMtx[,,drop=FALSE],cols=1:4) == colMaxs(countMtx[,,drop=FALSE],cols=6:9)

+		maxBaseIndex <- rowMaxs(countMtx[,1:4])/rowSums(countMtx[,1:4]) > 0.95 & 

+				rowMaxs(countMtx[,6:9])/rowSums(countMtx[,6:9]) > 0.95 & 

+				max.col(countMtx[,1:4]) == max.col(countMtx[,6:9])

 		## all 0 test

-		allZeroIndex <- rowSums(countMtx[,1:4,drop=FALSE]) == 0 | rowSums(countMtx[,6:9,drop=FALSE]) == 0

-#		allZeroIndex <- rowSums(countMtx[,,drop=FALSE],cols=1:4) == 0 | rowSums(countMtx[,,drop=FALSE],cols=6:9) == 0

-		

+		allZeroIndex <- rowSums(countMtx[,1:4]) == 0 | rowSums(countMtx[,6:9]) == 0		

 		## test index

 		testIndex <- !dupIndex & !maxBaseIndex & !allZeroIndex

 		if(sum(testIndex) > 0){ ## at least one candidate

-			countMtxTest <- countMtx[testIndex,,drop=FALSE]

+			countMtxTest <- countMtx[testIndex,]

-			testTmp <- lapply(1:nrow(countMtxTest),function(rowIndex){

+			testResultMtx <- do.call(rbind,llply(1:nrow(countMtxTest),function(rowIndex){

 						rowInfo <- countMtxTest[rowIndex,]

@@ -58,7 +75,7 @@ function(bam1,bam2,refFsa,regChr,regStart,regEnd,minBaseQuality = 13,

 						testMtx <- matrix(rowInfo[c(1:4,6:9)],nrow=2,byrow=TRUE)

 						## remove 0 columns

-						testMtx <- testMtx[,colSums(testMtx) > 0,drop=FALSE]

+						testMtx <- testMtx[,colSums(testMtx) > 0]

 						## get pValue and dist

 						pValue <- fisher.test(testMtx)$p.value

@@ -69,39 +86,24 @@ function(bam1,bam2,refFsa,regChr,regStart,regEnd,minBaseQuality = 13,

 						if(pValue <= as.numeric(pValueCutOff) & baseDist >= as.numeric(baseDistCutOff)){

-							as.character(c(names(x[['seqnames']]),rowInfo[11],

-											rowInfo[1],rowInfo[2],rowInfo[3],rowInfo[4],rowInfo[5],

-											rowInfo[6],rowInfo[7],rowInfo[8],rowInfo[9],rowInfo[10],

-											format(pValue,digits=3,scientific=TRUE),

-											format(baseDist,digits=3,scientific=TRUE)))

+							c(rowInfo[1],rowInfo[2],rowInfo[3],rowInfo[4],rowInfo[5],

+											rowInfo[6],rowInfo[7],rowInfo[8],rowInfo[9],rowInfo[10],rowInfo[11],

+											pValue,baseDist)

 						}

-					})

-			

-			testTmp2 <- do.call(rbind,testTmp[,drop=FALSE])

-			testMtx <- testTmp2[,-1,drop=FALSE]

-			class(testMtx) <- 'numeric'

-			testDf <- data.frame(chr = testTmp2[,1],testMtx,stringsAsFactors=FALSE)

-			

+					}))

+						

 			## for genotype duplicated positions, trace back

-			countMtxDup <- countMtx[!maxBaseIndex & !allZeroIndex,,drop=FALSE]

-			countMtxDupDf <- as.data.frame(countMtxDup[duplicated(countMtxDup[,-11,drop=FALSE]),,drop=FALSE],stringsAsFactors=FALSE)

+			countMtxDup <- countMtx[!maxBaseIndex & !allZeroIndex,]

+			countMtxDupIndex <- duplicated(countMtxDup[,-11])

-			if(nrow(countMtxDupDf) > 0 & nrow(testDf) > 0){

-				

-				countMtxDupDfFull <- cbind(countMtxDupDf,testDf[match(apply(countMtxDupDf[,-11,drop=FALSE],1,paste,collapse=' '),apply(testDf[,3:12],1,paste,collapse=' ')),c('X12','X13')])

-				countMtxDupDfFull <- countMtxDupDfFull[!is.na(countMtxDupDfFull$X12),,drop=FALSE]

-				

-				if(nrow(countMtxDupDfFull) > 0){

-					countMtxDupDfFull[,'X14'] <- unique(testDf$chr)

-					countMtxDupDfFull <- countMtxDupDfFull[,c(14,11,1:10,12,13)]

-					colnames(countMtxDupDfFull) <- colnames(testDf)

-					testDf <- rbind(testDf,countMtxDupDfFull)

-					testDf <- testDf[order(testDf[,'X1']),]

-				}

+			if(sum(countMtxDupIndex) > 0 & !is.null(testResultMtx)){

+				countMtxDupFull <- cbind(countMtxDup,testResultMtx[match(apply(countMtxDup[,1:10],1,paste,collapse='_'),apply(testResultMtx[,1:10],1,paste,collapse='_')),c(12,13)])

+				countMtxDupFull <- countMtxDupFull[!is.na(countMtxDupFull[,12]),]

+				testResultMtx <- unique(rbind(testResultMtx,countMtxDupFull))

 			}

-			return(testDf)          

+			return(testResultMtx)          

 		}

 	}

@@ -111,12 +113,8 @@ function(bam1,bam2,refFsa,regChr,regStart,regEnd,minBaseQuality = 13,

 	## build bam files for pileup

 	ppFiles <- PileupFiles(c(bam1,bam2))

-	## get region info

-	regStart <- as.numeric(regStart)

-	regEnd <- as.numeric(regEnd) + 1

-	

 	## split into 10M block

-	regSplit <- shift(as(breakInChunks(regEnd - regStart,1e7L),'IRanges'),regStart-1)

+	regSplit <- shift(as(breakInChunks(regEnd - regStart + 1,1e7L),'IRanges'),regStart-1)

 	## for each region

 	regSplitDiffFunc <- function(splitIndex){

@@ -127,7 +125,7 @@ function(bam1,bam2,refFsa,regChr,regStart,regEnd,minBaseQuality = 13,

 		}

 		regSplitParam <- ApplyPileupsParam(flag = scanBamFlag(isDuplicate=FALSE),

-				which=GRanges(regChr, IRanges(start=start(regSplit)[splitIndex], end=end(regSplit)[splitIndex])), 

+				which=GRanges(regChr, IRanges(start=start(regSplit[splitIndex]), end=end(regSplit[splitIndex]))), 

 				what='seq',minDepth = 0L,minBaseQuality, minMapQuality,maxDepth = 2500L, yieldBy = 'range')

 		applyPileups(ppFiles,calcInfoRange,param=regSplitParam)[[1]]

@@ -138,23 +136,21 @@ function(bam1,bam2,refFsa,regChr,regStart,regEnd,minBaseQuality = 13,

 	splitIndex <- NULL

 	## range out

-	rangeOut <- foreach(splitIndex=1:length(regSplit)) %dopar% regSplitDiffFunc(splitIndex)

-	

-	## remove NULL

-	rangeOut <- rangeOut[!sapply(rangeOut,is.null)]

-	

-	## combine the results

-	rangeOut <- do.call(rbind,rangeOut)

-	

-	if(!is.null(rangeOut)){

+	rangeOut <- ldply(1:length(regSplit),regSplitDiffFunc,.parallel=TRUE)

+		

 		if(nrow(rangeOut) > 0){

+			colnames(rangeOut) <- c('A1','C1','G1','T1','N1','A2','C2','G2','T2','N2','pos','p.value','d.value')

+			rangeOut$chr <- regChr

+			

 			## get ref by Biostrings

-			rangeOut$ref <- as.character(getSeq(FaFile(refFsa),GRanges(seqnames=rangeOut[,1],IRanges(start=rangeOut[,2],width=1))))

+			rangeOut$ref <- as.character(getSeq(FaFile(refFsa),GRanges(seqnames=rangeOut$chr,IRanges(start=rangeOut$pos,width=1))))

+			rangeOut <- rangeOut[order(rangeOut$pos),]

 			## format output    

-			rangeOut <- rangeOut[,c(1,2,15,3:14),drop=FALSE]

-			colnames(rangeOut) <- c('chr','pos','ref','A1','C1','G1','T1','N1','A2','C2','G2','T2','N2','p.value','d.value')

-			return(rangeOut[,,drop=FALSE])

-		}}

+			rangeOut <- rangeOut[,c('chr','pos','ref','A1','C1','G1','T1','N1','A2','C2','G2','T2','N2','p.value','d.value')]

+			return(rangeOut)

+		}else{

+			NULL

+		}

 }