... | ... |
@@ -1243,153 +1243,8 @@ construct.Illumina <- function(sampleSheet=NULL, |
1243 | 1243 |
return(cnSet) |
1244 | 1244 |
} |
1245 | 1245 |
|
1246 |
-genotype.Illumina <- function(sampleSheet=NULL, |
|
1247 |
- arrayNames=NULL, |
|
1248 |
- ids=NULL, |
|
1249 |
- path=".", |
|
1250 |
- arrayInfoColNames=list(barcode="SentrixBarcode_A", position="SentrixPosition_A"), |
|
1251 |
- highDensity=FALSE, |
|
1252 |
- sep="_", |
|
1253 |
- fileExt=list(green="Grn.idat", |
|
1254 |
- red="Red.idat"), |
|
1255 |
- cdfName, |
|
1256 |
- copynumber=TRUE, |
|
1257 |
- batch, |
|
1258 |
- fns, |
|
1259 |
- saveDate=TRUE, |
|
1260 |
- stripNorm=TRUE, |
|
1261 |
- useTarget=TRUE, |
|
1262 |
- mixtureSampleSize=10^5, |
|
1263 |
- eps=0.1, |
|
1264 |
- verbose=TRUE, |
|
1265 |
- seed=1, |
|
1266 |
- sns, |
|
1267 |
- probs=rep(1/3, 3), |
|
1268 |
- DF=6, |
|
1269 |
- SNRMin=5, |
|
1270 |
- recallMin=10, |
|
1271 |
- recallRegMin=1000, |
|
1272 |
- gender=NULL, |
|
1273 |
- returnParams=TRUE, |
|
1274 |
- badSNP=0.7) { |
|
1275 |
- is.lds <- ifelse(isPackageLoaded("ff"), TRUE, FALSE) |
|
1276 |
- if(missing(cdfName)) stop("must specify cdfName") |
|
1277 |
- if(!isValidCdfName(cdfName)) stop("cdfName not valid. see validCdfNames") |
|
1278 |
- callSet <- construct.Illumina(sampleSheet=sampleSheet, arrayNames=arrayNames, |
|
1279 |
- ids=ids, path=path, arrayInfoColNames=arrayInfoColNames, |
|
1280 |
- highDensity=highDensity, sep=sep, fileExt=fileExt, |
|
1281 |
- cdfName=cdfName, copynumber=copynumber, verbose=verbose, batch=batch, |
|
1282 |
- fns=fns, saveDate=saveDate) |
|
1283 |
- open(callSet) |
|
1284 |
- mixtureParams <- matrix(NA, 4, nrow(callSet)) |
|
1285 |
- is.snp <- isSnp(callSet) |
|
1286 |
- snp.index <- which(is.snp) |
|
1287 |
- |
|
1288 |
- RG = readIdatFiles(sampleSheet=sampleSheet, arrayNames=arrayNames, |
|
1289 |
- ids=ids, path=path, arrayInfoColNames=arrayInfoColNames, |
|
1290 |
- highDensity=highDensity, sep=sep, fileExt=fileExt, saveDate=saveDate) |
|
1291 |
- |
|
1292 |
- XY = RGtoXY(RG, chipType=cdfName) |
|
1293 |
- rm(RG); gc() |
|
1294 |
- if (missing(sns)) { sns = sampleNames(XY) #subsns = sampleNames(XY) |
|
1295 |
- } # else subsns = sns[j] |
|
1296 |
- |
|
1297 |
- res = preprocessInfinium2(XY, mixtureSampleSize=mixtureSampleSize, fitMixture=TRUE, verbose=verbose, |
|
1298 |
- seed=seed, eps=eps, cdfName=cdfName, sns=sns, stripNorm=stripNorm, useTarget=useTarget) |
|
1299 |
- # save.it=save.it, snpFile=snpFile, cnFile=cnFile) |
|
1300 |
- rm(XY); gc() |
|
1301 |
- if(verbose) message("Finished preprocessing.") |
|
1302 |
- np.index <- which(!is.snp) |
|
1303 |
- if(is.lds){ |
|
1304 |
- open(res[["A"]]) |
|
1305 |
- open(res[["B"]]) |
|
1306 |
- open(res[["SNR"]]) |
|
1307 |
- open(res[["mixtureParams"]]) |
|
1308 |
- bb = ocProbesets()*ncol(A)*8 |
|
1309 |
- ffrowapply(A(callSet)[i1:i2, ] <- res[["A"]][i1:i2, ], X=res[["A"]], BATCHBYTES=bb) |
|
1310 |
- ffrowapply(B(callSet)[i1:i2, ] <- res[["B"]][i1:i2, ], X=res[["B"]], BATCHBYTES=bb) |
|
1311 |
- if(length(np.index)>0) { |
|
1312 |
- for (j in 1:ncol(callSet)) { |
|
1313 |
- A(callSet)[np.index, j] <- res[["cnAB"]]$A[,j] |
|
1314 |
- B(callSet)[np.index, j] <- res[["cnAB"]]$B[,j] |
|
1315 |
- } |
|
1316 |
- } |
|
1317 |
- |
|
1318 |
- } else{ |
|
1319 |
- A(callSet)[snp.index, ] <- res[["A"]] |
|
1320 |
- B(callSet)[snp.index, ] <- res[["B"]] |
|
1321 |
- if(length(np.index)>0) { |
|
1322 |
- for (j in 1:ncol(callSet)) { |
|
1323 |
- A(callSet)[np.index, ] <- res[["cnAB"]]$A |
|
1324 |
- B(callSet)[np.index, ] <- res[["cnAB"]]$B |
|
1325 |
- } |
|
1326 |
- } |
|
1327 |
- } |
|
1328 |
- pData(callSet)$SKW <- res[["SKW"]] |
|
1329 |
- pData(callSet)$SNR <- res[["SNR"]] |
|
1330 |
- mixtureParams <- res[["mixtureParams"]] |
|
1331 |
- |
|
1332 |
-# if(verbose) message("Normalizing nonpolymorphic markers") |
|
1333 |
-# FUN <- ifelse(is.lds, "cnrma2", "cnrma") |
|
1334 |
- ## main purpose is to update 'alleleA' |
|
1335 |
-# cnrmaFxn <- function(FUN,...){ |
|
1336 |
-# switch(FUN, |
|
1337 |
-# cnrma=cnrma(...), |
|
1338 |
-# cnrma2=cnrma2(...)) |
|
1339 |
-# } |
|
1340 |
- ## consider passing only A for NPs. |
|
1341 |
-# AA <- cnrmaFxn(FUN, A=A(callSet), |
|
1342 |
-# filenames=filenames, |
|
1343 |
-# row.names=featureNames(callSet)[np.index], |
|
1344 |
-# cdfName=cdfName, |
|
1345 |
-# sns=sns, |
|
1346 |
-# seed=seed, |
|
1347 |
-# verbose=verbose) |
|
1348 |
-# if(!is.lds) A(callSet) <- AA |
|
1349 |
-# rm(AA) |
|
1350 |
- FUN <- ifelse(is.lds, "crlmmGT2", "crlmmGT") |
|
1351 |
- ## genotyping |
|
1352 |
- crlmmGTfxn <- function(FUN,...){ |
|
1353 |
- switch(FUN, |
|
1354 |
- crlmmGT2=crlmmGT2(...), |
|
1355 |
- crlmmGT=crlmmGT(...)) |
|
1356 |
- } |
|
1357 |
- tmp <- crlmmGTfxn(FUN, |
|
1358 |
- A=res[["A"]], |
|
1359 |
- B=res[["B"]], |
|
1360 |
- SNR=res[["SNR"]], |
|
1361 |
- mixtureParams=res[["mixtureParams"]], |
|
1362 |
- cdfName=cdfName, |
|
1363 |
- row.names=NULL, |
|
1364 |
- col.names=sampleNames(callSet), |
|
1365 |
- probs=probs, |
|
1366 |
- DF=DF, |
|
1367 |
- SNRMin=SNRMin, |
|
1368 |
- recallMin=recallMin, |
|
1369 |
- recallRegMin=recallRegMin, |
|
1370 |
- gender=gender, |
|
1371 |
- verbose=verbose, |
|
1372 |
- returnParams=returnParams, |
|
1373 |
- badSNP=badSNP) |
|
1374 |
- rm(res); gc() |
|
1375 |
- if(verbose) message("Genotyping finished. Updating container with genotype calls and confidence scores.") |
|
1376 |
- if(is.lds){ |
|
1377 |
- open(tmp[["calls"]]) |
|
1378 |
- open(tmp[["confs"]]) |
|
1379 |
- ffrowapply(snpCall(callSet)[i1:i2, ] <- tmp[["calls"]][i1:i2, ], X=tmp[["calls"]], BATCHBYTES=bb) |
|
1380 |
- ffrowapply(snpCallProbability(callSet)[i1:i2, ] <- tmp[["confs"]][i1:i2, ], X=tmp[["confs"]], BATCHBYTES=bb) |
|
1381 |
- close(tmp[["calls"]]) |
|
1382 |
- close(tmp[["confs"]]) |
|
1383 |
- } else { |
|
1384 |
- calls(callSet)[snp.index, ] <- tmp[["calls"]] |
|
1385 |
- snpCallProbability(callSet)[snp.index, ] <- tmp[["confs"]] |
|
1386 |
- } |
|
1387 |
- callSet$gender <- tmp$gender |
|
1388 |
- close(callSet) |
|
1389 |
- return(callSet) |
|
1390 |
-} |
|
1391 | 1246 |
|
1392 |
-genotype.Illumina2 <- function(sampleSheet=NULL, |
|
1247 |
+genotype.Illumina <- function(sampleSheet=NULL, |
|
1393 | 1248 |
arrayNames=NULL, |
1394 | 1249 |
ids=NULL, |
1395 | 1250 |
path=".", |