Extracting allele tables using pmotools-python

pmotools-python extract_allele_table

To extract allele table information from a PMO the command line interactive script with pmotools-python extract_allele_table can be used

pmotools-python extract_allele_table
Required arguments
- --file - the PMO file to extract from
- --output - the output stub of the files to be created
Optional arguments

By default only 3 fields are extracted by this extractor, 1) sampleID (library_sample_sample_name), 2) locus (target_name), and 3) allele (microhaplotype_id) with those default column names. This can be controlled by --default_base_col_names and if you supply 3 comma separated values you can change the default header.

You can also add to the table any values from the other portions of the PMO file by using the following arguments

adding fields arguments
- --specimen_info_meta_fields - Meta Fields if any to include from the specimen table
- --library_sample_info_meta_fields - Meta Fields if any to include from the library_sample table
- --microhap_fields - additional optional fields from the detected microhaplotype object to include
- --representative_haps_fields - additional optional fields from the detected representative object to include

Other optional arguments have to do with the ouput file over writing and delimiter being used, use -h to see all arguments

Code

pmotools-python extract_allele_table -h

usage: pmotools-python extract_allele_table [-h] --file FILE
                                            [--jsonschema JSONSCHEMA]
                                            [--delim DELIM] --output OUTPUT
                                            [--overwrite]
                                            [--allele_freqs_output ALLELE_FREQS_OUTPUT]
                                            [--specimen_info_meta_fields SPECIMEN_INFO_META_FIELDS]
                                            [--library_sample_info_meta_fields LIBRARY_SAMPLE_INFO_META_FIELDS]
                                            [--microhap_fields MICROHAP_FIELDS]
                                            [--representative_haps_fields REPRESENTATIVE_HAPS_FIELDS]
                                            [--default_base_col_names DEFAULT_BASE_COL_NAMES]

options:
  -h, --help            show this help message and exit
  --file FILE           PMO file
  --jsonschema JSONSCHEMA
                        the jsonschema to check the PMO against
  --delim DELIM         the delimiter of the input text file, examples
                        tab,comma
  --output OUTPUT       Output allele table file name path
  --overwrite           If output file exists, overwrite it
  --allele_freqs_output ALLELE_FREQS_OUTPUT
                        if also writing out allele frequencies, write to this
                        file
  --specimen_info_meta_fields SPECIMEN_INFO_META_FIELDS
                        Meta Fields if any to include from the specimen table
  --library_sample_info_meta_fields LIBRARY_SAMPLE_INFO_META_FIELDS
                        Meta Fields if any to include from the library sample
                        table
  --microhap_fields MICROHAP_FIELDS
                        additional optional fields from the detected
                        microhaplotype object to include
  --representative_haps_fields REPRESENTATIVE_HAPS_FIELDS
                        additional optional fields from the detected
                        representative object to include
  --default_base_col_names DEFAULT_BASE_COL_NAMES
                        default base column names, must be length 3

The python code for extract_allele_table script is below

Code

pmotools-python/src/pmotools/scripts/extractors_from_pmo/extract_allele_table.py

#!/usr/bin/env python3

import argparse
import json
import os

from pmotools.pmo_engine.pmo_reader import PMOReader
from pmotools.utils.small_utils import Utils
from pmotools.pmo_engine.pmo_checker import PMOChecker
from pmotools.pmo_engine.pmo_processor import PMOProcessor
from pmotools import __version__ as __pmotools_version__


