added minIndelCount (short: minCnt) command line argument. The call is made only if the number of reads supporting the consensus indel is equal or greater than the specified value (default: 0, so only minFraction filter is on in default runs!)

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@1320 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
asivache 2009-07-27 15:22:51 +00:00
parent 73ddf21bb7
commit 2499c09256
1 changed files with 14 additions and 6 deletions

View File

@ -45,12 +45,17 @@ public class IndelGenotyperWalker extends ReadWalker<Integer,Integer> {
doc="used only with --somatic; normal sample must have at least minNormalCoverage or more reads to call germline/somatic indel", required=false)
int minNormalCoverage = 4;
@Argument(fullName="minFraction", shortName="minFraction",
doc="minimum fraction of reads with CONSENSUS indel at a site, out of all reads covering the site, required for a consensus call"+
doc="Minimum fraction of reads with CONSENSUS indel at a site, out of all reads covering the site, required for a consensus call"+
" (fraction of non-consensus indels at the site is not considered here, see minConsensusFraction)", required=false)
double minFraction = 0.3;
@Argument(fullName="minConsensusFraction", shortName="minConsensusFraction",
doc="Minimum fraction of CONSENSUS indel observations at a site wrt all indel observations at the site required to make the call", required=false)
double minConsensusFraction = 0.7;
@Argument(fullName="minIndelCount", shortName="minCnt",
doc="Minimum count of reads supporting consensus indel required for making the call. "+
" This filter supercedes minFraction, i.e. indels with acceptable minFraction at low coverage "+
"(minIndelCount not met) will not pass.", required=false)
int minIndelCount = 0;
@Argument(fullName="refseq", shortName="refseq",
doc="Name of RefSeq transcript annotation file. If specified, indels will be annotated as GENOMIC/UTR/INTRON/CODING", required=false)
String RefseqFileName = null;
@ -275,13 +280,16 @@ public class IndelGenotyperWalker extends ReadWalker<Integer,Integer> {
}
/** Returns true if consensus (specified by the pair) should be considered a call given current values
* of the fraction cutoffs
* @param p
* @param coverage
* of the cutoffs.
* @param p pair with first element being the consensus indel variant, the second element being the <i>total</i> (consensus+others)
* count of indels at the site.
* @param coverage total coverage (number of spanning reads, including those with indel(s)) at the site.
* @return
*/
private boolean isCall(Pair<IndelVariant,Integer> p, int coverage) {
return ( (double)p.first.getCount() > minFraction * coverage && (double) p.first.getCount() > minConsensusFraction*p.second );
return ( p.first.getCount() >= minIndelCount &&
(double)p.first.getCount() > minFraction * coverage &&
(double) p.first.getCount() > minConsensusFraction*p.second );
}
/** Build output line for bed file and write it to the specified output writer if the latter is not null;
@ -502,7 +510,7 @@ public class IndelGenotyperWalker extends ReadWalker<Integer,Integer> {
}
message += "\tSOMATIC\t0/"+normal_cov;
} else {
Pair<IndelVariant,Integer> p_normal = findConsensus(tumor_variants);
Pair<IndelVariant,Integer> p_normal = findConsensus(normal_variants);
message += "\tGERMLINE\t"+p_normal.second+"/"+normal_cov;
}