Hi I am wondering which way could allow me to retrieve the first exons from genome with the information of the exon length. the purpose is to retrieve the first exons for gRNA design. Thanks, Xin
Hi,
How about extracting the exons from a TxDb object and pick up the "first exons"? (whatever that means):
library(TxDb.Dmelanogaster.UCSC.dm3.ensGene)
txdb <- TxDb.Dmelanogaster.UCSC.dm3.ensGene
ex <- exons(txdb)
ex
# GRanges object with 76920 ranges and 1 metadata column:
# seqnames ranges strand | exon_id
# <Rle> <IRanges> <Rle> | <integer>
# [1] chr2L 7529-8116 + | 1
# [2] chr2L 8193-8589 + | 2
# [3] chr2L 8193-9484 + | 3
# [4] chr2L 8229-9484 + | 4
# [5] chr2L 8668-9484 + | 5
# ... ... ... ... . ...
# [76916] chrYHet 327915-328204 - | 76916
# [76917] chrYHet 328270-328489 - | 76917
# [76918] chrUextra 523024-523048 - | 76918
# [76919] chrUextra 523024-523086 - | 76919
# [76920] chrUextra 523060-523086 - | 76920
# -------
# seqinfo: 15 sequences (1 circular) from dm3 genome
The exons returned by exons() are guaranteed to be sorted first by chromosome, then by strand, then by ascending position on the chromosome.
H.
Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
version 2.3.6
Hi Hervé, Thank you very much. Sorry for the unclear description. Actually, I want to retrieve the first exons for each gene where RNA polymerase begins transcribing. The reason that I will use the retrieved exon for CRISPR gRNA design as targeting on these regions would give high efficient knockout.
A gene has one TSS per transcript. You can easily get the 1st exon of each known transcript by extracting the exons grouped by transcript:
ex_by_txis a named GRangesList object with 1 GRanges object per known transcript. The names onex_by_txare the transcript internal ids (these ids are internal to the TxDb object and have no meaning outside it). Each GRanges object inex_by_txcontains the exons for the corresponding transcript. The exons are ordered by rank with respect to their transcript. See?exonsByfor more information.To keep the 1st exon for each transcript we can use
heads(), a version ofhead()that works on a list-like object and keeps the firstnelements within each list element (note thatheads()is a Bioconductor extension and is not part of base R):We can turn this into a GRanges object with
unlist():and add some useful metadata columns to it:
This final version of
first_exonscontains the genomic location of the 1st exon of each known transcript and each exon is mapped to the corresponding transcript name, gene name, and transcript length.Hope this helps,
H.
I should add that the fact that some genomic ranges are repeated multiple times in
first_exons(e.g. the first 3 exons are the same) simply reflects the fact that several transcripts in a gene can share the same first exon. So if you're interested in knowing the location of the first exon for each gene, you can start by getting rid of the repeated exons with:Even after doing this
first_exonscan still contain more than one "first exon" per gene (e.g. gene FBgn0067779 contains 4 "first exons"). To keep only one "first exon" per gene, I guess the natural thing to do is to go after "the most upstream first exon" for each gene, that is, the exon with the smaller start position for genes on the + strand and the exon with the greatest end position for genes on the minus strand. Here is one way to do this:H.
One more thing (and I'll stop replying to myself).
The very first exon in each gene can actually be obtained more simply with:
This method is actually doing the right thing if the goal is to extract the most upstream exon for each gene in the TxDb object. Unlike in
very_first_exon_in_gene, all the known genes (i.e. all the genes in the TxDb object) are represented invery_first_exon_in_each_gene.The previous method was loosing some genes along the way due to the fact that some genes can share the same most upstream exon and that we were calling
unique()at some point. For example genes FBgn0000055 and FBgn0000056 share the same most upstream exon:which is why FBgn0000056 didn't make it to
very_first_exon_in_gene:H.
Thank you very much, Hervé. Really helpful!
Hi Hervé, I am trying to design the genome-wide gRNA libraries for several different species. One paper described a pipeline for the design. Basically, they first select the exonic guide sites fitting the pattern G(N16–19)NGG, and then annotate these guide sites as targeting Ensembl GRCh37 genes models to generate candidate guides towards each gene. which R package would you recommend for this purpose. Thanks
I don't know. This sounds like is a slightly different topic though so I would recommend that maybe you try to ask this as a new question.