def parse_args_extract_for_allele_table():
    parser = argparse.ArgumentParser()
    parser.add_argument("--file", type=str, required=True, help="PMO file")
    parser.add_argument(
        "--jsonschema",
        default=os.path.join(
            os.path.dirname(
                os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
            ),
            "schemas/",
            f"portable_microhaplotype_object_v{__pmotools_version__}.schema.json",
        ),
        type=str,
        required=False,
        help="the jsonschema to check the PMO against",
    )

    parser.add_argument(
        "--delim",
        default="tab",
        type=str,
        required=False,
        help="the delimiter of the input text file, examples tab,comma",
    )
    parser.add_argument(
        "--output", type=str, required=True, help="Output allele table file name path"
    )
    parser.add_argument(
        "--overwrite", action="store_true", help="If output file exists, overwrite it"
    )
    parser.add_argument(
        "--allele_freqs_output",
        type=str,
        help="if also writing out allele frequencies, write to this file",
    )

    parser.add_argument(
        "--specimen_info_meta_fields",
        type=str,
        required=False,
        help="Meta Fields if any to include from the specimen table",
    )
    parser.add_argument(
        "--library_sample_info_meta_fields",
        type=str,
        required=False,
        help="Meta Fields if any to include from the library sample table",
    )
    parser.add_argument(
        "--microhap_fields",
        type=str,
        required=False,
        help="additional optional fields from the detected microhaplotype object to include",
    )
    parser.add_argument(
        "--representative_haps_fields",
        type=str,
        required=False,
        help="additional optional fields from the detected representative object to include",
    )
    parser.add_argument(
        "--default_base_col_names",
        type=str,
        required=False,
        default="library_sample_name,target_name,mhap_id",
        help="default base column names, must be length 3",
    )

    return parser.parse_args()


def extract_for_allele_table():
    args = parse_args_extract_for_allele_table()
    compressed_output = (
        "." not in args.output and args.file.endswith(".gz")
    ) or args.output.endswith(".gz")

    output_delim, output_extension = Utils.process_delimiter_and_output_extension(
        args.delim, gzip=compressed_output
    )

    allele_per_sample_table_out_fnp = (
        args.output
        if "STDOUT" == args.output
        else Utils.appendStrAsNeeded(args.output, output_extension)
    )
    Utils.inputOutputFileCheck(
        args.file, allele_per_sample_table_out_fnp, args.overwrite
    )

    allele_freq_output = ""
    if args.allele_freqs_output is not None:
        allele_freq_output = Utils.appendStrAsNeeded(
            args.allele_freqs_output, output_extension
        )
        Utils.inputOutputFileCheck(args.file, allele_freq_output, args.overwrite)

    pmodata = PMOReader.read_in_pmo(args.file)
    with open(args.jsonschema, "r") as f:
        schema_dict = json.load(f)
        checker = PMOChecker(schema_dict)
        # make sure PMO is valid
        checker.validate_pmo_json(pmodata)

    if args.specimen_info_meta_fields is not None:
        args.specimen_info_meta_fields = Utils.parse_delimited_input_or_file(
            args.specimen_info_meta_fields, ","
        )
    if args.microhap_fields is not None:
        args.microhap_fields = Utils.parse_delimited_input_or_file(
            args.microhap_fields, ","
        )
    if args.library_sample_info_meta_fields is not None:
        args.library_sample_info_meta_fields = Utils.parse_delimited_input_or_file(
            args.library_sample_info_meta_fields, ","
        )
    if args.representative_haps_fields is not None:
        args.representative_haps_fields = Utils.parse_delimited_input_or_file(
            args.representative_haps_fields, ","
        )

    allele_table = PMOProcessor.extract_alleles_per_sample_table(
        pmodata,
        additional_specimen_info_fields=args.specimen_info_meta_fields,
        additional_library_sample_info_fields=args.library_sample_info_meta_fields,
        additional_microhap_fields=args.microhap_fields,
        additional_representative_info_fields=args.representative_haps_fields,
        default_base_col_names=args.default_base_col_names.split(","),
    )
    with Utils.smart_open_write(allele_per_sample_table_out_fnp) as f:
        allele_table.to_csv(f, sep=output_delim, index=False)

    if args.allele_freqs_output is not None:
        allele_freqs = PMOProcessor.extract_allele_counts_freq_from_pmo(pmodata)
        with Utils.smart_open_write(allele_freq_output) as f:
            allele_freqs.to_csv(f, sep=output_delim, index=False)


