Skip to content

Home

Jump to bottom
Olivo Miotto edited this page Oct 22, 2016 · 43 revisions

Welcome to the documentation for the GeneticReportCard programs.

Introduction

This project contains a set of programs that are used to inspect BAM files and derive genotypes relevant to the Genetic Report Cards (GRC). The tools have been initially developed for the analysis of aligned P. falciparum Illumina paired short reads files, but will be extended to other organisms in the future. Most of the organism-specific parameters are encoded in configuration files for flexibility.

There are two main set of tools at present:

  • The GRC tools produce genotypes for the drug-resistance loci that are included in the GRC. These genotypes are (currently) expressed in terms of amino acid alleles
  • The species tools search for reads likely to originate from species other than P. falciparum, to determine the presence of co-infecting strains in the samples

Both tools take as primary input BAM files from alignments against the V3 P. falciparum 3D7 reference genome. At a future stage, the choice of reference genome will be configurable. The tools can be invoked either on an individual sample, or on a set of samples. Genotyping tools are provided to produce one file per sample containing the genotyping data. There is also an aggregation tool that merges all results, producing a genotype file for a whole sample set.

Concepts

The tools share a common model for how they determine their results. The process is controlled by a Task Configuration File, which specifies a set of loci where genotyping is to be performed. A locus is a region of the reference genome upon which reads covering the positions to be genotyped are likely to be mapped. Typically, for genotyping a single codon, we may specify a locus which contains all positions within +/- 200bp on either side.

For each locus we specify one or more anchors, i.e. nucleotide sequence patterns that must be matched in each read to be used for genotyping. An anchor is specified as a regular expression (regex), and the position of the first nucleotide in the expression. The tools will attempt to match the anchor regex in both reads that are mapped to the locus, and also unmapped reads. If an anchor is matched, the read will be aligned at the anchor’s starting position for further processing. No gaps are introduced in this alignment step.

If genotypes are to be produced at some position in the locus (this is the case for the GRC tools, but not for the species tools which simply match alleles), then one or more targets must be specified for the locus. A target is specified as a names interval of nucleotide positions, currently assumed to be codons (hence the target length must be a multiple of 3). The tools will establish the codon(s) at the target location for all reads aligned against the anchor(s) in order to determine the genotype for the sample.

In the species tools, targets are used somewhat differently: for each target (which need not contain an exact number of codons), a set of alleles is associated with a class (in this case, the name of a species) which is assigned to the reads. In other words, here the targets are used for exact sequence matching against a known set of alleles.

Building the tools

Clone the GeneticReportCard repository. This contains two projects:

  • AnalysisCommon contains library that are shared by multiple projects (and which will probably be packaged separately eventually).
  • SequencingReadsAnalysis contains the sources for the reads analysis tools (grc and species)

Ensure that Java JRE 1.8 or higher is installed; test using

java –version

Set working directory to the root folder of SequencingReadsAnalysis project and type

ant build

The Ant task should build both projects.

Running the tools

Running the tools requires Java JRE 1.8 or higher (see above). The tools use several third-party libraries, which are stored in the /lib folders of each project. These must be added to the Java class path, as must the /bin folders of both projects. The following is a set of commands that sets the classpath:

GRCC =/path/to/clone/of/GeneticReportCard/AnalysisCommon
GRCA =/path/to/clone/of/GeneticReportCard/SequencingReadsAnalysis
CLASSPATH=$GRCA/bin:$GRCC/bin:\
$GRCC/lib/commons-logging-1.1.1.jar:\
$GRCA/lib/apache-ant-1.8.2-bzip2.jar:\
$GRCA/lib/commons-compress-1.4.1.jar:\
$GRCA/lib/commons-jexl-2.1.1.jar:\
$GRCA/lib/htsjdk-2.1.0.jar:\
$GRCA/lib/ngs-java-1.2.2.jar:\
$GRCA/lib/snappy-java-1.0.3-rc3.jar:\
$GRCA/lib/xz-1.5.jar

You need to ensure that the Java tools can allocate sufficient memory. Unfortunately I have not been able to do rigorous testing, and memory usage seems to be largely determined by the Broad BAM processing libraries. An allocation of 2GB seems to work; set this with

java -Xms512m -Xmx2000m <class> <params...>

Genetic Report Card genotyping

Single Sample Genotyping Tool

Invoking the genotyping tool for a single sample:

 java -Xms512m -Xmx2000m 'org.cggh.bam.grc.GrcAnalysis$SingleSample' \
     <CONFIG_FILE> <SAMPLE_NAME> <BAM_FILE> <CHR_MAP_NAME> \
     <REF_FASTA_FILE> <CHR_MAP_FILE> <OUT_DIR>

