Entering edit mode
7.3 years ago
longoka
▴
40
I'm trying to create the plot below, but with only a single transcript (uc062uto.1) displayed, not all of the variants (see below).
What am I missing?
The code I'm using to make this image:
# simple genomic/transcript map
rm(list = ls()) # clear workspace
Gene <- "HTT"
require(biomaRt)
require(ggbio)
require(Biostrings)
require(ShortRead)
require(BSgenome)
require(GenomicFeatures)
## get gene ranges, strand, based on gene name (biomaRt functions)
mart <- useMart("ENSEMBL_MART_ENSEMBL", host = "useast.ensembl.org")
datasets <- listDatasets(mart)
mart <- useDataset("hsapiens_gene_ensembl",mart)
r1 <- getBM(filters="external_gene_name",
values = Gene,
attributes=c("external_gene_name",
"chromosome_name",
"start_position",
"end_position",
"strand",
"refseq_mrna",
"entrezgene",
"ucsc",
"ensembl_gene_id",
"hamap"),
mart = mart)
r1 <- r1[r1$refseq_mrna != "",]
CHR <- paste("chr",r1$chromosome_name[1], sep = "")
START <- r1$start_position[1]
END <- r1$end_position[1]
if (r1$strand[1] == 1) {
STRAND = "+"
} else {
STRAND = "-"
}
library(TxDb.Hsapiens.UCSC.hg38.knownGene) # knowngene db for hg38
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
gr <- GRanges(CHR, IRanges(START, END), strand = STRAND)
p <- autoplot(txdb,
which = gr,
ylab = paste(Gene," [",CHR,":",START,"-",END,"]", sep = ""),
geom = "full",
coord = "genome",
layout = "linear",
strand = STRAND)
p
Works like a charm; I think I need to verse myself a bit better on all the functions in the GenomicFeatures package. Thank you!