if __name__ == "__main__":
    extract_for_allele_table()

Can download example PMOs here

Code

wget https://plasmogenepi.github.io/PMO_Docs/format/moz2018_PMO.json.gz 

wget https://plasmogenepi.github.io/PMO_Docs/format/PathWeaverHeome1_PMO.json.gz

Code

mkdir -p example
cd example 

# default 
pmotools-python extract_allele_table --file ../../format/moz2018_PMO.json.gz --output extraction --overwrite

Code

cd example 

# changing default column names  
pmotools-python extract_allele_table --file ../../format/moz2018_PMO.json.gz --output extraction --overwrite --default_base_col_names sample,target,hapid

Changing the output file delimiter

Code

cd example 

pmotools-python extract_allele_table --file ../../format/moz2018_PMO.json.gz --output extraction --overwrite --delim ,

Adding on additional columns from the specimen_infos section

Code

cd example 

# adding other PMO fields 
pmotools-python extract_allele_table --file ../../format/moz2018_PMO.json.gz --output STDOUT --specimen_info_meta_fields collection_country,collection_date | head

bioinformatics_run_name library_sample_name target_name mhap_id collection_country  collection_date
Mozambique2018-SeekDeep 8025875168  t99 3   NA  NA
Mozambique2018-SeekDeep 8025875168  t99 0   NA  NA
Mozambique2018-SeekDeep 8025875168  t97 0   NA  NA
Mozambique2018-SeekDeep 8025875168  t97 3   NA  NA
Mozambique2018-SeekDeep 8025875168  t94 1   NA  NA
Mozambique2018-SeekDeep 8025875168  t94 0   NA  NA
Mozambique2018-SeekDeep 8025875168  t90 2   NA  NA
Mozambique2018-SeekDeep 8025875168  t90 0   NA  NA
Mozambique2018-SeekDeep 8025875168  t85 0   NA  NA

Can continue to add on more columns from other sections

Code

cd example 

# adding other PMO fields including seq field 
pmotools-python extract_allele_table --file ../../format/moz2018_PMO.json.gz --output STDOUT --specimen_info_meta_fields collection_country,collection_date --representative_haps_fields seq | head

