+ * After counting covariates in either the initial BAM File or again in the recalibrated BAM File, an analysis tool is available which + * reads the .csv file and outputs several PDF (and .dat) files for each read group in the given BAM. These PDF files graphically + * show the various metrics and characteristics of the reported quality scores (often in relation to the empirical qualities). + * In order to show that any biases in the reported quality scores have been generally fixed through recalibration one should run + * CountCovariates again on a bam file produced by TableRecalibration. In this way users can compare the analysis plots generated + * by pre-recalibration and post-recalibration .csv files. Our usual chain of commands that we use to generate plots of residual + * error is: CountCovariates, TableRecalibrate, samtools index on the recalibrated bam file, CountCovariates again on the recalibrated + * bam file, and then AnalyzeCovariates on both the before and after recal_data.csv files to see the improvement in recalibration. + * + *
+ * The color coding along with the RMSE is included in the plots to give some indication of the number of observations that went into + * each of the quality score estimates. It is defined as follows for N, the number of observations: + * + *
+ * NOTE: For those running this tool externally from the Broad, it is crucial to note that both the -Rscript and -resources options + * must be changed from the default. -Rscript needs to point to your installation of Rscript (this is the scripting version of R, + * not the interactive version) while -resources needs to point to the folder holding the R scripts that are used. For those using + * this tool as part of the Binary Distribution the -resources should point to the resources folder that is part of the tarball. + * For those using this tool by building from the git repository the -resources should point to the R/ subdirectory of the Sting checkout. + * + *
+ * See the GATK wiki for a tutorial and example recalibration accuracy plots. + * http://www.broadinstitute.org/gsa/wiki/index.php/Base_quality_score_recalibration + * + *
+ * The recalibration table file in CSV format that was generated by the CountCovariates walker. + *
+ * + *+ * java -Xmx4g -jar AnalyzeCovariates.jar \ + * -recalFile /path/to/recal.table.csv \ + * -outputDir /path/to/output_dir/ \ + * -resources resources/ \ + * -ignoreQ 5 + *+ * */ +@DocumentedGATKFeature( + groupName = "AnalyzeCovariates", + summary = "Package to plot residual accuracy versus error covariates for the base quality score recalibrator") public class AnalyzeCovariates extends CommandLineProgram { ///////////////////////////// // Command Line Arguments ///////////////////////////// - + /** + * After the header, data records occur one per line until the end of the file. The first several items on a line are the + * values of the individual covariates and will change depending on which covariates were specified at runtime. The last + * three items are the data- that is, number of observations for this combination of covariates, number of reference mismatches, + * and the raw empirical quality score calculated by phred-scaling the mismatch rate. + */ @Input(fullName = "recal_file", shortName = "recalFile", doc = "The input recal csv file to analyze", required = false) private String RECAL_FILE = "output.recal_data.csv"; @Argument(fullName = "output_dir", shortName = "outputDir", doc = "The directory in which to output all the plots and intermediate data files", required = false) @@ -67,11 +120,20 @@ public class AnalyzeCovariates extends CommandLineProgram { private int IGNORE_QSCORES_LESS_THAN = 5; @Argument(fullName = "numRG", shortName = "numRG", doc = "Only process N read groups. Default value: -1 (process all read groups)", required = false) private int NUM_READ_GROUPS_TO_PROCESS = -1; // -1 means process all read groups + + /** + * Combinations of covariates in which there are zero mismatches technically have infinite quality. We get around this situation + * by capping at the specified value. We've found that Q40 is too low when using a more completely database of known variation like dbSNP build 132 or later. + */ @Argument(fullName="max_quality_score", shortName="maxQ", required = false, doc="The integer value at which to cap the quality scores, default is 50") private int MAX_QUALITY_SCORE = 50; + + /** + * This argument is useful for comparing before/after plots and you want the axes to match each other. + */ @Argument(fullName="max_histogram_value", shortName="maxHist", required = false, doc="If supplied, this value will be the max value of the histogram plots") private int MAX_HISTOGRAM_VALUE = 0; - @Argument(fullName="do_indel_quality", shortName="indels", required = false, doc="If supplied, this value will be the max value of the histogram plots") + @Argument(fullName="do_indel_quality", shortName="indels", required = false, doc="If supplied, do indel quality plotting") private boolean DO_INDEL_QUALITY = false; diff --git a/public/java/src/org/broadinstitute/sting/analyzecovariates/package-info.java b/public/java/src/org/broadinstitute/sting/analyzecovariates/package-info.java new file mode 100644 index 000000000..9350e4a66 --- /dev/null +++ b/public/java/src/org/broadinstitute/sting/analyzecovariates/package-info.java @@ -0,0 +1,4 @@ +/** + * Package to plot residual accuracy versus error covariates for the base quality score recalibrator. + */ +package org.broadinstitute.sting.analyzecovariates; \ No newline at end of file diff --git a/public/java/src/org/broadinstitute/sting/commandline/Advanced.java b/public/java/src/org/broadinstitute/sting/commandline/Advanced.java new file mode 100644 index 000000000..7aeefe261 --- /dev/null +++ b/public/java/src/org/broadinstitute/sting/commandline/Advanced.java @@ -0,0 +1,40 @@ +/* + * Copyright (c) 2011, The Broad Institute + * + * Permission is hereby granted, free of charge, to any person + * obtaining a copy of this software and associated documentation + * files (the "Software"), to deal in the Software without + * restriction, including without limitation the rights to use, + * copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following + * conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES + * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT + * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, + * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR + * OTHER DEALINGS IN THE SOFTWARE. + */ + +package org.broadinstitute.sting.commandline; + +import java.lang.annotation.*; + +/** + * Indicates that a walker argument should is considered an advanced option. + * + * @author Mark DePristo + * @version 0.1 + */ +@Documented +@Inherited +@Retention(RetentionPolicy.RUNTIME) +@Target({ElementType.TYPE,ElementType.FIELD}) +public @interface Advanced { +} diff --git a/public/java/src/org/broadinstitute/sting/commandline/ArgumentSource.java b/public/java/src/org/broadinstitute/sting/commandline/ArgumentSource.java index e0e2ac378..8ec0d650a 100644 --- a/public/java/src/org/broadinstitute/sting/commandline/ArgumentSource.java +++ b/public/java/src/org/broadinstitute/sting/commandline/ArgumentSource.java @@ -151,6 +151,14 @@ public class ArgumentSource { return field.isAnnotationPresent(Hidden.class) || field.isAnnotationPresent(Deprecated.class); } + /** + * Is the given argument considered an advanced option when displaying on the command-line argument system. + * @return True if so. False otherwise. + */ + public boolean isAdvanced() { + return field.isAnnotationPresent(Advanced.class); + } + /** * Is this command-line argument dependent on some primitive argument types? * @return True if this command-line argument depends on other arguments; false otherwise. diff --git a/public/java/src/org/broadinstitute/sting/commandline/ArgumentTypeDescriptor.java b/public/java/src/org/broadinstitute/sting/commandline/ArgumentTypeDescriptor.java index d1d4ff914..ff992d77d 100644 --- a/public/java/src/org/broadinstitute/sting/commandline/ArgumentTypeDescriptor.java +++ b/public/java/src/org/broadinstitute/sting/commandline/ArgumentTypeDescriptor.java @@ -325,7 +325,7 @@ class RodBindingArgumentTypeDescriptor extends ArgumentTypeDescriptor { @Override public Object createTypeDefault(ParsingEngine parsingEngine, ArgumentSource source, Type type) { - Class parameterType = getParameterizedTypeClass(type); + Class parameterType = JVMUtils.getParameterizedTypeClass(type); return RodBinding.makeUnbound((Class)parameterType); } @@ -338,6 +338,8 @@ class RodBindingArgumentTypeDescriptor extends ArgumentTypeDescriptor { public Object parse(ParsingEngine parsingEngine, ArgumentSource source, Type type, ArgumentMatches matches) { ArgumentDefinition defaultDefinition = createDefaultArgumentDefinition(source); String value = getArgumentValue( defaultDefinition, matches ); + Class parameterType = JVMUtils.getParameterizedTypeClass(type); + try { String name = defaultDefinition.fullName; String tribbleType = null; @@ -372,19 +374,19 @@ class RodBindingArgumentTypeDescriptor extends ArgumentTypeDescriptor { FeatureManager.FeatureDescriptor featureDescriptor = manager.getByFiletype(file); if ( featureDescriptor != null ) { tribbleType = featureDescriptor.getName(); - logger.warn("Dynamically determined type of " + file + " to be " + tribbleType); + logger.info("Dynamically determined type of " + file + " to be " + tribbleType); } } + + if ( tribbleType == null ) + throw new UserException.CommandLineException( + String.format("No tribble type was provided on the command line and the type of the file could not be determined dynamically. " + + "Please add an explicit type tag :NAME listing the correct type from among the supported types:%n%s", + manager.userFriendlyListOfAvailableFeatures(parameterType))); } } - if ( tribbleType == null ) // error handling - throw new UserException.CommandLineException( - String.format("Could not parse argument %s with value %s", - defaultDefinition.fullName, value)); - Constructor ctor = (makeRawTypeIfNecessary(type)).getConstructor(Class.class, String.class, String.class, String.class, Tags.class); - Class parameterType = getParameterizedTypeClass(type); RodBinding result = (RodBinding)ctor.newInstance(parameterType, name, value, tribbleType, tags); parsingEngine.addTags(result,tags); parsingEngine.addRodBinding(result); @@ -395,20 +397,10 @@ class RodBindingArgumentTypeDescriptor extends ArgumentTypeDescriptor { value, source.field.getName())); } catch (Exception e) { throw new UserException.CommandLineException( - String.format("Failed to parse value %s for argument %s.", - value, source.field.getName())); + String.format("Failed to parse value %s for argument %s. Message: %s", + value, source.field.getName(), e.getMessage())); } } - - private Class getParameterizedTypeClass(Type t) { - if ( t instanceof ParameterizedType ) { - ParameterizedType parameterizedType = (ParameterizedType)t; - if ( parameterizedType.getActualTypeArguments().length != 1 ) - throw new ReviewedStingException("BUG: more than 1 generic type found on class" + t); - return (Class)parameterizedType.getActualTypeArguments()[0]; - } else - throw new ReviewedStingException("BUG: could not find generic type on class " + t); - } } /** diff --git a/public/java/src/org/broadinstitute/sting/commandline/Output.java b/public/java/src/org/broadinstitute/sting/commandline/Output.java index 22565dbf5..f8aef0355 100644 --- a/public/java/src/org/broadinstitute/sting/commandline/Output.java +++ b/public/java/src/org/broadinstitute/sting/commandline/Output.java @@ -55,7 +55,7 @@ public @interface Output { * --help argument is specified. * @return Doc string associated with this command-line argument. */ - String doc() default "An output file presented to the walker. Will overwrite contents if file exists."; + String doc() default "An output file created by the walker. Will overwrite contents if file exists"; /** * Is this argument required. If true, the command-line argument system will diff --git a/public/java/src/org/broadinstitute/sting/gatk/CommandLineExecutable.java b/public/java/src/org/broadinstitute/sting/gatk/CommandLineExecutable.java index 32132c7ca..32002e093 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/CommandLineExecutable.java +++ b/public/java/src/org/broadinstitute/sting/gatk/CommandLineExecutable.java @@ -96,24 +96,23 @@ public abstract class CommandLineExecutable extends CommandLineProgram { loadArgumentsIntoObject(walker); argumentSources.add(walker); - Collection
+ * It allows the user to clip bases in reads with poor quality scores, that match particular + * sequences, or that were generated by particular machine cycles. + * + *
+ * Any number of BAM files. + *
+ * + *+ * A new BAM file containing all of the reads from the input BAMs with the user-specified clipping + * operation applied to each read. + *
+ *+ *
+ * Number of examined reads 13 + * Number of clipped reads 13 + * Percent of clipped reads 100.00 + * Number of examined bases 988 + * Number of clipped bases 126 + * Percent of clipped bases 12.75 + * Number of quality-score clipped bases 126 + * Number of range clipped bases 0 + * Number of sequence clipped bases 0 + *+ * + * + *
+ *
+ * 314KGAAXX090507:1:19:1420:1123#0 16 chrM 3116 29 76M * * * + * TAGGACCCGGGCCCCCCTCCCCAATCCTCCAACGCATATAGCGGCCGCGCCTTCCCCCGTAAATGATATCATCTCA + * #################4?6/?2135;;;'1/=/<'B9;12;68?A79@,@==@9?=AAA3;A@B;A?B54;?ABA + *+ * + * If we are clipping reads with -QT 10 and -CR WRITE_NS, we get: + * + *
+ * 314KGAAXX090507:1:19:1420:1123#0 16 chrM 3116 29 76M * * * + * NNNNNNNNNNNNNNNNNTCCCCAATCCTCCAACGCATATAGCGGCCGCGCCTTCCCCCGTAAATGATATCATCTCA + * #################4?6/?2135;;;'1/=/<'B9;12;68?A79@,@==@9?=AAA3;A@B;A?B54;?ABA + *+ * + * Whereas with -CR WRITE_Q0S: + *
+ * 314KGAAXX090507:1:19:1420:1123#0 16 chrM 3116 29 76M * * * + * TAGGACCCGGGCCCCCCTCCCCAATCCTCCAACGCATATAGCGGCCGCGCCTTCCCCCGTAAATGATATCATCTCA + * !!!!!!!!!!!!!!!!!4?6/?2135;;;'1/=/<'B9;12;68?A79@,@==@9?=AAA3;A@B;A?B54;?ABA + *+ * + * Or -CR SOFTCLIP_BASES: + *
+ * 314KGAAXX090507:1:19:1420:1123#0 16 chrM 3133 29 17S59M * * * + * TAGGACCCGGGCCCCCCTCCCCAATCCTCCAACGCATATAGCGGCCGCGCCTTCCCCCGTAAATGATATCATCTCA + * #################4?6/?2135;;;'1/=/<'B9;12;68?A79@,@==@9?=AAA3;A@B;A?B54;?ABA + *+ * + * + *
+ * -T ClipReads -I my.bam -I your.bam -o my_and_your.clipped.bam -R Homo_sapiens_assembly18.fasta \ + * -XF seqsToClip.fasta -X CCCCC -CT "1-5,11-15" -QT 10 + *+ + * @author Mark DePristo + * @since 2010 */ @Requires({DataSource.READS}) public class ClipReadsWalker extends ReadWalker
+ * PrintReads can dynamically merge the contents of multiple input BAM files, resulting + * in merged output sorted in coordinate order. Can also optionally filter reads based on the + * --read_filter command line argument. + * + *
+ * One or more bam files. + *
+ * + *+ * A single processed bam file. + *
+ * + *+ * java -Xmx2g -jar GenomeAnalysisTK.jar \ + * -R ref.fasta \ + * -T PrintReads \ + * -o output.bam \ + * -I input1.bam \ + * -I input2.bam \ + * --read_filter MappingQualityZero + *+ * */ @BAQMode(QualityMode = BAQ.QualityMode.ADD_TAG, ApplicationTime = BAQ.ApplicationTime.ON_OUTPUT) @Requires({DataSource.READS, DataSource.REFERENCE}) public class PrintReadsWalker extends ReadWalker
+ * VariantAnnotator is a GATK tool for annotating variant calls based on their context. + * The tool is modular; new annotations can be written easily without modifying VariantAnnotator itself. + * + *
+ * A variant set to annotate and optionally one or more BAM files. + *
+ * + *+ * An annotated VCF. + *
+ * + *+ * java -Xmx2g -jar GenomeAnalysisTK.jar \ + * -R ref.fasta \ + * -T VariantAnnotator \ + * -I input.bam \ + * -o output.vcf \ + * -A DepthOfCoverage + * --variant input.vcf \ + * --dbsnp dbsnp.vcf + *+ * */ @Requires(value={}) @Allows(value={DataSource.READS, DataSource.REFERENCE}) @@ -69,8 +96,6 @@ public class VariantAnnotator extends RodWalker
This walker is intended to be run after Beagle has successfully executed. The full calling sequence for using Beagle along with the GATK is:
+ * + *1. Run ProduceBeagleInputWalker.
+ *2. Run Beagle
+ *3. Uncompress output files
+ *4. Run BeagleOutputToVCFWalker.
+ * + * + * Note that this walker requires all input files produced by Beagle. + * + * + *+ * java -Xmx4000m -jar dist/GenomeAnalysisTK.jar \ + * -R reffile.fasta -T BeagleOutputToVCF \ + * -V input_vcf.vcf \ + * -beagleR2:BEAGLE /myrun.beagle_output.r2 \ + * -beaglePhased:BEAGLE /myrun.beagle_output.phased \ + * -beagleProbs:BEAGLE /myrun.beagle_output.gprobs \ + * -o output_vcf.vcf + *+ +
Note that Beagle produces some of these files compressed as .gz, so gunzip must be run on them before walker is run in order to decompress them
+ */ public class BeagleOutputToVCFWalker extends RodWalker+ * + *
+ * A VCF with variants to convert to Beagle format + *
+ * + *+ * A single text file which can be fed to Beagle + *
+ *+ * Optional: A file with a list of markers + *
+ * + *+ * java -Xmx2g -jar dist/GenomeAnalysisTK.jar -L 20 \ + * -R reffile.fasta -T ProduceBeagleInput \ + * -V path_to_input_vcf/inputvcf.vcf -o path_to_beagle_output/beagle_output + *+ * */ + public class ProduceBeagleInputWalker extends RodWalker
+ * A very common question about a NGS set of reads is what areas of the genome are considered callable. The system + * considers the coverage at each locus and emits either a per base state or a summary interval BED file that + * partitions the genomic intervals into the following callable states: + *
+ * A BAM file containing exactly one sample. + *
+ * + *+ *
+ * -T CallableLociWalker \ + * -I my.bam \ + * -summary my.summary \ + * -o my.bed + *+ * + * would produce a BED file (my.bed) that looks like: + * + *
+ * 20 10000000 10000864 CALLABLE + * 20 10000865 10000985 POOR_MAPPING_QUALITY + * 20 10000986 10001138 CALLABLE + * 20 10001139 10001254 POOR_MAPPING_QUALITY + * 20 10001255 10012255 CALLABLE + * 20 10012256 10012259 POOR_MAPPING_QUALITY + * 20 10012260 10012263 CALLABLE + * 20 10012264 10012328 POOR_MAPPING_QUALITY + * 20 10012329 10012550 CALLABLE + * 20 10012551 10012551 LOW_COVERAGE + * 20 10012552 10012554 CALLABLE + * 20 10012555 10012557 LOW_COVERAGE + * 20 10012558 10012558 CALLABLE + * et cetera... + *+ * as well as a summary table that looks like: + * + *
+ * state nBases + * REF_N 0 + * CALLABLE 996046 + * NO_COVERAGE 121 + * LOW_COVERAGE 928 + * EXCESSIVE_COVERAGE 0 + * POOR_MAPPING_QUALITY 2906 + *+ * + * @author Mark DePristo + * @since May 7, 2010 */ @By(DataSource.REFERENCE) public class CallableLociWalker extends LocusWalker
+ * One or more bam files (with proper headers) to be analyzed for coverage statistics + * (Optional) A REFSEQ Rod to aggregate coverage to the gene level + *
+ * + *+ * Tables pertaining to different coverage summaries. Suffix on the table files declares the contents: + * - no suffix: per locus coverage + * - _summary: total, mean, median, quartiles, and threshold proportions, aggregated over all bases + * - _statistics: coverage histograms (# locus with X coverage), aggregated over all bases + * - _interval_summary: total, mean, median, quartiles, and threshold proportions, aggregated per interval + * - _interval_statistics: 2x2 table of # of intervals covered to >= X depth in >=Y samples + * - _gene_summary: total, mean, median, quartiles, and threshold proportions, aggregated per gene + * - _gene_statistics: 2x2 table of # of genes covered to >= X depth in >= Y samples + * - _cumulative_coverage_counts: coverage histograms (# locus with >= X coverage), aggregated over all bases + * - _cumulative_coverage_proportions: proprotions of loci with >= X coverage, aggregated over all bases + *
+ * + *+ * java -Xmx2g -jar GenomeAnalysisTK.jar \ + * -R ref.fasta \ + * -T VariantEval \ + * -o file_name_base \ + * -I input_bams.list + * [-geneList refSeq.sorted.txt] \ + * [-pt readgroup] \ + * [-ct 4 -ct 6 -ct 10] \ + * [-L my_capture_genes.interval_list] + ** - * @Author chartl - * @Date Feb 22, 2010 */ // todo -- cache the map from sample names to means in the print functions, rather than regenerating each time // todo -- support for granular histograms for total depth; maybe n*[start,stop], bins*sqrt(n) diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/coverage/GCContentByIntervalWalker.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/coverage/GCContentByIntervalWalker.java index a4944e939..5c2a967b9 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/coverage/GCContentByIntervalWalker.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/coverage/GCContentByIntervalWalker.java @@ -38,12 +38,32 @@ import java.util.List; /** * Walks along reference and calculates the GC content for each interval. + * + * + *
+ * One or more BAM files. + *
+ * + *+ * GC content calculations per interval. + *
+ * + *+ * java -Xmx2g -jar GenomeAnalysisTK.jar \ + * -R ref.fasta \ + * -T GCContentByInterval \ + * -o output.txt \ + * -I input.bam \ + * -L input.intervals + *+ * */ @Allows(value = {DataSource.REFERENCE}) @Requires(value = {DataSource.REFERENCE}) - @By(DataSource.REFERENCE) - public class GCContentByIntervalWalker extends LocusWalker
+ * Given variant ROD tracks, it replaces the reference bases at variation sites with the bases supplied by the ROD(s). + * Additionally, allows for a "snpmask" ROD to set overlapping bases to 'N'. + * + *
+ * The reference, requested intervals, and any number of variant rod files. + *
+ * + *+ * A fasta file representing the requested intervals. + *
+ * + *+ * java -Xmx2g -jar GenomeAnalysisTK.jar \ + * -R ref.fasta \ + * -T FastaAlternateReferenceMaker \ + * -o output.fasta \ + * -L input.intervals \ + * --variant input.vcf \ + * [--snpmask mask.vcf] + *+ * */ @WalkerName("FastaAlternateReferenceMaker") @Reference(window=@Window(start=-1,stop=50)) @Requires(value={DataSource.REFERENCE}) public class FastaAlternateReferenceWalker extends FastaReferenceWalker { + /** + * Variants from these input files are used by this tool to construct an alternate reference. + */ @Input(fullName = "variant", shortName = "V", doc="variants to model", required=false) - public List
+ * The output format can be partially controlled using the provided command-line arguments. + * + *
+ * The reference and requested intervals. + *
+ * + *+ * A fasta file representing the requested intervals. + *
+ * + *+ * java -Xmx2g -jar GenomeAnalysisTK.jar \ + * -R ref.fasta \ + * -T FastaReference \ + * -o output.fasta \ + * -L input.intervals + *+ * */ @WalkerName("FastaReferenceMaker") public class FastaReferenceWalker extends RefWalker
+ * VariantFiltration is a GATK tool for hard-filtering variant calls based on certain criteria. + * Records are hard-filtered by changing the value in the FILTER field to something other than PASS. + * + *
+ * A variant set to filter. + *
+ * + *+ * A filtered VCF. + *
+ * + *+ * java -Xmx2g -jar GenomeAnalysisTK.jar \ + * -R ref.fasta \ + * -T VariantFiltration \ + * -o output.vcf \ + * --variant input.vcf \ + * --filterExpression "AB < 0.2 || MQ0 > 50" \ + * --filterName "Nov09filters" \ + * --mask mask.vcf \ + * --maskName InDel + *+ * */ @Reference(window=@Window(start=-50,stop=50)) public class VariantFiltrationWalker extends RodWalker
+ * The GATK Unified Genotyper is a multiple-sample, technology-aware SNP and indel caller. It uses a Bayesian genotype + * likelihood model to estimate simultaneously the most likely genotypes and allele frequency in a population of N samples, + * emitting an accurate posterior probability of there being a segregating variant allele at each locus as well as for the + * genotype of each sample. The system can either emit just the variant sites or complete genotypes (which includes + * homozygous reference calls) satisfying some phred-scaled confidence value. The genotyper can make accurate calls on + * both single sample data and multi-sample data. + * + *
+ * The read data from which to make variant calls. + *
+ * + *+ * A raw, unfiltered, highly specific callset in VCF format. + *
+ * + *+ * java -jar GenomeAnalysisTK.jar \ + * -R resources/Homo_sapiens_assembly18.fasta \ + * -T UnifiedGenotyper \ + * -I sample1.bam [-I sample2.bam ...] \ + * --dbsnp dbSNP.vcf \ + * -o snps.raw.vcf \ + * -stand_call_conf [50.0] \ + * -stand_emit_conf 10.0 \ + * -dcov [50] \ + * [-L targets.interval_list] + *+ * + *
+ * The above command will call all of the samples in your provided BAM files [-I arguments] together and produce a VCF file + * with sites and genotypes for all samples. The easiest way to get the dbSNP file is from the GATK resource bundle. Several + * arguments have parameters that should be chosen based on the average coverage per sample in your data. See the detailed + * argument descriptions below. + *
+ * + *+ * java -jar /path/to/GenomeAnalysisTK.jar \ + * -l INFO \ + * -R resources/Homo_sapiens_assembly18.fasta \ + * -T UnifiedGenotyper \ + * -I /DCC/ftp/pilot_data/data/NA12878/alignment/NA12878.SLX.maq.SRP000031.2009_08.bam \ + * -o my.vcf \ + * --output_mode EMIT_ALL_SITES + *+ * + *
+ * The local realignment tool is designed to consume one or more BAM files and to locally realign reads such that the number of mismatching bases + * is minimized across all the reads. In general, a large percent of regions requiring local realignment are due to the presence of an insertion + * or deletion (indels) in the individual's genome with respect to the reference genome. Such alignment artifacts result in many bases mismatching + * the reference near the misalignment, which are easily mistaken as SNPs. Moreover, since read mapping algorithms operate on each read independently, + * it is impossible to place reads on the reference genome such at mismatches are minimized across all reads. Consequently, even when some reads are + * correctly mapped with indels, reads covering the indel near just the start or end of the read are often incorrectly mapped with respect the true indel, + * also requiring realignment. Local realignment serves to transform regions with misalignments due to indels into clean reads containing a consensus + * indel suitable for standard variant discovery approaches. Unlike most mappers, this walker uses the full alignment context to determine whether an + * appropriate alternate reference (i.e. indel) exists. Following local realignment, the GATK tool Unified Genotyper can be used to sensitively and + * specifically identify indels. + *
+ *
+ * An important note: the input bam(s), reference, and known indel file(s) should be the same ones used for the RealignerTargetCreator step. + *
+ * Another important note: because reads produced from the 454 technology inherently contain false indels, the realigner will not currently work with them + * (or with reads from similar technologies). + * + *
+ * One or more aligned BAM files and optionally one or more lists of known indels. + *
+ * + *+ * A realigned version of your input BAM file(s). + *
+ * + *+ * java -Xmx4g -jar GenomeAnalysisTK.jar \ + * -I input.bam \ + * -R ref.fasta \ + * -T IndelRealigner \ + * -targetIntervals intervalListFromRTC.intervals \ + * -o realignedBam.bam \ + * [--known /path/to/indels.vcf] \ + * [-compress 0] (this argument recommended to speed up the process *if* this is only a temporary file; otherwise, use the default value) + *+ * + * @author ebanks */ -//Reference(window=@Window(start=-30,stop=30)) @BAQMode(QualityMode = BAQ.QualityMode.ADD_TAG, ApplicationTime = BAQ.ApplicationTime.ON_OUTPUT) public class IndelRealigner extends ReadWalker
+ * LeftAlignIndels is a tool that takes a bam file and left-aligns any indels inside it. The same indel can often be + * placed at multiple positions and still represent the same haplotype. While a standard convention is to place an + * indel at the left-most position this doesn't always happen, so this tool can be used to left-align them. + * + *
+ * A bam file to left-align. + *
+ * + *+ * A left-aligned bam. + *
+ * + *+ * java -Xmx3g -jar GenomeAnalysisTK.jar \ + * -R ref.fasta \ + * -T LeftAlignIndels \ + * -I input.bam \ + * -o output.vcf + *+ * */ public class LeftAlignIndels extends ReadWalker
+ * The local realignment tool is designed to consume one or more BAM files and to locally realign reads such that the number of mismatching bases + * is minimized across all the reads. In general, a large percent of regions requiring local realignment are due to the presence of an insertion + * or deletion (indels) in the individual's genome with respect to the reference genome. Such alignment artifacts result in many bases mismatching + * the reference near the misalignment, which are easily mistaken as SNPs. Moreover, since read mapping algorithms operate on each read independently, + * it is impossible to place reads on the reference genome such at mismatches are minimized across all reads. Consequently, even when some reads are + * correctly mapped with indels, reads covering the indel near just the start or end of the read are often incorrectly mapped with respect the true indel, + * also requiring realignment. Local realignment serves to transform regions with misalignments due to indels into clean reads containing a consensus + * indel suitable for standard variant discovery approaches. Unlike most mappers, this walker uses the full alignment context to determine whether an + * appropriate alternate reference (i.e. indel) exists. Following local realignment, the GATK tool Unified Genotyper can be used to sensitively and + * specifically identify indels. + *
+ *
+ * An important note: the input bam(s), reference, and known indel file(s) should be the same ones to be used for the IndelRealigner step. + *
+ * Another important note: because reads produced from the 454 technology inherently contain false indels, the realigner will not currently work with them + * (or with reads from similar technologies). This tool also ignores MQ0 reads and reads with consecutive indel operators in the CIGAR string. + * + *
+ * One or more aligned BAM files and optionally one or more lists of known indels. + *
+ * + *+ * A list of target intervals to pass to the Indel Realigner. + *
+ * + *+ * java -Xmx2g -jar GenomeAnalysisTK.jar \ + * -I input.bam \ + * -R ref.fasta \ + * -T RealignerTargetCreator \ + * -o forIndelRealigner.intervals \ + * [--known /path/to/indels.vcf] + *+ * + * @author ebanks */ -@ReadFilters({Platform454Filter.class, MappingQualityZeroReadFilter.class, BadCigarFilter.class}) +@ReadFilters({Platform454Filter.class, MappingQualityZeroFilter.class, BadCigarFilter.class}) @Reference(window=@Window(start=-1,stop=50)) @Allows(value={DataSource.READS, DataSource.REFERENCE}) @By(DataSource.REFERENCE) @BAQMode(ApplicationTime = BAQ.ApplicationTime.FORBIDDEN) public class RealignerTargetCreator extends RodWalker
+ * Performs physical phasing of SNP calls, based on sequencing reads. + *
+ * + *+ * VCF file of SNP calls, BAM file of sequence reads. + *
+ * + *+ * Phased VCF file. + *
+ * + *+ * java + * -jar GenomeAnalysisTK.jar + * -T ReadBackedPhasing + * -R reference.fasta + * -I reads.bam + * --variant:vcf SNPs.vcf + * -BTI variant + * -BTIMR INTERSECTION + * -o phased_SNPs.vcf + * --phaseQualityThresh 20.0 + *+ * + * @author Menachem Fromer + * @since July 2010 */ @Allows(value = {DataSource.READS, DataSource.REFERENCE}) @Requires(value = {DataSource.READS, DataSource.REFERENCE}) @By(DataSource.READS) -@ReadFilters({MappingQualityZeroReadFilter.class}) // Filter out all reads with zero mapping quality +@ReadFilters({MappingQualityZeroFilter.class}) public class ReadBackedPhasingWalker extends RodWalker
* Simplest example of a locus walker. + * + * + *
+ * One or more BAM files. + *
+ * + *+ * Number of loci traversed. + *
+ * + *+ * java -Xmx2g -jar GenomeAnalysisTK.jar \ + * -R ref.fasta \ + * -T CountLoci \ + * -o output.txt \ + * -I input.bam \ + * [-L input.intervals] + *+ * */ public class CountLociWalker extends LocusWalker
+ * One or more bam files. + *
+ * + *+ * Number of pairs seen. + *
+ * + *+ * java -Xmx2g -jar GenomeAnalysisTK.jar \ + * -R ref.fasta \ + * -T CountPairs \ + * -o output.txt \ + * -I input.bam + *+ * * @author mhanna */ public class CountPairsWalker extends ReadPairWalker
+ * One or more rod files. + *
+ * + *+ * Number of rods seen. + *
+ * + *+ * java -Xmx2g -jar GenomeAnalysisTK.jar \ + * -R ref.fasta \ + * -T CountRODsByRef \ + * -o output.txt \ + * --rod input.vcf + *+ * */ -public class CountRodByRefWalker extends RefWalker
+ * One or more rod files. + *
+ * + *+ * Number of rods seen. + *
+ * + *+ * java -Xmx2g -jar GenomeAnalysisTK.jar \ + * -R ref.fasta \ + * -T CountRODs \ + * -o output.txt \ + * --rod input.vcf + *+ * */ -public class CountRodWalker extends RodWalker
* Can also count the number of reads matching a given criterion using read filters (see the * --read-filter command line argument). Simplest example of a read-backed analysis. + * + * + *
+ * One or more BAM files. + *
+ * + *+ * Number of reads seen. + *
+ * + *+ * java -Xmx2g -jar GenomeAnalysisTK.jar \ + * -R ref.fasta \ + * -T CountReads \ + * -o output.txt \ + * -I input.bam \ + * [-L input.intervals] + *+ * */ @Requires({DataSource.READS, DataSource.REFERENCE}) public class CountReadsWalker extends ReadWalker
Body test
+ */ +public class DocumentationTest extends RodWalker+ * This walker is designed to work as the first pass in a two-pass processing step. It does a by-locus traversal operating + * only at sites that are not in dbSNP. We assume that all reference mismatches we see are therefore errors and indicative + * of poor base quality. This walker generates tables based on various user-specified covariates (such as read group, + * reported quality score, cycle, and dinucleotide). Since there is a large amount of data one can then calculate an empirical + * probability of error given the particular covariates seen at this site, where p(error) = num mismatches / num observations. + * The output file is a CSV list of (the several covariate values, num observations, num mismatches, empirical quality score). + *
+ * Note: ReadGroupCovariate and QualityScoreCovariate are required covariates and will be added for the user regardless of whether or not they were specified. * - * Note: ReadGroupCovariate and QualityScoreCovariate are required covariates and will be added for the user regardless of whether or not they were specified - * Note: This walker is designed to be used in conjunction with TableRecalibrationWalker. + *
+ * See the GATK wiki for a tutorial and example recalibration accuracy plots. + * http://www.broadinstitute.org/gsa/wiki/index.php/Base_quality_score_recalibration + * + *
+ * The input read data whose base quality scores need to be assessed. + *
+ * A database of known polymorphic sites to skip over. + *
+ * + *+ * A recalibration table file in CSV format that is used by the TableRecalibration walker. + *
+ * + *+ * java -Xmx4g -jar GenomeAnalysisTK.jar \ + * -R resources/Homo_sapiens_assembly18.fasta \ + * -knownSites bundle/hg18/dbsnp_132.hg18.vcf \ + * -knownSites another/optional/setOfSitesToMask.vcf \ + * -I my_reads.bam \ + * -T CountCovariates \ + * -cov ReadGroupCovariate \ + * -cov QualityScoreCovariate \ + * -cov CycleCovariate \ + * -cov DinucCovariate \ + * -recalFile my_reads.recal_data.csv + ** - * @author rpoplin - * @since Nov 3, 2009 */ @BAQMode(ApplicationTime = BAQ.ApplicationTime.FORBIDDEN) @By( DataSource.READS ) // Only look at covered loci, not every loci of the reference file -@ReadFilters( {MappingQualityZeroReadFilter.class, MappingQualityUnavailableReadFilter.class} ) // Filter out all reads with zero or unavailable mapping quality +@ReadFilters( {MappingQualityZeroFilter.class, MappingQualityUnavailableFilter.class} ) // Filter out all reads with zero or unavailable mapping quality @Requires( {DataSource.READS, DataSource.REFERENCE, DataSource.REFERENCE_BASES} ) // This walker requires both -I input.bam and -R reference.fasta @PartitionBy(PartitionType.LOCUS) public class CountCovariatesWalker extends LocusWalker
+ * This walker is designed to work as the second pass in a two-pass processing step, doing a by-read traversal. For each + * base in each read this walker calculates various user-specified covariates (such as read group, reported quality score, + * cycle, and dinuc). Using these values as a key in a large hashmap the walker calculates an empirical base quality score + * and overwrites the quality score currently in the read. This walker then outputs a new bam file with these updated (recalibrated) reads. * - * For each base in each read this walker calculates various user-specified covariates (such as read group, reported quality score, cycle, and dinuc) - * Using these values as a key in a large hashmap the walker calculates an empirical base quality score and overwrites the quality score currently in the read. - * This walker then outputs a new bam file with these updated (recalibrated) reads. + *
+ * See the GATK wiki for a tutorial and example recalibration accuracy plots. + * http://www.broadinstitute.org/gsa/wiki/index.php/Base_quality_score_recalibration * - * Note: This walker expects as input the recalibration table file generated previously by CovariateCounterWalker. - * Note: This walker is designed to be used in conjunction with CovariateCounterWalker. + *
+ * The input read data whose base quality scores need to be recalibrated. + *
+ * The recalibration table file in CSV format that was generated by the CountCovariates walker. + *
+ * + *+ * A bam file in which the quality scores in each read have been recalibrated. + *
+ * + *+ * java -Xmx4g -jar GenomeAnalysisTK.jar \ + * -R resources/Homo_sapiens_assembly18.fasta \ + * -I my_reads.bam \ + * -T TableRecalibration \ + * -o my_reads.recal.bam \ + * -recalFile my_reads.recal_data.csv + ** - * @author rpoplin - * @since Nov 3, 2009 */ @BAQMode(QualityMode = BAQ.QualityMode.ADD_TAG, ApplicationTime = BAQ.ApplicationTime.ON_OUTPUT) @@ -79,24 +100,54 @@ public class TableRecalibrationWalker extends ReadWalker