Note that the qualified class name HAS to be between single quotes. The parameters are as follows:

  • CONFIG_FILE This is a configuration file, defining loci and targets, in Java .properties format. The format is specified in section “GRC Configuration”
  • SAMPLE_NAME The identifier of the sample, e.g. “PA0001-C”
  • BAM_FILE The path of the BAM file for the sample
  • CHR_MAP_NAME The name of the chromosome name map to be used for this sample (e.g. “default”). See section “Chromosome Maps”
  • REF_FASTA_FILE The path to a FASTA file containing the reference sequence used for alignment.
  • CHR_MAP_FILE The path to a file specifying the chromosome name maps available. See section “Chromosome Maps”.
  • OUT_DIR The path to the folder where the output files will be written to.

Multiple Samples Genotyping Tool

Invoking the genotyping tool for multiple samples:

java -Xms512m -Xmx2000m 'org.cggh.bam.grc.GrcAnalysis$MultiSample' \
     <CONFIG_FILE> <SAMPLE_LIST_FILE> \
     <REF_FASTA_FILE> <CHR_MAP_FILE> <OUT_DIR>

Note that the qualified class name HAS to be between single quotes. The parameters are as follows:

  • CONFIG_FILE This is a configuration file, defining loci and targets, in Java .properties format. The format is specified in section “GRC Configuration”

  • SAMPLE_LIST_FILE The path of tab-separated text file, with one line (record) per sample to be processed. Each record contains three fields:

    • SAMPLE_NAME The identifier of the sample, e.g. “PA0001-C”
    • BAM_FILE The path of the BAM file for the sample
    • CHR_MAP_NAME The name of the chromosome name map to be used for this sample (e.g. “default”). See section “Chromosome Maps”

  • REF_FASTA_FILE The path to a FASTA file containing the reference sequence used for alignment.

  • CHR_MAP_FILE The path to a file specifying the chromosome name maps available. See section “Chromosome Maps”.

  • OUT_DIR The path to the folder where the output files will be written to.

Genotyping Tools Outputs

The genotyping tools output two tab-separated text files for each sample into the output folder:

  • <SAMPLE_NAME>.alleles.tab contains one line (record) for each allele genotyped at each target. If at a given target multiple alleles are observed, the file will contain multiple records for the target, one per allele. Each record has the following fields:

    • Num A sequential record number (can be ignored)
    • Sample The identifier of the sample, e.g. “PA0001-C”
    • Locus The name of the target’s locus, as specified in the configuration file
    • Target The name of the target, as specified in the configuration file
    • Allele Allele (nucleotide sequence) found
    • Amino Amino acid translation of allele
    • Count Number of supporting reads for this allele

  • <SAMPLE_NAME>.locusCoverage.tab contains one line (record) for each locus, detailing the coverage (number of reads used in the analysis). Each record has the following fields:

    • Num A sequential record number (can be ignored)
    • Sample The identifier of the sample, e.g. “PA0001-C”
    • Locus The name of the target’s locus, as specified in the configuration file
    • Aligned The number of reads that were aligned by anchor matching
    • Misaligned The number of reads where an anchor was matched, but a high number of point differences (currently >10) was found between the read sequence and the consensus sequence of the aligned reads. These reads are discarded from further analysis. A high proportion of misaligned reads may indicate the presence of indels or repeat polymorphisms at the locus

Result Aggregation Tool

The Results Aggregation Tool uses the files produced by the Genotyping Tool and aggregates them, producing sampleset-wide outputs. The aggregation tool performs some statistics across the whole sample set, applying a filtering scheme similar to MinAlt: an allele is discarded if there is not a single sample in which the allele is homozygous, or in which it is supported by at least 10 reads. In addition, target alleles supported by singleton reads are disregarded.

Invoking the result aggregation tool:

java -Xms512m -Xmx2000m 'org.cggh.bam.grc.GrcAnalysis$MergeResults' \
     <CONFIG_FILE> <SAMPLE_LIST_FILE> \
     <REF_FASTA_FILE> <CHR_MAP_FILE> <OUT_DIR>

