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

 Package: SICtools

 Type: Package

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

-        relationship

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

 Version: 1.1.0

 Date: 2014-12-11

 Author: Xiaobin Xing, Wu Wei

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

         span no less than 2bp on both sides of indel region.

 License: GPL (>=2)

 LazyLoad: Yes

-Depends: R (>= 3.0.0), methods, Rsamtools (>= 1.18.1), doMC (>= 1.3.3),

-        Biostrings (>= 2.32.1), stringr (>= 0.6.2), matrixStats (>=

-        0.10.0)

+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)

 Imports: foreach, IRanges, GenomicRanges

 Suggests: knitr, RUnit, BiocGenerics

-biocViews: Alignment, Sequencing, Coverage, SequenceMatching,

-        QualityControl, DataImport, Software, SNP, VariantDetection

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

 VignetteBuilder: knitr

+Packaged: 2016-01-04 15:11:21 UTC; xing

@@ -1,5 +1,5 @@

 import(Rsamtools)

-import(doMC)

+import(doParallel)

 import(methods)

 import(stringr)

 import(Biostrings)

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

 	gtDistCutOff <- as.numeric(gtDistCutOff)

 	## number of cores used

-	registerDoMC(cores=nCores)

+	registerDoParallel(cores=nCores)

 	## get mpileupPlus result; the samtools will calculate BAQ for SNP around indel

 	pathSICtools <- system.file(package = "SICtools","etc","samtools2SIC")

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

 	minMapQuality <- as.numeric(minMapQuality)

 	## number of cores used

-	registerDoMC(cores=nCores)

+	registerDoParallel(cores=nCores)

 	################################ function to call difference ##############################

 	calcInfoRange <- function(x){

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

 				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)

+		

 		## 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

+		

 		## test index

 		testIndex <- !dupIndex & !maxBaseIndex & !allZeroIndex

@@ -75,6 +81,22 @@ function(bam1,bam2,refFsa,regChr,regStart,regEnd,minBaseQuality = 13,

 			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,]

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

+			

+			if(nrow(countMtxDupDf) > 0){

+				

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

+				countMtxDupDfFull <- countMtxDupDfFull[!is.na(countMtxDupDfFull$X12),]

+				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']),]

+			}

+			

 			return(testDf)          

 		}

 	}