viva

#!/usr/bin/env julia

println()
println("Welcome to VIVA.")
println()
println("Loading dependency packages:")
println()

#println("1/5 ArgParse")
#import ArgParse

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

println("1/4 PlotlyJS")
using PlotlyJS

println("2/4 GeneticVariation")
using GeneticVariation

println("3/4 DelimitedFiles")
using DelimitedFiles

println("4/4 VariantVisualization")

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 = VariantVisualization.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 = GeneticVariation.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"]

#=

#improvement note for developers: implement feature to display VCF summary stats when no plotting options are chosen. This was implemented with VCFTools.jl package which is not supported in Julia v1+. Would need to find new function, perhaps in GeneticVariation.jl
#nrecords() and nsamples() were from VCFTools.jl package which is no longer supported
#can also implement this as flag option to display number variants and samples before filtering

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

            #sample selection occurs in this case during sample grouping step.
            #when using sample grouping and selection options, sample metadata matric should have the same sample names as the selection list.
            #improvement note for developers: make function to check that sample ids in metadata matrix and sample selection list match. If not, print error and choose to go with metadata matrix sample ids.

        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
            #sample selection occurs in this case during sample grouping step.
            #when using sample grouping and selection options, sample metadata matric should have the same sample names as the selection list.
            #improvement note for developers: make function to check that sample ids in metadata matrix and sample selection list match. If not, print error and choose to go with metadata matrix sample ids.
        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_new_legend(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
            #sample selection occurs in this case during sample grouping step.
            #when using sample grouping and selection options, sample metadata matric should have the same sample names as the selection list.
            #improvement note for developers: make function to check that sample ids in metadata matrix and sample selection list match. If not, print error and choose to go with metadata matrix sample ids.
        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
            #sample selection occurs in this case during sample grouping step.
            #when using sample grouping and selection options, sample metadata matric should have the same sample names as the selection list.
            #improvement note for developers: make function to check that sample ids in metadata matrix and sample selection list match. If not, print error and choose to go with metadata matrix sample ids.

        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
            #sample selection occurs in this case during sample grouping step.
            #when using sample grouping and selection options, sample metadata matric should have the same sample names as the selection list.
            #improvement note for developers: make function to check that sample ids in metadata matrix and sample selection list match. If not, print error and choose to go with metadata matrix sample ids.
        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
            #sample selection occurs in this case during sample grouping step.
            #when using sample grouping and selection options, sample metadata matric should have the same sample names as the selection list.
            #improvement note for developers: make function to check that sample ids in metadata matrix and sample selection list match. If not, print error and choose to go with metadata matrix sample ids.

        end
                graphic = VariantVisualization.genotype_heatmap2_new_legend(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

            #sample selection occurs in this case during sample grouping step.
            #when using sample grouping and selection options, sample metadata matric should have the same sample names as the selection list.
            #improvement note for developers: make function to check that sample ids in metadata matrix and sample selection list match. If not, print error and choose to go with metadata matrix sample ids.

        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

            #sample selection occurs in this case during sample grouping step.
            #when using sample grouping and selection options, sample metadata matric should have the same sample names as the selection list.
            #improvement note for developers: make function to check that sample ids in metadata matrix and sample selection list match. If not, print error and choose to go with metadata matrix sample ids.
        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) #close the VCF.Reader object to remove it from memory - important for running analysis of same file again and clearing memory space when done