Note that the qualified class name HAS to be between single quotes. The parameters are as follows:

  • REF_FASTA_FILE
  • CONFIG_FILE This is a configuration file, defining loci and targets, in Java .properties format. The format is specified in section “GRC Configuration”
  • SAMPLE_LIST_FILE The path of tab-separated text file, with one line (record) per sample to be processed (note, only the first column is used here) . Each record contains three fields:
    • SAMPLE_NAME The identifier of the sample, e.g. “PA0001-C”
    • BAM_FILE The path of the BAM file for the sample
    • CHR_MAP_NAME The name of the chromosome name map to be used for this sample (e.g. “default”). See section “Chromosome Maps”
  • REF_FASTA_FILE The path to a FASTA file containing the reference sequence used for alignment.
  • CHR_MAP_FILE The path to a file specifying the chromosome name maps available. See section “Chromosome Maps”.
  • OUT_DIR The path to the folder where the output files will be written to.

Aggregation Tool Outputs

The aggregation tool outputs several tab-separated text files in output folder <OUT_DIR>

  • <SAMPLE_NAME>.alleles.tab

  • AllCallsBySample.tab This is the main output file, which summarizes the genotype calls for each sample. Each line (record) shows the genotypes for a sample, with the following fields:

    • Num A sequential record number (can be ignored)
    • Sample The identifier of the sample, e.g. “PA0001-C”
    • One column per target, containing the genotype. Each column has a header showing the target name, following by the 3D7 allele in square bracket. The value for each sample is the amino acid sequence of the allele; in heterozygous samples, a comma-separated list of all the alleles found is given.

  • AllCallsNrefBySample.tab Contains the same data as AllCallsBySample.tab, except that the reference allele is replaced by a dot (‘.’). As a result, this file is somewhat more informative at a glance if one wants to inspect the presence of non-reference alleles

  • AlleleSampleCount.final.<LOCUS>_<TARGET>.tab There is one of these files produced per target. It shows the counts of reads supporting each allele observed at this target, in each sample, with the following fields:

    • Num A sequential record number (can be ignored)
    • Sample The identifier of the sample, e.g. “PA0001-C”
    • One column per allele observed at this target in the sample set, containing the number of reads supporting the allele’s presence in the sample. The column header is the nucleotide sequence of the allele, followed by the amino acid sequence in square brackets.

  • AlleleSampleCount.beforeMinAlt.<LOCUS>_<TARGET>.tab This file is equivalent to the AlleleSampleCount.final.<LOCUS>_<TARGET>.tab file, but shows the alleles before they were filtered. It may be useful if trying to track anomalies.

  • AlleleStats.final.<LOCUS>_<TARGET>.tab There is one of these files produced per target. It shows some statistics for each of the observed alleles at this target in the sampleset. Each record represents an allele, with the following fields:

    • Num A sequential record number (can be ignored)
    • Allele The allele, represented as the nucleotide sequence of the allele, followed by the amino acid sequence in square brackets.
    • SampleCount The count of samples in which the allele was observed
    • MaxReads The highest read count seen for this allele in any one sample
    • MaxReadFraction The highest fraction of total reads seen for this allele in any one sample

  • AlleleStats.beforeMinAlt.<LOCUS>_<TARGET>.tab This file is equivalent to the AlleleStats.final.<LOCUS>_<TARGET>.tab file, but shows the statistics before the alleles were filtered. It may be useful if trying to track anomalies.

  • CallsBySample. final.<LOCUS>_<TARGET>.tab There is one of these files produced per target. It shows details of how each sample was called- essentially, a more detailed form of AllCallsBySample.tab. Each record represents a sample, with the following fields:

    • Num A sequential record number (can be ignored)
    • Sample The identifier of the sample, e.g. “PA0001-C”
    • Call The sample call: WT (reference allele), MU (mutant allele), HE (heterozygous), MI (missing call)
    • Alleles A comma-separated list of the amino acid alleles
    • AlleleReads A comma-separated list of the amino acid alleles; each allele is followed by a semicolon, and the count of reads supporting the allele
    • NtAlleles A comma-separated list of the nucleotide alleles
    • NtAlleleReads A comma-separated list of the nucleotide alleles; each allele is followed by a semicolon, and the count of reads supporting the allele

  • LocusCoverage.<LOCUS>.tab There is one of these files per locus. It shows counts of the reads that were used in the analysis of that locus. Each record represents a sample, with the following fields:

    • Num A sequential record number (can be ignored)
    • Sample The identifier of the sample, e.g. “PA0001-C”
    • Locus The name of the target’s locus, as specified in the configuration file
    • Aligned The number of reads that were aligned by anchor matching
    • Misaligned The number of reads where an anchor was matched, but a high number of point differences (currently >10) was found between the read sequence and the consensus sequence of the aligned reads.

Clone this wiki locally