RefNAAP_CLI.py

#!/usr/bin/env python

import sys
import os
import glob
import re
from datetime import date
import argparse
import subprocess
from pathlib import Path


def main():
    
    local_path = os.path.dirname(os.path.realpath(__file__))
    #print(local_path)
    data_path = f"{local_path}"
    scaffold_helper = f"{local_path}/scaffold_cutter.R"
    gapfixer_helper = f"{local_path}/gapfixer.R"
    now = date.today()
    home = str(Path.home())

    cli = argparse.ArgumentParser()
    cli.add_argument('-i', '--InputFolder', help="Folder containing barcoded fastq", required=True)
    cli.add_argument('-o', '--OutputFolder', help="Output Folder", required=False, default=f"{home}/refnaap_results/output_{now}")
    cli.add_argument('-r', '--RefFile', help="Reference File ", required=False, default=f'{local_path}/Americas2.fasta')

    cli.add_argument('--TopN', help="The top N reference sequences with the most depth are analyzed.", type=int, required=False, default=1)
    cli.add_argument('--MinCov', help="Amplicon regions need a minimum of this average coverage number", type=int, required=False, default=5)
    cli.add_argument('--Left', help="Bases to trim from left side of read", type=int, required=False, default=25)
    cli.add_argument('--Right', help="Bases to trim from right side of read", type=int, required=False, default=25)
    cli.add_argument('--Size', help="Filter reads less than this length", type=int, required=False, default=50)
    cli.add_argument('--threads', help="Number of threads. More is faster if your computer supports it", type=int, required=False, default=4)
    cli.add_argument('--verbose', help = "Keep Intermediate Files", required=False , action='store_true')
    cli.add_argument('--model', help="Basecall Model", required=False, type=str, default='r10_min_high_g303')
    args = cli.parse_args()


    #Run fastqc and multiqc on all the fastq/fastq.gz files in the folder
    subprocess.check_output(['python', local_path+'/fastqc_multiqc.py', '-i', args.InputFolder, '-o', args.OutputFolder+'/multiqc'])
    subprocess.check_output(['cp', args.OutputFolder+'/multiqc/multiqc_report.html', args.OutputFolder+'/multiqc_report.html'])

    #Interate over all the fastq/fastq.gz files
    files = sorted([f for f in glob.glob(args.InputFolder+"/**", recursive = True) if re.search(r'(.*)\.((fastq|fq)(|\.gz))$', f)])   
    print(files)
    OutputFolder = os.path.expanduser(args.OutputFolder)

    for i in range(0, len(files)):
        filec = files[i]

        base = os.path.splitext(os.path.basename(filec))[0]
        base = os.path.splitext(base)[0]
        print(base)

        filec2 = args.OutputFolder+'/'+"filtered/"+base+"_filtered.fastq"
        #Trim and filter the reads
        subprocess.check_output(['python', local_path+'/seqtk_sizefilter_trim.py', '-i', filec, '-o', filec2, '-l', str(args.Left), '-r', str(args.Right), '-s', str(args.Size)])


        #Get assembly
        subprocess.check_output(['python', local_path+'/refnaap_cli_helper.py', '-i', filec2, '-o', args.OutputFolder+'/assembly/'+base+"_assembly/", '-r', args.RefFile, '-t', str(args.threads), '--TopN', str(args.TopN), '--MinCov', str(args.MinCov)])
        subprocess.check_output(['cp', args.OutputFolder+'/assembly/'+base+"_assembly/final_scaffold.fasta", args.OutputFolder+"/"+base+"_final_scaffold.fasta"])



        print("progress: {}/{}".format(i+1, len(files)))

    

    if not args.verbose:
        subprocess.check_output(['rm', '-rf', args.OutputFolder+'/assembly'])
        subprocess.check_output(['rm', '-rf', args.OutputFolder+'/filtered'])
        subprocess.check_output(['rm', '-rf', args.OutputFolder+'/multiqc'])

if __name__ == "__main__":
    sys.exit(main())