Entering edit mode
rop.jesse
•
0
@ropjesse-14657
Last seen 6.9 years ago
I am trying to use the trackViewer package to create a lolliplot plot to portray the variants in the fut3 gene from the 1000 genomes .popvcf file for this locus. I am able to replicate the example given in the trackViewer vignette but when I alter the code to plot fut3 variants I get this error, "Error: allwidthSNP.gr[[i]]) == 1) is not TRUE"
Below is the code from the trackViewer vignette that works fine.
> library(VariantAnnotation)
> library(TxDb.Hsapiens.UCSC.hg19.knownGene)
> library(org.Hs.eg.db)
> fl <- system.file("extdata", "chr22.vcf.gz", package="VariantAnnotation")
> gr <- GRanges("22", IRanges(50968014, 50970514, names="TYMP"))
> tab <- TabixFile(fl)
> vcf <- readVcf(fl, "hg19", param=gr)
> mmutation.frequency <- rowRanges(vcf)
> mcols(mmutation.frequency) <- cbind(mcols(mmutation.frequency),
+ VariantAnnotation::info(vcf))
> mmutation.frequency$border <- "gray30"
> mmutation.frequency$color <-
+ ifelse(grepl("^rs", names(mmutation.frequency)), "lightcyan", "lavender")
> ## plot Global Allele Frequency based on AC/AN
> mmutation.frequency$score <- mmutation.frequency$AF*100
> seqlevelsStyle(gr) <- seqlevelsStyle(mmutation.frequency) <- "UCSC"
> trs <- geneModelFromTxdb(TxDb.Hsapiens.UCSC.hg19.knownGene,
+ org.Hs.eg.db,
+ gr=gr)
'select()' returned 1:1 mapping between keys and columns
There were 15 warnings (use warnings() to see them)
> features <- c(range(trs[[1]]@dat), range(trs[[5]]@dat))
> names(features) <- c(trs[[1]]@name, trs[[5]]@name)
> features$fill <- c("lightblue", "mistyrose")
> features$height <- c(.02, .04)
> lolliplot(mmutation.frequency, features, ranges=gr)
>
Below is the code after I have made adjustments to plot the lolliplot for variants in the fut3 gene.
> fut3_gr = GRanges("19", IRanges(5842040,5852344, names="FUT3"))
> bgzip( "F:/WT 18 month project/Analysis/fut3.popvcf")
Error in .zip(.bgzip, file, dest, overwrite) : 'dest' exists:
dest: F:/WT 18 month project/Analysis/fut3.popvcf.bgz
> indexTabix ( "F:/WT 18 month project/Analysis/fut3.popvcf.bgz", format = "vcf4")
[1] "F:/WT 18 month project/Analysis/fut3.popvcf.bgz.tbi"
> fut3_path = "F:/WT 18 month project/Analysis/fut3.popvcf.bgz"
> fut3_tab = TabixFile(fut3_path)
> fut3_vcf = readVcf(fut3_path, "hg19", param = fut3_gr)
> mutation.frequency <- rowRanges(fut3_vcf)
> mcols(mutation.frequency) <- cbind(mcols(mutation.frequency), VariantAnnotation::info(fut3_vcf))
> mutation.frequency$border <- "gray30"
> mutation.frequency$color <- ifelse(grepl("^rs", names(mutation.frequency)), "lightcyan", "lavender")
> mutation.frequency$score <- mutation.frequency$AF*100
> seqlevelsStyle(fut3_gr) <- seqlevelsStyle(mutation.frequency) <- "UCSC"
> fut3_trs <- geneModelFromTxdb (TxDb.Hsapiens.UCSC.hg19.knownGene, org.Hs.eg.db, gr=fut3_gr)
'select()' returned 1:1 mapping between keys and columns
There were 15 warnings (use warnings() to see them)
> fut3_features <- c(range(fut3_trs[[1]]@dat))
> names(fut3_features) <- c(fut3_trs[[1]]@name)
> fut3_features$fill <- c("lightblue")
> fut3_features$height <- c(.02)
> lolliplot(mutation.frequency, fut3_features, ranges=fut3_gr)
Error: all(width(SNP.gr[[i]]) == 1) is not TRUE
Below is the output from the traceback() command
> traceback()
3: stop(sprintf(ngettext(length(r), "%s is not TRUE", "%s are not all TRUE"),
ch), call. = FALSE, domain = NA)
2: stopifnot(all(width(SNP.gr[[i]]) == 1))
1: lolliplot(mutation.frequency, fut3_features, ranges = fut3_gr)
Below is the output from sessionInfo() command
> sessionInfo() R version 3.3.3 (2017-03-06) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows >= 8 x64 (build 9200) locale: [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C [5] LC_TIME=English_United States.1252 attached base packages: [1] grid stats4 parallel stats graphics grDevices utils datasets methods [10] base other attached packages: [1] dplyr_0.7.4 org.Hs.eg.db_3.4.0 [3] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2 GenomicFeatures_1.26.4 [5] AnnotationDbi_1.36.2 trackViewer_1.10.2 [7] VariantAnnotation_1.20.3 Rsamtools_1.26.2 [9] Biostrings_2.42.1 XVector_0.14.1 [11] SummarizedExperiment_1.4.0 Biobase_2.34.0 [13] GenomicRanges_1.26.4 GenomeInfoDb_1.10.3 [15] IRanges_2.8.2 S4Vectors_0.12.2 [17] BiocGenerics_0.20.0 loaded via a namespace (and not attached): [1] bitops_1.0-6 matrixStats_0.52.2 bit64_0.9-7 [4] RColorBrewer_1.1-2 httr_1.3.1 tools_3.3.3 [7] backports_1.1.2 R6_2.2.2 rpart_4.1-11 [10] Hmisc_4.0-3 DBI_0.7 lazyeval_0.2.1 [13] Gviz_1.18.2 colorspace_1.3-2 nnet_7.3-12 [16] gridExtra_2.3 bit_1.1-12 htmlTable_1.11.0 [19] grImport_0.9-0 rtracklayer_1.34.2 scales_0.5.0 [22] checkmate_1.8.5 pbapply_1.3-3 stringr_1.2.0 [25] digest_0.6.13 foreign_0.8-69 base64enc_0.1-3 [28] dichromat_2.0-0 pkgconfig_2.0.1 htmltools_0.3.6 [31] ensembldb_1.6.2 BSgenome_1.42.0 htmlwidgets_0.9 [34] rlang_0.1.4 rstudioapi_0.7 RSQLite_2.0 [37] BiocInstaller_1.24.0 shiny_1.0.5 bindr_0.1 [40] BiocParallel_1.8.2 acepack_1.4.1 RCurl_1.95-4.8 [43] magrittr_1.5 Formula_1.2-2 Matrix_1.2-12 [46] Rcpp_0.12.14 munsell_0.4.3 stringi_1.1.6 [49] yaml_2.1.16 zlibbioc_1.20.0 plyr_1.8.4 [52] AnnotationHub_2.6.5 blob_1.1.0 lattice_0.20-35 [55] splines_3.3.3 knitr_1.17 biomaRt_2.30.0 [58] XML_3.98-1.9 glue_1.2.0 biovizBase_1.22.0 [61] latticeExtra_0.6-28 data.table_1.10.4-3 httpuv_1.3.5 [64] gtable_0.2.0 purrr_0.2.4 tidyr_0.7.2 [67] assertthat_0.2.0 ggplot2_2.2.1 mime_0.5 [70] xtable_1.8-2 survival_2.41-3 tibble_1.3.4 [73] GenomicAlignments_1.10.1 memoise_1.1.0 bindrcpp_0.2 [76] cluster_2.0.6 interactiveDisplayBase_1.12.0 >
Many thanks!!!
