#!/usr/bin/env julia
println()
println("Welcome to VIVA.")
println()
println("Loading packages:")
println()

println("1/5 ArgParse")
using ArgParse

println("2/5 VariantVisualization")
using Pkg #remove after VariantVisualization is deployed
Pkg.activate(@__DIR__) #remove after VariantVisualization is deployed
using VariantVisualization

println("3/5 PlotlyJS")
using PlotlyJS

println("4/5 GeneticVariation")
using GeneticVariation

println("5/5 DelimitedFiles")
using DelimitedFiles


println()
println("Finished loading packages!")
println()

function test_parse_main(ARGS::Vector{String})

    # initialize the settings (the description prints text when help is called)
    s = ArgParseSettings(

    description = "VIVA VCF Visualization Tool is a tool for creating publication quality plots of data contained within VCF files. For a complete description of features with examples read the docs here https://github.com/compbiocore/VariantVisualization.jl",
    suppress_warnings = true,
    epilog = "Thank you for using VIVA. Please submit any bugs to https://github.com/compbiocore/VariantVisualization.jl/issues "
    )

    @add_arg_table s begin
        "--vcf_file", "-f"
        help = "vcf filename in format: file.vcf"
        arg_type = String
        required = true

        "--output_directory", "-o"
        help =" function checks if directory exists and saves there, if not creates and saves here"
        arg_type = String
        default = "output"

        "--save_format", "-s"
        help = "file format you wish to save graphics as (eg. pdf, html, png). Defaults to html"
        arg_type = String
        default = "html"

        "--genomic_range", "-r"
        help = "select rows within a given chromosome range. Provide chromosome range after this flag in format chr4:20000000-30000000."
        arg_type = String

        "--pass_filter", "-p"
        help = "select rows with PASS in the FILTER field."
        action = :store_true

        "--positions_list", "-l"
        help = "select variants matching list of chromosomal positions. Provide filename of text file formatted with two columns in csv format: 1,2000345."
        arg_type = String

        "--group_samples", "-g"
        help = "group samples by common trait using user generated matrix key of traits and sample names following format guidelines in documentation. Provide file name of .csv file"
        nargs = 2
        arg_type = String

        "--select_samples"
        help = "select samples to include in visualization by providing tab delimited list of sample names (eg. samplenames.txt). Works for heatmap visualizations and numeric array generation only (not average dp plots)"
        arg_type = String

        "--heatmap", "-m"
        help = "genotype field to visualize (eg. genotype, read_depth, or 'genotype,read_depth' to visualize each separately)"
        arg_type = String
        default = "genotype,read_depth"

        "--y_axis_labels", "-y"
        help = "specify whether to label y-axis with all chromosome positions (options = positions / chromosome) separators. Defaults to chromosome separators."
        default = "chromosomes"
        arg_type = String

        "--x_axis_labels", "-x"
        help = "flag to specify whether to label x-axis with sample ids from vcf file. Defaults to FALSE."
        action = :store_true

        "--num_array", "-n"
        help = "flag to save numeric array of categorical genotype values or read depth values before heatmap plotting. Must be used with --heatmap set."
        action = :store_true

        "--heatmap_title", "-t"
        help = "Specify filename for heatmap with underscores for spaces."
        arg_type = String

        "--avg_dp"
        help = "visualize average read depths as line chart. Options: average sample read depth, average variant read depth, or both. eg. =sample, =variant, =sample,variant"

        "--save_remotely"
        help = "Save html support files online rather than locally so files can be shared between systems. Files saved in this way require internet access to open."
        action = :store_true

        #nargs = 2
        #metavar = ["avg_option", "markers_or_lines"]
        #default = ["variant,sample", "markers"]
        #default = "sample,variant" #turn on when plotly working

    end

    parsed_args = parse_args(s)
    # can turn off printing parsed args after development"
    #println("Parsed args:")

#activate block to show all argument keys

