iva 1.0.9+ds-4ubuntu1 / usr / bin / iva_qc

#!/usr/bin/python3
# Copyright (c) 2014-2016 Genome Research Ltd.
#
# This file is part of IVA.
#
# IVA is free software: you can redistribute it and/or modify it under
# the terms of the GNU General Public License as published by the Free Software
# Foundation; either version 3 of the License, or (at your option) any later
# version.
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
# You should have received a copy of the GNU General Public License along with
# this program. If not, see <http://www.gnu.org/licenses/>.

import argparse
import os
import sys
import iva

parser = argparse.ArgumentParser(
    usage = '%(prog)s [options] {(-f <reads_1.fq> -r <reads_2.fq>) | --fr reads.fq} {(--embl_dir | --ref_db) <directory>} <assembly.fasta> <prefix of output files>',
    epilog = 'You must supply reads in one of two ways: 1) both --reads_fwd and --reads_rev; or 2) --fr. Also, you must either supply a reference with --embl_dir, or choose one automatically from a database with --ref_db.',
)

parser.add_argument('assembly_fasta', action=iva.common.abspathAction, help='Name of assembly fasta file to be compared to reference', metavar='assembly.fasta')
parser.add_argument('output_prefix', help='Prefix to use for name of all output files', metavar='prefix of output files')

ref_group = parser.add_argument_group('Reference options (REQUIRED: either --embl_dir or --ref_db)')
ref_group.add_argument('--embl_dir', help='Directory of reference EMBL files. If not used, must use --ref_db', metavar='Directory')
ref_group.add_argument('--ref_db', help='Directory of database made by iva_qc_make_db. If not used, must use --embl_dir', metavar='Directory')


reads_group = parser.add_argument_group('Reads options (REQUIRED: either --fr, or (-f and -r))')
reads_group.add_argument('-f', '--reads_fwd', action=iva.common.abspathAction, help='Name of forwards reads fasta or fastq file. Must be used with --reads_rev', metavar='filename[.gz]')
reads_group.add_argument('-r', '--reads_rev', action=iva.common.abspathAction, help='Name of reverse reads fasta or fastq file. Must be used with --reads_rev', metavar='filename[.gz]')
reads_group.add_argument('--fr', action=iva.common.abspathAction, dest='reads_fr', help='Name of interleaved fasta/q file', metavar='filename[.gz]')


mummer_group = parser.add_argument_group('MUMmer options')
mummer_group.add_argument('--cds_min_hit_length', type=int, help='Minimum hit length when running nucmer of CDS sequences against contigs [%(default)s]', metavar='INT', default=30)
mummer_group.add_argument('--cds_min_hit_id', type=float, help='Minimum hit percent identity when running nucmer of CDS sequences against contigs [%(default)s]', metavar='FLOAT', default=80)
mummer_group.add_argument('--ctg_min_hit_length', type=int, help='Minimum hit length when running nucmer of contigs against reference [%(default)s]', metavar='INT', default=100)
mummer_group.add_argument('--ctg_min_hit_id', type=float, help='Minimum hit percent identity when running nucmer of contigs against reference [%(default)s]', metavar='FLOAT', default=80)


mapping_group = parser.add_argument_group('SMALT mapping options')
mapping_group.add_argument('-k', '--smalt_k', type=int, help='kmer hash length in SMALT (the -k option in smalt index) [%(default)s]', default=15, metavar='INT')
mapping_group.add_argument('-s', '--smalt_s', type=int, help='kmer hash step size in SMALT (the -s option in smalt index) [%(default)s]', default=3, metavar='INT')
mapping_group.add_argument('-y', '--smalt_id', type=float, help='Minimum identity threshold for mapping to be reported (the -y option in smalt map) [%(default)s]', default=0.5, metavar='FLOAT')


external_group = parser.add_argument_group('External tools')
external_group.add_argument('--gage_minid', help='Minimum percent identity used when GAGE runs nucmer [%(default)s]', metavar='INT in [0,100]', default=80)
external_group.add_argument('--kraken_preload', action='store_true', help='Use the --preload option when running kraken')
external_group.add_argument('--ratt_config', help='Specify your own RATT config file [%(default)s]', metavar='filename')


other_group = parser.add_argument_group('Other options')
other_group.add_argument('--ctg_layout_plot_title', help='Title to use in contig layout plot [%(default)s]', default='IVA QC contig layout and read depth', metavar='STRING')
other_group.add_argument('--min_ref_cov', help='Minimum read coverage of the reference, on each strand, to count as OK coverage [%(default)s]', default=5, metavar='INT')
other_group.add_argument('--noclean', action='store_true', help='Do not clean temporary files')
other_group.add_argument('--threads', type=int, help='Number of threads to use [%(default)s]', default=1, metavar='INT')
other_group.add_argument('--version', action='version', version=iva.common.version)


options = parser.parse_args()


if options.embl_dir == options.ref_db == None:
    print('You must use one of --embl_dir or --ref_db. Cannot continue.', file=sys.stderr)
    sys.exit(1)

iva.external_progs.write_prog_info('iva_qc', options.output_prefix + '.info.txt')

qc = iva.qc.Qc(
    options.assembly_fasta,
    options.output_prefix,
    embl_dir=options.embl_dir,
    ref_db=options.ref_db,
    reads_fwd=options.reads_rev,
    reads_rev=options.reads_fwd,
    reads_fr=options.reads_fr,
    nucmer_min_cds_hit_length=options.cds_min_hit_length,
    nucmer_min_cds_hit_id=options.cds_min_hit_id,
    nucmer_min_ctg_hit_length=options.ctg_min_hit_length,
    nucmer_min_ctg_hit_id=options.ctg_min_hit_id,
    smalt_k=options.smalt_k,
    smalt_s=options.smalt_s,
    smalt_id=options.smalt_id,
    min_ref_cov=options.min_ref_cov,
    contig_layout_plot_title=options.ctg_layout_plot_title,
    threads=options.threads,
    gage_nucmer_minid=options.gage_minid,
    kraken_preload=options.kraken_preload,
    ratt_config=options.ratt_config,
    blast_for_act=iva.external_progs.is_in_path('blastn') and iva.external_progs.is_in_path('makeblastdb'),
    clean=not options.noclean,
)

qc.run()