Parse GENCODE data to summarise some key features of the human genome
Go on the latest human release of GENCODE .
Download the Comprehensive gene annotation as Content and CHR as Regions. Of note, between GTF and GFF3 formats, pick the GFF3 one.
rtracklayer to parse GFF3 filesrtracklayer::import(). Assign to the name gencode.This steps requires some large memory, such as ~ 4 GB. So your computer needs at least 8 GB to perform this project. The object gencode below is 1.60 GB alone. Thus, please download in your data sub-folder the lighter file which contains only 14 chromosomes gencode_light.rds (21 MB). The compressed tsv file can be then loaded using the following command, this object is a tibble, no need to do the import step, go directly to the Explore section.
gencode <- read_tsv("data/gencode_light.tsv.gz")rtracklayer is a Bioconductor package. See its homepage for installation instructions.
import can read directly an URL.
import() returns an object of class GRanges. Convert to a tibble
Of note: This operation takes about 2-3 minutes. You may want to save the object while testing your workflow.
That is way too many lines for the number of genes we expect for the human genome (i.e ~ 20,000 genes).
type. What kind of data is this column?typetype column.forcats::infreq() is your friendy aesthetic to type.stat = "count" in geom_text().theme_classic() is usually betterlabs() give some contextscales::label_comma() does it for you# y mapping to get flipped axis
ggplot(palmerpenguins::penguins, aes(y = species)) +
# geom_bar does count how many penguins are observed for each species
geom_bar(fill = "pink") +
# add count in bars, label is counting, and horizontal adjust is 1.1 to be just before bar height
geom_text(aes(label = scales::comma(after_stat(count))), hjust = 1.1, stat = "count") +
# light theme
theme_classic() +
labs(y = NULL)
The best columns for transcripts and genes are transcripts_id and gene_id. By distribution, we expect the density or binned histogram of the univariate number of transcript per gene. Depending on the distribution shape, log-transforming the counts is an option
Now, let’s look at genes only.
gencode to restrict rows where type is gene, save as genesgene_type for genes in the same way you did the type. Sorting the counts.You can restrict to counts > 500 to avoid showing all the smaller categories
We can observe quite many unique gene_type, even restriction to the ones that are present > 500 times.
gene_type to create 4 meta-categoriesAs follows:
fct_collapse() combine with fct_other() from the forcats package are good helpers. See below toy example of collapsing levels with “RNA”, then “c” and “d” and finally group the rest in a “rest” level.
x <- tibble(type = c("rRNA", "tRNA", "tRNA", "c", "c", "d", "e", "f"))
# str_subset, will subset a string vector based on a pattern, here: contains "RNA"
x_rna <- str_subset(x$type, "RNA")
x |>
# collapse the rna together (tRNA, rRNA)
mutate(meta_type = fct_collapse(type, RNA = x_rna,
# collapse c and d
cd = c("c", "d")) |>
# the rest is lumped into a "rest" category
fct_other(keep = c("ab", "cd"),
other_level = "rest"))# A tibble: 8 × 2
type meta_type
<chr> <fct>
1 rRNA rest
2 tRNA rest
3 tRNA rest
4 c cd
5 c cd
6 d cd
7 e rest
8 f rest Are the genes uniformly distributed across the genome? To address this question, we want to count the number of genes per chromosome and normalizing by the chromosome length.
You can use this file from USCS to get the chromosome sizes.
chr <- c("chr1", "chr2", "MT", "chr6_GL000255v2_alt", "chr16_KI270728v1_random")
# accessory chromosomes have an underscore
str_subset(chr, "_")[1] "chr6_GL000255v2_alt" "chr16_KI270728v1_random"# use the negate argument to keep the classic chrs
str_subset(chr, "_", negate = TRUE)[1] "chr1" "chr2" "MT" How the gene sizes distribution look like? Here we are several questions, we are need to split sequentially. Also, we keep the 4 meta-categories previously obtained since they represent different kind of genes.
width is directly giving you the length in base pairs.annotation_logticks()alpha parameterNow we want to compute the gene sizes, excluding the introns. Thus, we need to sum up the exons sizes. Of note, exon in the type terminology includes also the UTR. To make things easier, we ignore this inclusion and work on the exon whatever it encompasses.
typetranscript_lengthtranscript_lengthWhat can make the protein coding bigger then? Maybe the introns. This computation is challenging since we don’t have an easy intron item.
In the type column associated with the width column that directly provide the length. The idea is then to work per transcript, and subtract the sum of all exons to the transcript length.
Your workflow could be:
transcript_id)cowplot::plot_grid() or patchwork are 2 possibilities to compose ggplots in one figure. They could share the same legend, as the filling of densities per meta-categories is identical.