#=
    for (key,val) in parsed_args
        println("  $key  =>  $(repr(val))")
    end
    =#

    return parsed_args
end

parsed_args = test_parse_main(ARGS)

#filter vcf and load matrix of filtered vcf records

vcf_filename = (parsed_args["vcf_file"])
println("Reading $vcf_filename ...")
println()

reader = VCF.Reader(open(vcf_filename, "r"))
sample_names = get_sample_names(reader)

save_ext = parsed_args["save_format"]

remote_option = parsed_args["save_remotely"]

VariantVisualization.checkfor_outputdirectory(parsed_args["output_directory"])
output_directory=parsed_args["output_directory"]

#=
if parsed_args["avg_dp"] == nothing && parsed_args["heatmap"] == nothing
    number_records = nrecords((parsed_args["vcf_file"]))
    number_samples = nsamples((parsed_args["vcf_file"]))

    println("_______________________________________________")
    println()
    println("Summary Statistics of $(parsed_args["vcf_file"])")
    println()
    println("number of records: $number_records")
    println("number of samples: $number_samples")
    println("_______________________________________________")
    println()
    println("No plotting options specified. Plot data with --heatmap or --avg_dp_plot options")
    println()
end

=#

if parsed_args["x_axis_labels"] == true
    x_axis_label_option = true
else
    x_axis_label_option = false
end

#pass_filter
if parsed_args["pass_filter"] == true && parsed_args["genomic_range"] == nothing && parsed_args["positions_list"] == nothing
    println("Only pass filter is applied. Large vcf files with many PASS variants will take a long time to process and heatmap visualizations will lose resolution at this scale unless viewed in interactive html for zooming.")
    println()
    sub = VariantVisualization.io_pass_filter(vcf_filename)
    number_rows = size(sub,1)
    println("Selected $number_rows variants with Filter status: PASS")
    heatmap_input = "pass_filtered"
end

#chr_range
if parsed_args["genomic_range"] != nothing && parsed_args["pass_filter"] == false && parsed_args["positions_list"] == nothing
    sub = VariantVisualization.io_genomic_range_vcf_filter(parsed_args["genomic_range"],vcf_filename)
    number_rows = size(sub,1)
    println("Selected $number_rows variants within chromosome range: $(parsed_args["genomic_range"])")
    heatmap_input = "range_filtered"
end

#list
if parsed_args["positions_list"] != nothing && parsed_args["genomic_range"] == nothing && parsed_args["pass_filter"] == false
    sig_list =  load_siglist(parsed_args["positions_list"])
    sub = VariantVisualization.io_sig_list_vcf_filter(sig_list,vcf_filename)
    number_rows = size(sub,1)
    println("Selected $number_rows variants that match list of chromosome positions of interest")
    heatmap_input = "positions_filtered"
end

#pass_filter and chr_range and list
if parsed_args["pass_filter"] == true && parsed_args["genomic_range"] != nothing && parsed_args["positions_list"] != nothing
    sig_list =  load_siglist(parsed_args["positions_list"])
    sub = VariantVisualization.pass_genomic_range_siglist_filter(vcf_filename, sig_list, parsed_args["genomic_range"])
    number_rows = size(sub,1)
    println("Selected $number_rows variants with Filter status: PASS, that match list of chromosome positions of interest, and are within chromosome range: $(parsed_args["genomic_range"])")
end

#pass_filter and chr_range
if parsed_args["pass_filter"] == true && parsed_args["genomic_range"] != nothing && parsed_args["positions_list"] == nothing
    sub = VariantVisualization.pass_genomic_range_filter(reader, parsed_args["genomic_range"],vcf_filename)
    number_rows = size(sub,1)
    println("Selected $number_rows variants with Filter status: PASS and are within chromosome range: $(parsed_args["genomic_range"])")
end