bioinformatics_run_name library_sample_name target_name mhap_id collection_country  collection_date seq
Mozambique2018-SeekDeep 8025875168  t99 3   NA  NA  ATGGAAAAATGGAATATGAAGTATTAAGTGATGATAAAATAGTGTATGAAAATATACAACATGATTTATTAAAAACAATAGAAGATGGTGAAGAAATGTTAAAAGGAACTGAAAGGAAGGATAATATAGATATACTGAGGACTCCTGGAAGGGGAGAATATAATATGTGGTCTACTTCTGGACTAGGGTTCTATGAAT
Mozambique2018-SeekDeep 8025875168  t99 0   NA  NA  ATGGAAAAATGGAATATGAAGTATTAAGTGATGATAAAATAGTGTATGAAAATATACAACATGATTTATTAAAAACAATAGAAGATGATGACGAAATGTTAAAAGGAACTGAAAGGAAGGATAATATAGATATACTGAGGACTCCTGGAAGGGGAGAATATAATATGTGGTCTACTTCTGGACTAGGGTTCTATGAAT
Mozambique2018-SeekDeep 8025875168  t97 0   NA  NA  TAAACACCAGTACCATTTTTTTCTGATAAATTAATATTTTTTTGTATAACATCATATTTATCCCTTTTCGTGGTAAGTGCAGTATCCTGTTTTATTATTATATTATCGAATTCATCATGGTGTATATTTCTTTCT
Mozambique2018-SeekDeep 8025875168  t97 3   NA  NA  TAATCACCAGTACCATTTTTTTCTGATAAATTAATATTTTTTTGTATAACATCATATTTATCCCTTTTCGTGGTAAGTGCAGTATCCTGTTTTATTATTATATTATCGAATTCATCATGGTGTATATTTCTTTCA
Mozambique2018-SeekDeep 8025875168  t94 1   NA  NA  TATTAAAACTTTTTTTTCTTTCTGTAAAGTTTGTACATTATGTTTTGATGAGTTTTTATTATCTTCATAAAACTTTATATATTTATAAAAATTATTTTGTATAAAATCATTTAATAAAGGTAACATAATTTTTTTAGCTTGATTCAATTCACTACATGAATGTATAT
Mozambique2018-SeekDeep 8025875168  t94 0   NA  NA  TATTAAAACTTTTTTTTCTTTCTGTAAAGTTTGTACATTATGTTTTGATGAGTTTTGATTATCTTCATAAAACTTTATATATTTATAAAAATTATTTTGTATAAAATCATTTAATAAAGGTAACATAATTTTTTTAGCTTGATTCAATTCACTACATGAATGTATAT
Mozambique2018-SeekDeep 8025875168  t90 2   NA  NA  TTACAATGTTCTTCGCATTCGAAATTTTTTTCAGGATTACTTGAAAAGCCTTGTGGACAATTACAATATTCATATCCATGAGCATTCTTACAAACACCTTTTCCACAATTTAAAAAACATTTTTCTTCATTTAA
Mozambique2018-SeekDeep 8025875168  t90 0   NA  NA  TTACAATGTTCTTCGCATTCGAAATTTTTTTCAGGATTACTTGAAAAGCCTTCTGGACAATTACAATATTCATATCCATGAACATTCTTACAAACACCTTTTCCACAATTTAAAAAACATTTTTCTTCATTTAA
Mozambique2018-SeekDeep 8025875168  t85 0   NA  NA  AACATTTTTTTAACATCTTTACCTTTTTGACTTGGTTCTTCATCATAATTCTGTTGTTCTGCAGAATCAGCATTTACTTCAGTTTCTTCTTTATTTTGAAAAGTGTTTGATTCTACATGAGAATTGGAAGATGAACGTCTATGTTTTACTTCTGTATAACTAGTACGTTCTCCAGTATGATGAGCCTTAAGGTTTACATCTTCAGTTTCTT

Creating output for MOIRE

MOIRE is a program that can be used to estimate COI and other population estimates from a population. See it’s github for full usage.

Code

mkdir -p example
cd example 

# default table is all moire needs 
pmotools-python extract_allele_table --file ../../format/moz2018_PMO.json.gz --output extraction --overwrite

Code

df <- read.csv("example/extraction_allele_table.tsv", sep = "\t")

data <- load_long_form_data(df)

# With data in appropriate format, run MCMC as follows
mcmc_results <- moire::run_mcmc(data, is_missing = data$is_missing)

Creating output for dcifer

dcifer is a program that can estimate IBD even from mixed infections. See it’s github for full usage

Code

mkdir -p example
cd example 

# default 
pmotools-python extract_allele_table --file ../../format/moz2018_PMO.json.gz --output extraction --overwrite --delim ,

# dcifer can calculate allele frequencies if not provided or you can have extract_allele_table write them as well 
pmotools-python extract_allele_table --file ../../format/moz2018_PMO.json.gz --output extraction --overwrite --allele_freqs_output allele_freqs_extraction --delim ,

Code

dsmp <- readDat("example/extraction_allele_table.csv", svar = "sampleID", lvar = "locus", avar = "allele")

lrank <- 2
coi   <- getCOI(dsmp, lrank = lrank)

afreq <- calcAfreq(dsmp, coi, tol = 1e-5) 

dres0 <- ibdDat(dsmp, coi, afreq, pval = TRUE, confint = TRUE, rnull = 0, 
                alpha = 0.05, nr = 1e3)