Entering edit mode
stephen sefick
▴
60
@stephen-sefick-3922
Last seen 9.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))