#pass_filter and list
if parsed_args["pass_filter"] == true && parsed_args["genomic_range"] == nothing && parsed_args["positions_list"] != nothing
    sig_list =  load_siglist(parsed_args["positions_list"])
    sub = VariantVisualization.pass_siglist_filter(vcf_filename, sig_list)
    number_rows = size(sub,1)
    println("Selected $number_rows variants with Filter status: PASS and that match list of chromosome positions of interest")
end

#chr_range and list
if parsed_args["pass_filter"] == false && parsed_args["genomic_range"] != nothing && parsed_args["positions_list"] != nothing
    sig_list =  load_siglist(parsed_args["positions_list"])
    sub = VariantVisualization.genomic_range_siglist_filter(vcf_filename, sig_list, parsed_args["genomic_range"])
    number_rows = size(sub,1)
    println("Selected $number_rows variants that are within chromosome range: $(parsed_args["genomic_range"]) and that match list of chromosome positions of interest")
end

#no filters
if parsed_args["pass_filter"] == false && parsed_args["genomic_range"] == nothing && parsed_args["positions_list"] == nothing

    println("No filters applied. Large vcf files will take a long time to process and heatmap visualizations will lose resolution at this scale unless viewed in interactive html for zooming.")
    println()
    println("Loading VCF file into memory for visualization")

    sub = Array{Any}(undef,0)

    for record in reader
        push!(sub,record)
    end
    number_rows = size(sub,1)
    println("Selected $number_rows variants with no filters applied")


end

y_axis_label_option = parsed_args["y_axis_labels"]

number_rows = size(sub,1)

#convert to numeric array for plotting and generate/save heatmaps and scatter plots

if parsed_args["heatmap"] == "genotype"
    gt_num_array,gt_chromosome_labels = combined_all_genotype_array_functions(sub)

    if parsed_args["heatmap_title"] != nothing
        title = "Genotype_$(parsed_args["heatmap_title"])"
    else
        bn = Base.Filesystem.basename(parsed_args["vcf_file"])
        title = "Genotype_$bn"
    end

    chr_pos_tuple_list = generate_chromosome_positions_for_hover_labels(gt_chromosome_labels)

    chrom_label_info = VariantVisualization.chromosome_label_generator(gt_chromosome_labels[:,1])

    if length(parsed_args["group_samples"]) == 2

        if parsed_args["select_samples"] != nothing && length(parsed_args["group_samples"]) != 2
            println("Selecting samples listed in $(parsed_args["select_samples"])")
            gt_num_array,col_selectedcolumns = select_columns(parsed_args["select_samples"],
                                          gt_num_array,
                                          sample_names)
            sample_names = col_selectedcolumns



        elseif parsed_args["select_samples"] != nothing && length(parsed_args["group_samples"]) == 2
            println("Can not select samples with phenotype matrix provided. Please include same sample ids in phenotype matrix as in list of sample names to select.")

        end

        group_trait_matrix_filename=((parsed_args["group_samples"])[1])
        trait_to_group_by = ((parsed_args["group_samples"])[2])
        println()
        println("Grouping samples by $trait_to_group_by")
        println()

        ordered_num_array,group_label_pack,pheno,id_list,trait_labels = sortcols_by_phenotype_matrix(group_trait_matrix_filename, trait_to_group_by, gt_num_array, sample_names)

        if parsed_args["num_array"] == true
        save_numerical_array(ordered_num_array,sample_names,chr_pos_tuple_list,title,output_directory)
        end

        pheno_num_array,trait_label_array,chrom_label_info=add_pheno_matrix_to_gt_data_for_plotting(ordered_num_array,pheno,trait_labels,chrom_label_info,number_rows)

        graphic = VariantVisualization.genotype_heatmap_with_groups(pheno_num_array,title,chrom_label_info,group_label_pack,id_list,chr_pos_tuple_list,y_axis_label_option,trait_label_array,x_axis_label_option,number_rows)

    else

        if parsed_args["select_samples"] != nothing && length(parsed_args["group_samples"]) != 2
            println("Selecting samples listed in $(parsed_args["select_samples"])")

            gt_num_array,col_selectedcolumns = select_columns(parsed_args["select_samples"], gt_num_array, sample_names)
            sample_names=col_selectedcolumns

        elseif parsed_args["select_samples"] != nothing && length(parsed_args["group_samples"]) == 2
            println("Can not select samples with phenotype matrix provided. Please include same sample ids in phenotype matrix as in list of sample names to select.")

        end

        if parsed_args["num_array"] == true
        save_numerical_array(gt_num_array,sample_names,chr_pos_tuple_list,title,output_directory)
        end

        graphic = VariantVisualization.genotype_heatmap2(gt_num_array,title,chrom_label_info,sample_names,chr_pos_tuple_list,y_axis_label_option,x_axis_label_option)

    end

    println("Saving genotype heatmap")

    save_graphic(graphic,output_directory,save_ext,title,remote_option)

