Entering edit mode
stephen sefick
▴
60
@stephen-sefick-3922
Last seen 8.0 years ago
Hello:
I am having a problem with filterVcf. I am passing 3 filters defined
as fxns to FilterRules, and I am having the Error below. All of these
filters work individually, but when passed together in a list they do
not work. I have included the relevant code below. If it is helpful, I
can provide full script and test data. Please let me know.
Error:
starting filter filtering 10 records Error in value[[3L]](cond) : Filter 'b' failed: dim(X) must have a positive length Calls: filterVcf ... tryCatch -> tryCatchList -> tryCatchOne -> <Anonymous> Execution halted
Relevant Code:
#filters #snps #"QD < 2.0"; "MQ < 40.0"; "FS > 60.0"; "SOR > 3.0"; "MQRankSum < -12.5"; "ReadPosRankSum < -8.0"; "DP < 5"; "DP > 120" #indels #"FS > 200.0"; "ReadPosRankSum < -20.0"; "SOR > 10.0"; "DP < 5"; "DP > 120" ###################################################################################################################### #RGQ RGQ_filter <- function(x, RGQvalue=RGQ_value, minDP=min_DP, maxDP=max_DP) { ## Filter on RGQ > value; default is 30 #out <- apply(geno(x)$RGQ>=value, 1, function(x)sum(x, na.rm=TRUE))==4 out <- apply(geno(x)$RGQ >= RGQvalue & geno(x)$DP > minDP & geno(x)$DP < maxDP, 1, function(x)sum(x, na.rm=TRUE))==4 snps_indels <- isIndel(x) | isSNV(x) out. <- ifelse(snps_indels==TRUE, snps_indels, out) out.[is.na(out.)] <- TRUE return(out.) } #ADDED 20161111 #GATK HARD FILTERS #snps #"QD < 2.0"; "MQ < 40.0"; "FS > 60.0"; "SOR > 3.0"; "MQRankSum < -12.5"; "ReadPosRankSum < -8.0"; "DP < 5"; "DP > 120" snps <- function(x, minDP.=min_DP, maxDP.=max_DP, GQ.=GQ){ out. <- apply(info(x)$QD > 2 & info(x)$MQ > 40 & info(x)$FS < 60 & info(x)$SOR < 3 & info(x)$MQRankSum > -12.5 & info(x)$ReadPosRankSum > -8 & geno(x)$DP >= minDP. & geno(x)$DP <= maxDP. & geno(x)$GQ >= GQ., 1, function(x)sum(x, na.rm=TRUE))==4 out.. <- ifelse(!isSNV(x)==TRUE, TRUE, out.) out..[is.na(out..)] <- TRUE return(out..) } #"FS > 200.0"; "ReadPosRankSum < -20.0"; "SOR > 10.0"; "DP < 5"; "DP > 120" indels <- function(x, minDP..=min_DP, maxDP..=max_DP, GQ..=GQ){ out <- apply(info(x)$FS < 200 & info(x)$SOR < 3 & info(x)$ReadPosRankSum > -8 & geno(x)$DP>=minDP.. & geno(x)$DP<=maxDP.. & geno(x)$GQ>=GQ.., 1, function(x)sum(x, na.rm=TRUE))==4 out. <- ifelse(!isIndel(x)==TRUE, TRUE, out) out.[is.na(out.)] <- TRUE return(out.) } ############################################################## ############################################################## ############################################################## filters <- FilterRules(list(a=RGQ_filter, b=indels))