I am wondering whether Gviz can plot thin box(noncoding part) directly from GRange
Entering edit mode
Last seen 4 weeks ago

Hi Dr.Ivanek.I want to plot a Gene structure using GeneRegionTrack and I use the script and test gtf below.I can not plot the noncoding part(5UTR,3UTR).But I make it using the Txdb.In the Gviz Vignettes, I find the below sentence:

A nice bonus when building GeneRegionTracks from TxDb objects is that we get additional information about coding and non-coding regions of the transcripts, i.e., coordinates of the 5’ and 3’ UTRs and of the CDS regions. The class’ plotting method can use this information to distinguish between coding and non-coding regions based on the shape of the feature. All coding regions are plotted just as we have seen in the previous examples, whereas the non-coding regions are drawn as slightly thinner boxes. The distinction between coding and non-coding is made on the basis of the object’s feature values in combination with a special display parameter thinBoxFeature that enumerates all feature types that are to be treated as non-coding. Obviously this feature is available to all GeneRegionTracks, not only the ones that were build from TxDb objects. However, the coding information has to be added manually and the default value of the thinBoxFeature parameter may not be sufficient to cover all possible cases. It is up to the user to come up with a complete list of non-coding feature types depending on the source of the data.

But I did not know how to make it using the GeneRegionTracks made by GTF.

Here is the code

# loda package

# load gene annotation
At_gtf <- import.gff("test.gtf") 

At_gtf_new <- GRanges(At_gtf,
                      source = At_gtf$source,
                      feature = At_gtf$type,
                      transcript = At_gtf$transcript_id,
                      gene = At_gtf$gene_id)

# set gene name
gene_name <- "AT5G59340"

# get gene Grange location
gene_region_select <- subset(At_gtf_new,gene == gene_name)
gene_region_select$feature <- as.character(gene_region_select$feature)

# get gene start,end
region_start <- start(gene_region_select) %>% min() - 500
region_end <- end(gene_region_select) %>% max() + 500
region_Chr <- seqnames(gene_region_select)[1] %>% as.character()

region_region_select_total_GR <- GRanges(seqnames = region_Chr,
                                         ranges = IRanges(start = region_start,
                                                          end = region_end))

genemodel_track <- GeneRegionTrack(gene_region_select,
                                   collapseTranscripts = "longest",
                                   name = "gene_model",
                                   transcriptAnnotation = "gene")

Chr5    Araport11   3UTR    23933338    23933407    .   -   .   transcript_id "AT5G59340.1"; gene_id "AT5G59340";
Chr5    Araport11   exon    23933338    23933889    .   -   .   transcript_id "AT5G59340.1"; gene_id "AT5G59340";
Chr5    Araport11   CDS 23933408    23933889    .   -   2   transcript_id "AT5G59340.1"; gene_id "AT5G59340";
Chr5    Araport11   start_codon 23933887    23933889    .   -   .   transcript_id "AT5G59340.1"; gene_id "AT5G59340";
Chr5    Araport11   stop_codon  23934327    23934329    .   -   .   transcript_id "AT5G59340.1"; gene_id "AT5G59340";
Chr5    Araport11   CDS 23934327    23934627    .   -   0   transcript_id "AT5G59340.1"; gene_id "AT5G59340";
Chr5    Araport11   exon    23934327    23935050    .   -   .   transcript_id "AT5G59340.1"; gene_id "AT5G59340";
Chr5    Araport11   5UTR    23934628    23935050    .   -   .   transcript_id "AT5G59340.1"; gene_id "AT5G59340";
Gviz • 273 views
Entering edit mode
Robert Ivanek ▴ 700
Last seen 4 days ago

Dear Guandong Shang,

You can create TxDb object form your GTF file by using makeTxDbFromGFF function from GenomicFeatures package. Afterwards you can make GeneRegionTrack and plot it with Gviz.

Best Robert

TxDb <- makeTxDbFromGFF("test.gtf")
genemodel_track <- GeneRegionTrack(TxDb, collapseTranscripts = "longest",
                                  name = "gene_model", transcriptAnnotation = "gene")

Login before adding your answer.

Traffic: 378 users visited in the last hour
Help About
Access RSS

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6