end

if parsed_args["heatmap"] == "read_depth"

    dp_num_array,dp_chromosome_labels = combined_all_read_depth_array_functions(sub)

    if parsed_args["heatmap_title"] != nothing
        title = "Read_Depth_$(parsed_args["heatmap_title"])"
    else
        bn = Base.Filesystem.basename(parsed_args["vcf_file"])
        title = "Read_Depth_$bn"
    end

    chr_pos_tuple_list = generate_chromosome_positions_for_hover_labels(dp_chromosome_labels)

    chrom_label_info = VariantVisualization.chromosome_label_generator(dp_chromosome_labels[:,1])

    if length(parsed_args["group_samples"]) == 2

        if parsed_args["select_samples"] != nothing && length(parsed_args["group_samples"]) != 2
            println("Selecting samples listed in $(parsed_args["select_samples"])")
            dp_num_array,col_selectedcolumns = select_columns(parsed_args["select_samples"], dp_num_array, sample_names)
            sample_names=col_selectedcolumns


        elseif parsed_args["select_samples"] != nothing && length(parsed_args["group_samples"]) == 2
            println("Can not select samples with phenotype matrix provided. Please include same sample ids in phenotype matrix as in list of sample names to select.")
        end

        group_trait_matrix_filename=((parsed_args["group_samples"])[1])
        trait_to_group_by = ((parsed_args["group_samples"])[2])
        println()
        println("Grouping samples by $trait_to_group_by")
        println()

        ordered_dp_num_array,group_label_pack,pheno,id_list,trait_labels = sortcols_by_phenotype_matrix(group_trait_matrix_filename, trait_to_group_by, dp_num_array, sample_names)
        dp_num_array_limited=read_depth_threshhold(ordered_dp_num_array)

        if parsed_args["num_array"] == true
        save_numerical_array(ordered_dp_num_array,sample_names,chr_pos_tuple_list,title,output_directory)
        end

        pheno_num_array,trait_label_array,chrom_label_info = add_pheno_matrix_to_dp_data_for_plotting(dp_num_array_limited,pheno,trait_labels,chrom_label_info,number_rows)

        graphic = VariantVisualization.dp_heatmap2_with_groups(pheno_num_array,title,chrom_label_info,group_label_pack,id_list,chr_pos_tuple_list,y_axis_label_option,trait_label_array,x_axis_label_option,number_rows)

    else

        if parsed_args["select_samples"] != nothing && length(parsed_args["group_samples"]) != 2
            println("Selecting samples listed in $(parsed_args["select_samples"])")

            dp_num_array,col_selectedcolumns = select_columns(parsed_args["select_samples"], dp_num_array, sample_names)
            sample_names=col_selectedcolumns


        elseif parsed_args["select_samples"] != nothing && length(parsed_args["group_samples"]) == 2
            println("Can not select samples with phenotype matrix provided. Please include same sample ids in phenotype matrix as in list of sample names to select.")

        end

        if parsed_args["num_array"] == true
        save_numerical_array(dp_num_array,sample_names,chr_pos_tuple_list,title,output_directory)
        end

        dp_num_array_limited=read_depth_threshhold(dp_num_array)
        graphic = VariantVisualization.dp_heatmap2(dp_num_array, title, chrom_label_info, sample_names,chr_pos_tuple_list,y_axis_label_option,x_axis_label_option)
    end

    println("Saving read depth heatmap")

    save_graphic(graphic,output_directory,save_ext,title,remote_option)
