knausb / vcfR

Tools to work with variant call format files
248 stars 54 forks source link

Hide a panel from chromoqc or plot #210

Closed NowickiLab closed 1 year ago

NowickiLab commented 1 year ago

Hi Brian, Love your package here! Please advice how to hide a panel of my choice from chromoqc or plot(chrom) - as in, if I don't have MQ in vcf, it's silly to show empty panels.

Thank you for any advice (and apologies if this Q has been already answered)! -Marcin

knausb commented 1 year ago

Hi,

how about the following?

library(vcfR)
#> 
#>    *****       ***   vcfR   ***       *****
#>    This is vcfR 1.13.0.9999 
#>      browseVignettes('vcfR') # Documentation
#>      citation('vcfR') # Citation
#>    *****       *****      *****       *****
pkg <- "pinfsc50"
vcf_file <- system.file("extdata", "pinf_sc50.vcf.gz", package = pkg)
dna_file <- system.file("extdata", "pinf_sc50.fasta", package = pkg)
gff_file <- system.file("extdata", "pinf_sc50.gff", package = pkg)

vcf <- read.vcfR( vcf_file, verbose = FALSE )
dna <- ape::read.dna(dna_file, format = "fasta")
gff <- read.table(gff_file, sep="\t", quote="")

chrom <- create.chromR(name='Supercontig', vcf=vcf, seq=dna, ann=gff)
#> Names in vcf:
#>   Supercontig_1.50
#> Names of sequences:
#>   Supercontig_1.50 of Phytophthora infestans T30-4
#> Warning in create.chromR(name = "Supercontig", vcf = vcf, seq = dna, ann = gff): 
#>         Names in variant data and sequence data do not match perfectly.
#>         If you choose to proceed, we'll do our best to match the data.
#>         But prepare yourself for unexpected results.
#> Names in annotation:
#>   Supercontig_1.50
#> Initializing var.info slot.
#> var.info slot initialized.
chromoqc(chrom, dp.alpha=20)


# Set MQ to NULL.
chrom@var.info$MQ <- NULL
!is.null(chrom@var.info$MQ)
#> [1] FALSE
chromoqc(chrom, dp.alpha=20)


plot(chrom)


myDP <- extract.info(vcf, element = "DP", as.numeric = TRUE)
hist(myDP)

Created on 2023-02-03 with reprex v2.0.2

If we set the 'MQ' column to NULL it knows to omit it. The plot() function just pulls columns out of the INFO column. The function extract.info() should help you extract whatever is in your data so you can plot it however you want.

Does this provide you with a solution? Brian

NowickiLab commented 1 year ago

Perfection! Many thanks for your fast and detailed answer. Onwards with vcfR! -m