Bug Reports and Issue Tracking
mucor is software to aggregate variant information sourced from multiple VCF files (and some others, see below) into a variety of summary files with varying levels of detail and statistics. The outputs range from high-level summaries to full, detailed reports in text and Microsoft Excel formats. The intended audience is both biologists and bioinformaticians.
The following modules are required:
- numpy (http://www.numpy.org/)
- pandas (http://pandas.pydata.org/)
- HTSeq (http://www-huber.embl.de/HTSeq/)
The following modules are optional:
- pytabix (https://github.com/slowkow/pytabix)
- XlsxWriter (https://github.com/jmcnamara/XlsxWriter)
- xlwt (https://pypi.python.org/pypi/xlwt)
- bgzip and Tabix (optional, see Databases, below; https://github.com/samtools/htslib)
mucor requires a reference annotation file (in GTF or GFF3 format) for feature definition. Contig (chromosome) names must match between the input VCF and reference annotation file; for example, for Homo sapiens, UCSC uses the form chr17 while Ensembl uses the form 17.
Protip: If your human data are aligned and variants called with chr style contig names but you wish to use Ensembl genes, you can try the GENCODE annotation which contains most Ensembl genes but using chr style contig numbering.
mucor can optionally check variants against any number of supplied databases and report presence (including identifier) or absence in the database. Common choices would be one or more releases of dbSNP (e.g., dbSNP 137 and latest dbSNP; dbSNP clinvar and dbSNP all), 1000 Genomes, NHLBI 6500 Exomes, or COSMIC. Users may also supply their own custom databases.
Databses must meet the following requirements:
- Conform to VCF standard format
- be compressed with bgzip (not gzip; see below)
- be indexed with tabix (see below).
bgzip and tabix provide massive speedups looking up variants in these external files. The database lookup feature is not available if files are not correctly bgzipped or tabix indexed, or if the Pytabix module is not installed. bgzip and tabix are available as part of htslib (https://github.com/samtools/htslib)
Beginning with an example VCF dbSNP.vcf, execute the following:
$ bgzip dbSNP.vcf
$ tabix -p vcf dbSNP.vcf.gz
Note that bgzipped files can be treated as regular .gz files by other tools.
These data are the data that you wish to aggregate and summarize. In theory, mucor will work with any well-formed generic VCF, but it has been tested with and contains specific code to handle VCF files generated by the following tools:
- snpEff
- Ion Torrent PGM (default machine output)
- Illumina Miseq (default machine output)
- GATK
- GATK SomaticIndelDetector
- GATK HaplotypeCaller
- MuTect
- VarScan
- FreeBayes
- Samtools
In addition, mucor can read the more detailed .out files produced by MuTect.
clone from https://github.com/blachlylab
Running mucor to aggregate and summarize variants is a two step process:
-
Project setup / configuration
mucor_config.pycreates a JSON config file
-
Variant Aggregation!
mucor.py
The configuration step is completed using the provided mucor_config.py utility. It accepts the following command line arguments and creates a JSON file that will be passed to the main mucor script.
-ex, --example Print a valid, example JSON config file and exit.
Function that will write a template of the JSON config. It can be edited manually and supplied to mucor.
-g GFF, --gff GFF
Reference annotation GFF/GTF for feature binning. Required
-f FEATURETYPE, --featuretype FEATURETYPE
Feature type into which to bin. Gencode GTF example: gene_name, gene_id, transcript_name, transcript_id, etc. Required
-db DATABASES, --databases DATABASES
Colon delimited name and path to variant database in bgzipped VCF format. Can be declared >= 0 times. Ex: -db name1:/full/user/path/name1.vcf.gz. Optional
-s SAMPLES, --samples SAMPLES
Text file containing sample names. One sample per line. mucor_config.py attempts to guess which files belong with which sample IDs using globbing (wildcard filename matching). This means that the auto-configuration may be incorrect if any sample names are contained within another sample name. Ex: U-23 and U-238. U-23 would erroneously identify U-238 files, requiring manual modification of the JSON file. Sample IDs U-023 and U-238 would not exhibit this problem. Required
-d PROJECT_DIRECTORY, --project_directory PROJECT_DIRECTORY
Working/project directory, in which to find variant call files to aggregate. Variant calls can be in the provided directory, or any of its subdirectories. Default: current working directory
-vcff VCF_FILTERS, --vcf_filters VCF_FILTERS
Comma separated list of VCF filters to allow. Default: PASS
-a ARCHIVE_DIRECTORY, --archive_directory ARCHIVE_DIRECTORY
Specify directory in which to read/write archived annotations. This step will significantly speed up future runs that use the same annotation and feature type, even if the sample data changes. Undefined will prevent using the annotation archive features. Optional
-r REGIONS, --regions REGIONS
Comma separated list of bed regions and/or bed files by which to limit output. Bed regions can be specific positions, or entire chromosomes. Ex: chr1:10230-10240,chr2,my_regions.bed. Optional
-u, --union
Join all items with same ID for feature_type (specified by -f) into a single, continuous bin. For example, if you want intronic variants counted with a gene, use this option. WARNING, this will lead to spurious results for features that are duplicated on the same contig. When feature names are identical, the bin will range from the beginning of the first instance to the end of the last, even if they are several megabases apart. Refer to the documentation for a resolution using 'detect_union_bin_errors.py.' Optional.
-jco JSON_CONFIG_OUTPUT, --json_config_output JSON_CONFIG_OUTPUT
Name of JSON configuration output file. This is the configuration file fed into mucor. Required
-outd OUTPUT_DIRECTORY, --output_directory OUTPUT_DIRECTORY
Name of directory in which to write mucor output. Required
-outt OUTPUT_TYPE, --output_type OUTPUT_TYPE
Comma separated list of desired output types. Options include: counts, txt, longtxt, xls, longxls, bed, featXsamp, featmutXsamp, all. Default: counts,txt
(See the detailed description of Output File Formats, below)
Example:
python ./mucor_config.py -g ~/references/human/gencode/gencode.v19.annotation.gtf -f gene_name -s samples.txt -a ./fast -u -jco mucor_config.json -outd ./mucor_output -outt all -db 1000G:~/references/human/1000_genomes/chrm_1000Genomes.20130502.genotypes.vcf.gz -db dbsnp:~/references/human/dbsnp/common_all.hg19.sorted.leftalign.vcf.gz
mucor_config.py produces a JSON-formatted configuration file embodying the selected options for the subsequent mucor run. The configuration could be edited manually (e.g., to tweak a database name or path) or programmatically (e.g., if mucor_config.py's guesses about sample IDs were incomplete) at this stage.
Alternatively, mucor_config.py -ex produces a syntactically valid example JSON file for editing.
Producing a configuration file is intended to facilitate the following:
- Consistency between runs
- Documentation of settings
- Easier use with dozens to thousands of input files
Execution: mucor.py <config.json>
mucor is executed by launching the main mucor.py script, passing as a single parameter the name of the previously generated/edited JSON file.
Output files of the specified type are placed in the output directory specified during configuration.
Output types are specified at the time of configuration. The user may select any number and combination of output types from the list below.
all Execute all output types
counts
Print counts of mutations per feature. Output: counts.txt
txt
Print all information about each variant, with one-per-row, irrespective of how many samples in which it appears. Useful for variant-centric projects. Identical to xls in layout. Output: variant_details.txt
longtxt
Similar to txt above, but writes each instance of a variant to a new row. Each variant is written once per source file, instead of combining recurrent variations into 1 unique row. Identical to longxls in layout. Output: long_variant_details.txt
xls
Print all information about each variant, with one-per-row, irrespective of how many samples in which it appears. Useful for variant-centric projects. Identical to txt in layout. Output: variant_details.xls/xlsx
longxls
Similar to xls above, but writes each instance of a variant to a new row. Each variant is written once per source file, instead of combining recurrent variations into 1 unique row. Identical to longtxt in layout. Output: long_variant_details.xls/xlsx
NB: The XLS format has a hard limit of 2^16 rows; in long record format, a moderate sized study could exceed this (2,000 total variants/sample * 32 samples = 65,536 rows). mucor can use Python's xlwt module to write .xls format, but it is preferrable to have XlsxWriter or openPyxl installed for .xlsx support.
bed
Print bed file of the variant locations Output: variant_locations.bed
vcf
Print vcf file of the variant locations, features, depths, and variant frequencies. Output: variant_locations.vcf
featXsamp
Print table of mutation counts per feature per sample. Samples are in columns, while features are in rows. The count of unique mutations per sample per feature are the table values. This output is useful for examining patterns in variation across samples, for example, to look at combinatoric mutation status for selected recurrently mutated genes. This output could be used directly to make a heatmap. Output: feature_by_sample.txt
mutXsamp
Print table of mutations per sample. Unlike featXsamp, this differentiates among different variants within the same features. For example, in acute leukemia, the functional effect of mutations in DNMT3A depends on whether it is an R882 mutation or non-R882 mutation. As before, samples are in columns, with features in rows. However, rows 2-4 contain information about chromosome, position, ref, and alt. The table values are boolean: 1 for present mutation, 0 for missing mutation. This output could be used directly or with appropriate filtering to make an Oncoprint. Output: feature_and_mutation_by_sample.txt
mutXsampVAF
Print table of mutations per sample. Identical to mutXsamp except boolean values are replaced by the respective variant VAF. Output: feature_and_mutation_by_sample.txt
DepthGauge, a companion utility to mucor, measures coverage at regions of interest to increase confidence in variant calls. The purpose of this tool is not only to verify sufficient depth at mutations, but more importantly, to verify sufficient depth in places that are not called mutant. The lack of a mutation call may indicate no mutation, or simply that the region of interest had insufficient coverage for analysis. Like mucor, DepthGauge is dependent only on the JSON formatted configuration file which may be entirely hand created, entirely automatically generated by mucor_config.py, or automatically generated and hand-tuned. DepthGauge has three additional options, the latter of which can override the regions, if any, specified in the JSON configuration file.
depth_gauge.py config_file.json [-p] [-c] [-r REGIONS]
-p, --point Instead of reporting an average depth for every location within a window (default), take the middle coordinate within the range and calculate the depth at that point as a surrogate for the entire region
-c, --count Instead of reporting an average depth for every location within a window (default), or taking a point estimate from the middle (-p), instead count the total number of reads mapped within the region.
-r REGIONS This option specifies a list of regions to query for depth. If present, it overrides the region(s) specified in the JSON configuration file. As an alternative, the JSON configuration file could be edited before running DepthGauge.
Output: Depth_of_Coverage.xlsx
The --union feature will behave inappropriately when genomic feature names are duplicated on the same contig. For example, if gene "ABC" is duplicated on the beginning and the end of chromosome 1, the feature bin for gene "ABC" will cover the whole contig (from the beginning of the first copy of "ABC", to the end of the last copy). Users may select another feature type, such as gene_id, which is unique to every copy of a gene. Otherwise, users may run the included python script [detect_union_bin_errors.py] to detect potential bin errors. It accepts a feature_type, a GTF/GFF annotation, and an output directory. The output is a text file list of feature names likely to cause large bin errors. Place this text document into the working directory where mucor will be executed. mucor will automatically search for the text file by name ['union_incompatible_genes.txt'] in the current directory and print a message indicating when the file is found.
python ./detect_union_bin_errors.py -o ~/projects/mucor -g ~/references/human/gencode/gencode.v19.annotation.gtf -f gene_name
There is an issue when a variant file presents a contig that the pickled (archived) annotation does not have. This will throw a warning that shows how many contigs were unknown, and how many mutations were encountered on these contigs. The solution is to disable the archive feature by omitting the -a or --archive_directory option. This will permit the unknown contig in output, but all mutations on the unknown contig will be shown as having no feature.
*** WARNING: 18 Contigs and 39 mutations are in areas unknown to the genomic array of sets. If using --fast, perhaps try again without it. ***
The VCF files need to have columns #CHROM, POS … etc. The configuration script checks each VCF file for proper columns and will print a warning if any are missing or wrong. However, it does not halt execution and will include any malformed column VCF files and attempt to process them regardless. The main script may finish execution with the malformed VCF, but the output may be perturbed or useless.
File "mucor/inputs.py", in parse_VarScan
position = int(row[fieldId['POS']])
KeyError: 'POS'
Users may not supply vcf files that have inconsistent 'effect' and/or 'functional consequence' for the same variant. Presumably, if the variant has the same location and reference and alternate allele, the effect and functional consequences would be predicted to be the same. The problem lies in collapsing mutations in the variant_details output type(s). This issue may arise when different platforms or pipelines are used for samples containing a common variant within a single run of mucor. The current solution is to annotate all VCF files with the same functional consequence decoration software.
Check the issue tracker at the github repository. Please report bugs and request new features there for better tracking.
James S. Blachly, MD
[email protected]