end

if parsed_args["heatmap"] == "genotype,read_depth" || parsed_args["heatmap"] == "read_depth,genotype"

    gt_num_array,gt_chromosome_labels = combined_all_genotype_array_functions(sub)

    if parsed_args["heatmap_title"] != nothing
        title = "Genotype_$(parsed_args["heatmap_title"])"
    else
        bn = Base.Filesystem.basename(parsed_args["vcf_file"])
        title = "Genotype_$bn"
    end

    chr_pos_tuple_list = generate_chromosome_positions_for_hover_labels(gt_chromosome_labels)

    chrom_label_info = VariantVisualization.chromosome_label_generator(gt_chromosome_labels[:,1])

    if length(parsed_args["group_samples"]) == 2

        if parsed_args["select_samples"] != nothing && length(parsed_args["group_samples"]) != 2
            println("Selecting samples listed in $(parsed_args["select_samples"])")
            gt_num_array,col_selectedcolumns = select_columns(parsed_args["select_samples"],
                                          gt_num_array,
                                          sample_names)
            sample_names=col_selectedcolumns

        elseif parsed_args["select_samples"] != nothing && length(parsed_args["group_samples"]) == 2
            println("Can not select samples with phenotype matrix provided. Please include same sample ids in phenotype matrix as in list of sample names to select.")
        end

        group_trait_matrix_filename=((parsed_args["group_samples"])[1])
        trait_to_group_by = ((parsed_args["group_samples"])[2])
        println()
        println("Grouping samples by $trait_to_group_by")
        println()

        ordered_num_array,group_label_pack,pheno,id_list,trait_labels = sortcols_by_phenotype_matrix(group_trait_matrix_filename, trait_to_group_by, gt_num_array, sample_names)

        if parsed_args["num_array"] == true
        save_numerical_array(ordered_num_array,sample_names,chr_pos_tuple_list,title,output_directory)
        end

        pheno_num_array,trait_label_array,chrom_label_info=add_pheno_matrix_to_gt_data_for_plotting(ordered_num_array,pheno,trait_labels,chrom_label_info,number_rows)

        graphic = VariantVisualization.genotype_heatmap_with_groups(pheno_num_array,title,chrom_label_info,group_label_pack,id_list,chr_pos_tuple_list,y_axis_label_option,trait_label_array,x_axis_label_option,number_rows)

    else

        if parsed_args["select_samples"] != nothing && length(parsed_args["group_samples"]) != 2
            println("Selecting samples listed in $(parsed_args["select_samples"])")

            gt_num_array,col_selectedcolumns = select_columns(parsed_args["select_samples"], gt_num_array, sample_names)
            sample_names=col_selectedcolumns

        elseif parsed_args["select_samples"] != nothing && length(parsed_args["group_samples"]) == 2
            println("Can not select samples with phenotype matrix provided. Please include same sample ids in phenotype matrix as in list of sample names to select.")

        end
                graphic = VariantVisualization.genotype_heatmap2(gt_num_array,title,chrom_label_info,sample_names,chr_pos_tuple_list,y_axis_label_option,x_axis_label_option)
    end

    println("Saving genotype heatmap")

    if parsed_args["num_array"] == true
    save_numerical_array(gt_num_array,sample_names,chr_pos_tuple_list,title,output_directory)
    end

    save_graphic(graphic,output_directory,save_ext,title,remote_option)

    if parsed_args["heatmap_title"] != nothing
        title = "Read_Depth_$(parsed_args["heatmap_title"])"
    else
        bn = Base.Filesystem.basename(parsed_args["vcf_file"])
        title = "Read_Depth_$bn"
    end

    dp_num_array,dp_chromosome_labels = combined_all_read_depth_array_functions(sub)

    chr_pos_tuple_list = generate_chromosome_positions_for_hover_labels(dp_chromosome_labels)

    chrom_label_info = VariantVisualization.chromosome_label_generator(dp_chromosome_labels[:,1])

    if length(parsed_args["group_samples"]) == 2

        if parsed_args["select_samples"] != nothing && length(parsed_args["group_samples"]) != 2
            println("Selecting samples listed in $(parsed_args["select_samples"])")
            dp_num_array,col_selectedcolumns = select_columns(parsed_args["select_samples"], dp_num_array, sample_names)
            sample_names=col_selectedcolumns

        elseif parsed_args["select_samples"] != nothing && length(parsed_args["group_samples"]) == 2
        end

        group_trait_matrix_filename=((parsed_args["group_samples"])[1])
        trait_to_group_by = ((parsed_args["group_samples"])[2])

        ordered_dp_num_array,group_label_pack,pheno,id_list,trait_labels = sortcols_by_phenotype_matrix(group_trait_matrix_filename, trait_to_group_by, dp_num_array, sample_names)
        dp_num_array_limited=read_depth_threshhold(ordered_dp_num_array)

        if parsed_args["num_array"] == true
        save_numerical_array(ordered_dp_num_array,sample_names,chr_pos_tuple_list,title,output_directory)
        end

        pheno_num_array,trait_label_array,chrom_label_info = add_pheno_matrix_to_dp_data_for_plotting(dp_num_array_limited,pheno,trait_labels,chrom_label_info,number_rows)

        graphic = VariantVisualization.dp_heatmap2_with_groups(pheno_num_array,title,chrom_label_info,group_label_pack,id_list,chr_pos_tuple_list,y_axis_label_option,trait_label_array,x_axis_label_option,number_rows)

    else

        if parsed_args["select_samples"] != nothing && length(parsed_args["group_samples"]) != 2
            println("Selecting samples listed in $(parsed_args["select_samples"])")

            dp_num_array,col_selectedcolumns = select_columns(parsed_args["select_samples"], dp_num_array, sample_names)
            sample_names=col_selectedcolumns

        elseif parsed_args["select_samples"] != nothing && length(parsed_args["group_samples"]) == 2
        end

        if parsed_args["num_array"] == true
        save_numerical_array(dp_num_array,sample_names,chr_pos_tuple_list,title,output_directory)
        end

        dp_num_array_limited=read_depth_threshhold(dp_num_array)
        graphic = VariantVisualization.dp_heatmap2(dp_num_array, title, chrom_label_info, sample_names,chr_pos_tuple_list,y_axis_label_option,x_axis_label_option)

    end

    println("Saving read depth heatmap")

    save_graphic(graphic,output_directory,save_ext,title,remote_option)

