Skip to content

Heteroallelic Genotyping tools

Jump to bottom
Olivo Miotto edited this page Jun 13, 2017 · 11 revisions

Single Heteroallelic Genotyping Tool

Invoking the genotyping tool for a single sample:

 java -Xms512m -Xmx2000m 'org.cggh.bam.heteroallelic.HeteroallelicAnalysis$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 “Heteroallelic configuration file”
  • 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 Heteroallelic Genotyping Tool

Invoking the genotyping tool for multiple samples:

java -Xms512m -Xmx2000m 'org.cggh.bam.heteroallelic.HeteroallelicAnalysis$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 “Heteroallelic configuration file”

  • 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. These files are placed in an output folder <OUT_DIR>/<SUB> where <SUB> is a string consisting of the first 4 letters of the sample name. For example, file PA0001-C.alleles.tab will be written to folder <OUT_DIR>/PA00. This simple file hashing scheme, which is particularly suited for the MalariaGEN naming scheme, prevents thousands of files accumulating in a single folder, which may cause indexing problems especially on Linux.

  • <SAMPLE_NAME>.<LOCUS>.mutations.tab contains one line (record) for each codon within the locus in which a nonsynonymous mutation is found in this sample. If there are multiple mutations, the file will contain multiple records for the locus, one per mutation. 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 locus, as specified in the configuration file
    • Codon The number of the codon where the mutation is found
    • Call The call for this mutation ("MU" is a homozygous mutant, "HE" a heterozygous one)
    • Mutation The mutation observed (e.g. C580Y represents a C -> Y mutation at codon 580)
    • TotalReadCount The read coverage of this codon
    • MutantReadCount The number of codon reads supporting the mutation call
    • MutantReadProp The proportion of codon reads supporting the mutation call

  • <SAMPLE_NAME>.<LOCUS>.calls.tab contains one line (record) containing the composite call resulting from aggregating all mutation calls within the locus found in this sample. The 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 locus, as specified in the configuration file
    • Call The call for this sample ("WT" is wild-type, "MU" is a homozygous mutant, "HE" a heterozygous mutant, and "MI" is a missing call)
    • Mutation The mutation(s) observed, as a comma-separated list. Het mutations are marked by an asterisk
    • MissingCodonCallsProp The proportion of codons within the locus where a call could not be made (coverage too low)
    • MedianReadCount The median read count at codons across the locus
    • MeanReadCount The mean read count at codons across the locus

  • <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

Genotype Aggregation Tool

The Genotype 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.heteroallelic.HeteroallelicAnalysis$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:

  • CONFIG_FILE This is a configuration file, defining loci and targets, in Java .properties format. The format is specified in section “Heteroallelic configuration file”
  • 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.

Genotype Aggregation Tool Outputs

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

  • 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.

#Heteroallelic Configuration File The configuration file for the Heteroallelic genotyping tools is a file in Java .properties format (a text-based file of name/value pairs). The following properties are specified:

  • grc.loci A comma-separated list of locus names that will be used in the analysis

  • grc.remapAllReads If set to 'true' all BAM-mapped reads will be remapped using the anchors (see below) instead of using the BAM-assigned position. Default (when property is missing) is 'false'.

  • grc.locus.<LOCUS_NAME>.region The span of the locus, over which mapped reads will be analyzed, in the form <chr>:<startPos>-<endPos> (must be specified for each locus)

  • grc.locus.<LOCUS_NAME>.analyzeUnmappedReads If set to 'true' the unmapped reads in the BAM will be searched for anchor sequences and remapped. Default (when property is missing) is 'false'. This has a performance impact, so it should only be set to 'true' at loci thought to be problematic for alignments.

  • grc.locus.<LOCUS_NAME>.targets A comma-separated list of targets to be genotyped. Each target is specified in the format <targetName>@<targetStartPos>-<targetEndPos>

  • grc.locus.<LOCUS_NAME>.anchors A comma-separated list of the anchors (regex) that will be used to match the reads. Each anchor is specified in the format <anchorStartPos>@<regex>

The following is an example of a configuration file that will genotype the core (72-76) haplotype of pfcrt, and the arps10-127 mutation. Note that the arps10 anchors will match alternative alleles at some positions (“[CA]” will match either C or A), and skip testing some position (“.” is a wildcard)

grc.loci=crt_core,arps10
# CRT
grc.locus.crt_core.region=Pf3D7_07_v3:403500-403800
grc.locus.crt_core.targets=crt_72-76@403612-403626
grc.locus.crt_core.anchors=403593@TATTATTTATTTAAGTGTA,403627@ATTTTTGCTAAAAGAAC
# ARPS10
grc.locus.arps10.region=Pf3D7_14_v3:2480900-2481200
grc.locus.arps10.targets=arps10_127@2481070-2481072
grc.locus.arps10.anchors=2481045@ATTTAC[CA]TTTTTGCGATCTCCCCAT...[GC],2481079@GACAGT[AC]G[AG]GA[GA]CAATTCGAAATAAAAC 

Clone this wiki locally