From 1c127fe1dd099afeb714d5618b6b9b603fb5a015 Mon Sep 17 00:00:00 2001 From: Pavel Dimens Date: Mon, 8 Feb 2021 22:15:58 -0600 Subject: [PATCH 1/2] tweak VariantCall to be available in namespace --- src/PopGen.jl | 5 +- src/io/VariantCall.jl | 135 +---------------- ...{VariantCallGz.jl => VariantCallGzLazy.jl} | 0 src/io/VariantCallLazy.jl | 142 ++++++++++++++++++ 4 files changed, 147 insertions(+), 135 deletions(-) rename src/io/{VariantCallGz.jl => VariantCallGzLazy.jl} (100%) create mode 100755 src/io/VariantCallLazy.jl diff --git a/src/PopGen.jl b/src/PopGen.jl index 99f3cf186..fc0db7c79 100755 --- a/src/PopGen.jl +++ b/src/PopGen.jl @@ -56,11 +56,12 @@ include("io/Delimited.jl") include("io/Genepop.jl") include("io/Read.jl") include("io/Structure.jl") +include("io/VariantCall.jl") @init @require GeneticVariation="9bc6ac9d-e6b2-5f70-b0a8-242a01662520" begin - include("io/VariantCall.jl") + include("io/VariantCallLazy.jl") end @init @require GeneticVariation="9bc6ac9d-e6b2-5f70-b0a8-242a01662520" begin - @require GZip="92fee26a-97fe-5a0c-ad85-20a5f3185b63" include("io/VariantCallGz.jl") + @require GZip="92fee26a-97fe-5a0c-ad85-20a5f3185b63" include("io/VariantCallGzLazy.jl") end # example data include("io/Datasets.jl") diff --git a/src/io/VariantCall.jl b/src/io/VariantCall.jl index 393f424a7..efa17b979 100755 --- a/src/io/VariantCall.jl +++ b/src/io/VariantCall.jl @@ -1,22 +1,5 @@ -using .GeneticVariation - export bcf, vcf -""" - openvcf(::String) -Open VCF file (`.vcf/.gz`, or `.bcf/.gz`) and return an `IO` stream in reading mode `"r"`. -""" -function openvcf(infile::String) - if endswith(infile, ".vcf") || endswith(infile, ".bcf") - return open(infile, "r") - elseif endswith(infile, ".gz") - throw(ArgumentError("Please load in GZip.jl with \`using GZip\`")) - else - throw(ArgumentError("The filename must end with .vcf or .bcf")) - end -end - - """ bcf(infile::String; ; rename_snp::Bool, silent::Bool, allow_monomorphic::Bool) Load a BCF file into memory as a PopData object. Population information needs to be provided separately. @@ -44,64 +27,7 @@ julia> mydata.loci.genotype = mydata.loci.genotype |> Array{Union{Missing, NTup ``` """ function bcf(infile::String; rename_loci::Bool = false, silent::Bool = false, allow_monomorphic::Bool = false) - bases = (A = Int8(1), T = Int8(2), C = Int8(3), G = Int8(4), miss = Int8(0)) - stream = BCF.Reader(openvcf(infile)) - nmarkers = countlines(openvcf(infile)) - length(BCF.header(stream)) - 1 - sample_ID = header(stream).sampleID - nsamples = length(sample_ID) - loci_names = fill("marker", nmarkers) - geno_df = DataFrame(:name => sample_ID, :population => "missing") - if silent == false - @info "\n$(abspath(infile))\n$nsamples samples detected\n$nmarkers markers detected\npopulation info must be added <---" - end - for record in stream - ref_alt = Dict(-1 => "miss", 0 => BCF.ref(record), [i => j for (i,j) in enumerate(BCF.alt(record))]...) - raw_geno = BCF.genotype(record, 1:nsamples, "GT") - conv_geno = map(raw_geno) do rg - tmp = replace.(rg, "." => "-1") - ig = collect(parse.(Int8, split(tmp, r"\/|\|"))) - [bases[Symbol(ref_alt[i])] for i in ig] |> sort |> Tuple - end - insertcols!(geno_df, Symbol(BCF.chrom(record) * "_" * string(BCF.pos(record))) => conv_geno) - end - close(stream) - if rename_snp - rnm = append!([:name, :population], [Symbol.("snp_" * i) for i in string.(1:nmarkers)]) - rename!(geno_df, rnm) - end - stacked_geno_df = DataFrames.stack(geno_df, DataFrames.Not(1:2)) - rename!(stacked_geno_df, [:name, :population, :locus, :genotype]) - # set columns as PooledArrays - select!( - stacked_geno_df, - :name => PooledArray => :name, - :population => PooledArray => :population, - :locus => (i -> PooledArray(i |> Vector{String})) => :locus, - :genotype - ) - # replace missing genotypes as missing - stacked_geno_df.genotype = map(stacked_geno_df.genotype) do geno - if all(0 .== geno) - return missing - else - return geno - end - end - sort!(stacked_geno_df, [:name, :locus]) - #meta_df = generate_meta(stacked_geno_df) - # ploidy finding - meta_df = DataFrames.combine(DataFrames.groupby(stacked_geno_df, :name), - :genotype => (i -> Int8(length(first(skipmissing(i))))) => :ploidy - ) - insertcols!(meta_df, 2, :population => "missing") - insertcols!(meta_df, 4, :longitude => Vector{Union{Missing, Float32}}(undef, nsamples), :latitude => Vector{Union{Missing, Float32}}(undef, nsamples)) - - if allow_monomorphic - pd_out = PopData(meta_df, stacked_geno_df) - else - pd_out = drop_monomorphic!(PopData(meta_df, stacked_geno_df)) - end - return pd_out + error("Please load in GeneticVariation.jl with \`using GeneticVariation\`") end ### VCF parsing ### @@ -134,62 +60,5 @@ julia> mydata.loci.genotype = mydata.loci.genotype |> Array{Union{Missing, NTup ``` """ function vcf(infile::String; rename_snp::Bool = false, silent::Bool = false, allow_monomorphic::Bool = false) - bases = (A = Int8(1), T = Int8(2), C = Int8(3), G = Int8(4), miss = Int8(0)) - stream = VCF.Reader(openvcf(infile)) - nmarkers = countlines(openvcf(infile)) - length(VCF.header(stream)) - 1 - sample_ID = header(stream).sampleID - nsamples = length(sample_ID) - loci_names = fill("marker", nmarkers) - geno_df = DataFrame(:name => sample_ID, :population => "missing") - if silent == false - @info "\n$(abspath(infile))\n$nsamples samples detected\n$nmarkers markers detected\npopulation info must be added <---" - end - for record in stream - ref_alt = Dict(-1 => "miss", 0 => VCF.ref(record), [i => j for (i,j) in enumerate(VCF.alt(record))]...) - raw_geno = VCF.genotype(record, 1:nsamples, "GT") - conv_geno = map(raw_geno) do rg - tmp = replace.(rg, "." => "-1") - ig = collect(parse.(Int8, split(tmp, r"\/|\|"))) - [bases[Symbol(ref_alt[i])] for i in ig] |> sort |> Tuple - end - insertcols!(geno_df, Symbol(VCF.chrom(record) * "_" * string(VCF.pos(record))) => conv_geno) - end - close(stream) - if rename_snp - rnm = append!([:name, :population], [Symbol.("snp_" * i) for i in string.(1:nmarkers)]) - rename!(geno_df, rnm) - end - stacked_geno_df = DataFrames.stack(geno_df, DataFrames.Not(1:2)) - rename!(stacked_geno_df, [:name, :population, :locus, :genotype]) - # set columns as PooledArrays - select!( - stacked_geno_df, - :name => PooledArray => :name, - :population => PooledArray => :population, - :locus => (i -> PooledArray(i |> Vector{String})) => :locus, - :genotype - ) - # replace missing genotypes as missing - stacked_geno_df.genotype = map(stacked_geno_df.genotype) do geno - if all(0 .== geno) - return missing - else - return geno - end - end - sort!(stacked_geno_df, [:name, :locus]) - # ploidy finding - #meta_df = generate_meta(stacked_geno_df) - meta_df = DataFrames.combine(DataFrames.groupby(stacked_geno_df, :name), - :genotype => (i -> Int8(length(first(skipmissing(i))))) => :ploidy - ) - insertcols!(meta_df, 2, :population => "missing") - insertcols!(meta_df, 4, :longitude => Vector{Union{Missing, Float32}}(undef, nsamples), :latitude => Vector{Union{Missing, Float32}}(undef, nsamples)) - - if allow_monomorphic - pd_out = PopData(meta_df, stacked_geno_df) - else - pd_out = drop_monomorphic!(PopData(meta_df, stacked_geno_df)) - end - return pd_out + error("Please load in GeneticVariation.jl with \`using GeneticVariation\`") end diff --git a/src/io/VariantCallGz.jl b/src/io/VariantCallGzLazy.jl similarity index 100% rename from src/io/VariantCallGz.jl rename to src/io/VariantCallGzLazy.jl diff --git a/src/io/VariantCallLazy.jl b/src/io/VariantCallLazy.jl new file mode 100755 index 000000000..0d52871ee --- /dev/null +++ b/src/io/VariantCallLazy.jl @@ -0,0 +1,142 @@ +using .GeneticVariation + +export bcf, vcf + +""" + openvcf(::String) +Open VCF file (`.vcf/.gz`, or `.bcf/.gz`) and return an `IO` stream in reading mode `"r"`. +""" +function openvcf(infile::String) + if endswith(infile, ".vcf") || endswith(infile, ".bcf") + return open(infile, "r") + elseif endswith(infile, ".gz") + throw(ArgumentError("Please load in GZip.jl with \`using GZip\`")) + else + throw(ArgumentError("The filename must end with .vcf or .bcf")) + end +end + + +function bcf(infile::String; rename_loci::Bool = false, silent::Bool = false, allow_monomorphic::Bool = false) + bases = (A = Int8(1), T = Int8(2), C = Int8(3), G = Int8(4), miss = Int8(0)) + stream = BCF.Reader(openvcf(infile)) + nmarkers = countlines(openvcf(infile)) - length(BCF.header(stream)) - 1 + sample_ID = header(stream).sampleID + nsamples = length(sample_ID) + loci_names = fill("marker", nmarkers) + geno_df = DataFrame(:name => sample_ID, :population => "missing") + if silent == false + @info "\n$(abspath(infile))\n$nsamples samples detected\n$nmarkers markers detected\npopulation info must be added <---" + end + for record in stream + ref_alt = Dict(-1 => "miss", 0 => BCF.ref(record), [i => j for (i,j) in enumerate(BCF.alt(record))]...) + raw_geno = BCF.genotype(record, 1:nsamples, "GT") + conv_geno = map(raw_geno) do rg + tmp = replace.(rg, "." => "-1") + ig = collect(parse.(Int8, split(tmp, r"\/|\|"))) + [bases[Symbol(ref_alt[i])] for i in ig] |> sort |> Tuple + end + insertcols!(geno_df, Symbol(BCF.chrom(record) * "_" * string(BCF.pos(record))) => conv_geno) + end + close(stream) + if rename_snp + rnm = append!([:name, :population], [Symbol.("snp_" * i) for i in string.(1:nmarkers)]) + rename!(geno_df, rnm) + end + stacked_geno_df = DataFrames.stack(geno_df, DataFrames.Not(1:2)) + rename!(stacked_geno_df, [:name, :population, :locus, :genotype]) + # set columns as PooledArrays + select!( + stacked_geno_df, + :name => PooledArray => :name, + :population => PooledArray => :population, + :locus => (i -> PooledArray(i |> Vector{String})) => :locus, + :genotype + ) + # replace missing genotypes as missing + stacked_geno_df.genotype = map(stacked_geno_df.genotype) do geno + if all(0 .== geno) + return missing + else + return geno + end + end + sort!(stacked_geno_df, [:name, :locus]) + #meta_df = generate_meta(stacked_geno_df) + # ploidy finding + meta_df = DataFrames.combine(DataFrames.groupby(stacked_geno_df, :name), + :genotype => (i -> Int8(length(first(skipmissing(i))))) => :ploidy + ) + insertcols!(meta_df, 2, :population => "missing") + insertcols!(meta_df, 4, :longitude => Vector{Union{Missing, Float32}}(undef, nsamples), :latitude => Vector{Union{Missing, Float32}}(undef, nsamples)) + + if allow_monomorphic + pd_out = PopData(meta_df, stacked_geno_df) + else + pd_out = drop_monomorphic!(PopData(meta_df, stacked_geno_df)) + end + return pd_out +end + +### VCF parsing ### + +function vcf(infile::String; rename_snp::Bool = false, silent::Bool = false, allow_monomorphic::Bool = false) + bases = (A = Int8(1), T = Int8(2), C = Int8(3), G = Int8(4), miss = Int8(0)) + stream = VCF.Reader(openvcf(infile)) + nmarkers = countlines(openvcf(infile)) - length(VCF.header(stream)) - 1 + sample_ID = header(stream).sampleID + nsamples = length(sample_ID) + loci_names = fill("marker", nmarkers) + geno_df = DataFrame(:name => sample_ID, :population => "missing") + if silent == false + @info "\n$(abspath(infile))\n$nsamples samples detected\n$nmarkers markers detected\npopulation info must be added <---" + end + for record in stream + ref_alt = Dict(-1 => "miss", 0 => VCF.ref(record), [i => j for (i,j) in enumerate(VCF.alt(record))]...) + raw_geno = VCF.genotype(record, 1:nsamples, "GT") + conv_geno = map(raw_geno) do rg + tmp = replace.(rg, "." => "-1") + ig = collect(parse.(Int8, split(tmp, r"\/|\|"))) + [bases[Symbol(ref_alt[i])] for i in ig] |> sort |> Tuple + end + insertcols!(geno_df, Symbol(VCF.chrom(record) * "_" * string(VCF.pos(record))) => conv_geno) + end + close(stream) + if rename_snp + rnm = append!([:name, :population], [Symbol.("snp_" * i) for i in string.(1:nmarkers)]) + rename!(geno_df, rnm) + end + stacked_geno_df = DataFrames.stack(geno_df, DataFrames.Not(1:2)) + rename!(stacked_geno_df, [:name, :population, :locus, :genotype]) + # set columns as PooledArrays + select!( + stacked_geno_df, + :name => PooledArray => :name, + :population => PooledArray => :population, + :locus => (i -> PooledArray(i |> Vector{String})) => :locus, + :genotype + ) + # replace missing genotypes as missing + stacked_geno_df.genotype = map(stacked_geno_df.genotype) do geno + if all(0 .== geno) + return missing + else + return geno + end + end + sort!(stacked_geno_df, [:name, :locus]) + # ploidy finding + #meta_df = generate_meta(stacked_geno_df) + meta_df = DataFrames.combine(DataFrames.groupby(stacked_geno_df, :name), + :genotype => (i -> Int8(length(first(skipmissing(i))))) => :ploidy + ) + insertcols!(meta_df, 2, :population => "missing") + insertcols!(meta_df, 4, :longitude => Vector{Union{Missing, Float32}}(undef, nsamples), :latitude => Vector{Union{Missing, Float32}}(undef, nsamples)) + + if allow_monomorphic + pd_out = PopData(meta_df, stacked_geno_df) + else + pd_out = drop_monomorphic!(PopData(meta_df, stacked_geno_df)) + end + return pd_out +end From 15d5113e8b012393cea3d047f50e9d1b1337dbb1 Mon Sep 17 00:00:00 2001 From: Pavel Dimens Date: Mon, 8 Feb 2021 22:25:41 -0600 Subject: [PATCH 2/2] remove extension deps for openvcf() --- src/io/VariantCallGzLazy.jl | 8 +++----- src/io/VariantCallLazy.jl | 8 +++----- 2 files changed, 6 insertions(+), 10 deletions(-) diff --git a/src/io/VariantCallGzLazy.jl b/src/io/VariantCallGzLazy.jl index 161e5c3a6..4318003b0 100644 --- a/src/io/VariantCallGzLazy.jl +++ b/src/io/VariantCallGzLazy.jl @@ -2,11 +2,9 @@ using .GZip # if GZip is loaded in, overwrite openvcf to this method function openvcf(infile::String) - if endswith(infile, ".vcf") || endswith(infile, ".bcf") - return open(infile, "r") - elseif endswith(infile, ".vcf.gz") || endswith(infile, ".bcf.gz") + if endswith(infile, ".gz") return GZip.open(infile, "r") else - throw(ArgumentError("The filename must end with .vcf/.bcf or .vcf.gz/.bcf.gz")) + return open(infile, "r") end -end +end \ No newline at end of file diff --git a/src/io/VariantCallLazy.jl b/src/io/VariantCallLazy.jl index 0d52871ee..6baf5bcd4 100755 --- a/src/io/VariantCallLazy.jl +++ b/src/io/VariantCallLazy.jl @@ -4,15 +4,13 @@ export bcf, vcf """ openvcf(::String) -Open VCF file (`.vcf/.gz`, or `.bcf/.gz`) and return an `IO` stream in reading mode `"r"`. +Open VCF file (`.vcf(.gz)`, or `.bcf(.gz)`) and return an `IO` stream in reading mode `"r"`. """ function openvcf(infile::String) - if endswith(infile, ".vcf") || endswith(infile, ".bcf") - return open(infile, "r") - elseif endswith(infile, ".gz") + if endswith(infile, ".gz") throw(ArgumentError("Please load in GZip.jl with \`using GZip\`")) else - throw(ArgumentError("The filename must end with .vcf or .bcf")) + return open(infile, "r") end end