end

if parsed_args["avg_dp"] == "sample"

    dp_num_array,dp_chromosome_labels=combined_all_read_depth_array_functions_for_avg_dp(sub)

    chr_pos_tuple_list = generate_chromosome_positions_for_hover_labels(dp_chromosome_labels)

    dp_num_array_limited=read_depth_threshhold(dp_num_array)
    avg_list = VariantVisualization.avg_dp_samples(dp_num_array)

    list = VariantVisualization.list_sample_names_low_dp(avg_list, sample_names)
    #DelimitedFiles.writedlm(joinpath("$(parsed_args["output_directory"])","Samples_with_low_dp.csv"),list, ",") #turn on when create argument for this and add argument for cutoff e.g. list samples with avg dp <30.
    #println("The following samples have read depth of under 15: $list")
    graphic = avg_sample_dp_scatter(avg_list,sample_names,x_axis_label_option)

    bn = Base.Filesystem.basename(parsed_args["vcf_file"])
    title = "Average_Sample_Read_Depth_$bn"

    save_graphic(graphic,output_directory,save_ext,title,remote_option)

elseif parsed_args["avg_dp"] == "variant"

    dp_num_array,dp_chromosome_labels=combined_all_read_depth_array_functions_for_avg_dp(sub)

    chr_pos_tuple_list = generate_chromosome_positions_for_hover_labels(dp_chromosome_labels)
    chrom_label_info = VariantVisualization.chromosome_label_generator(dp_chromosome_labels[:,1])

    dp_num_array_limited=read_depth_threshhold(dp_num_array)
    avg_list = VariantVisualization.avg_dp_variant(dp_num_array)

    list = VariantVisualization.list_variant_positions_low_dp(avg_list, dp_chromosome_labels)
    #DelimitedFiles.writedlm(joinpath("$(parsed_args["output_directory"])","Variant_positions_with_low_dp.csv"),list,",") #turn on when create argument for this and add argument for cutoff e.g. list samples with avg dp <30.
    #println("The following variants have read depth of less than 15: $list") #turn on when create argument for this and add argument for cutoff e.g. list samples with avg dp <30.

    bn = Base.Filesystem.basename(parsed_args["vcf_file"])
    title = "Average_Variant_Read_Depth$bn"

    graphic = avg_variant_dp_line_chart(avg_list,chr_pos_tuple_list,y_axis_label_option,chrom_label_info)
    save_graphic(graphic,output_directory,save_ext,title,remote_option)

elseif parsed_args["avg_dp"] == "variant,sample" || parsed_args["avg_dp"] == "sample,variant"

    dp_num_array,dp_chromosome_labels=combined_all_read_depth_array_functions_for_avg_dp(sub)

    chr_pos_tuple_list = generate_chromosome_positions_for_hover_labels(dp_chromosome_labels)
    chrom_label_info = VariantVisualization.chromosome_label_generator(dp_chromosome_labels[:,1])

    dp_num_array_limited=read_depth_threshhold(dp_num_array)
    avg_list = VariantVisualization.avg_dp_variant(dp_num_array)

    list = VariantVisualization.list_variant_positions_low_dp(avg_list, dp_chromosome_labels)
   #DelimitedFiles.writedlm(joinpath("$(parsed_args["output_directory"])","Variant_positions_with_low_dp.csv"),list,",") #turn on when create argument for this and add argument for cutoff e.g. list samples with avg dp <30.
    #println("The following variants have read depth of less than 15: $list")

    bn = Base.Filesystem.basename(parsed_args["vcf_file"])
    title = "Average_Variant_Read_Depth$bn"

    graphic = avg_variant_dp_line_chart(avg_list,chr_pos_tuple_list,y_axis_label_option,chrom_label_info)
    save_graphic(graphic,output_directory,save_ext,title,remote_option)

    avg_list = VariantVisualization.avg_dp_samples(dp_num_array)
    list = VariantVisualization.list_sample_names_low_dp(avg_list, sample_names)
    #DelimitedFiles.writedlm(joinpath("$(parsed_args["output_directory"])","Samples_with_low_dp.csv"),list, ",") #turn on when create argument for this and add argument for cutoff e.g. list samples with avg dp <30.
    #println("The following samples have read depth of under 15: $list")

    bn = Base.Filesystem.basename(parsed_args["vcf_file"])
    title = "Average_Sample_Read_Depth_$bn"

    graphic = avg_sample_dp_scatter(avg_list,sample_names,x_axis_label_option)

    save_graphic(graphic,output_directory,save_ext,title,remote_option)

elseif parsed_args["avg_dp"] != nothing
    println("--avg_dp flag did not find expected arguments: sample, variant, or sample,variant. Average read depth plot not saved.")
end

close(reader)
