RNA-seq Analysis in R
Annotation and Visualisation of RNA-seq results
Mark Dunning
February 2020
Original Authors: Belinda Phipson, Anna Trigos, Matt Ritchie, Maria Doyle, Harriet Dashnow, Charity Law, Stephane Ballereau, Oscar Rueda, Ashley Sawle Based on the course RNAseq analysis in R delivered on May 11/12th 2016
Resources and data files
This material has been created using the following resources:
http://www.statsci.org/smyth/pubs/QLedgeRPreprint.pdf [@Lun2016]
http://monashbioinformaticsplatform.github.io/RNAseq-DE-analysis-with-R/99-RNAseq_DE_analysis_with_R.html
http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html https://bioconductor.github.io/BiocWorkshops/rna-seq-data-analysis-with-deseq2.html
Before starting this section, we will make sure we have all the relevant objects from the Differential Expression analysis present.
suppressPackageStartupMessages(library(DESeq2))
load("Robjects/DE.Rdata")
load("Robjects/preprocessing.Rdata")Overview
- Visualising DE results
- Getting annotation using Bioconductor databases
- Getting annotation using BiomaRt
- Customising RNA-seq plots with
ggplot2 - Retrieving gene models
We can now have a list of genes ordered according to their evidence for being differentially-expressed.
library(dplyr)
library(tibble)
results.status <- as.data.frame(results(de.mf,contrast=c("Status","lactation","virgin"))) %>%
rownames_to_column("ENSEMBL")
results.ordered <- arrange(results.status, padj)
head(results.ordered) ENSEMBL baseMean log2FoldChange lfcSE stat
1 ENSMUSG00000000381 398943.44 9.779005 0.4261732 22.94608
2 ENSMUSG00000061937 1134857.22 8.204288 0.4158950 19.72683
3 ENSMUSG00000026417 11320.43 6.589672 0.3634068 18.13304
4 ENSMUSG00000022491 514850.04 9.589474 0.5336395 17.96995
5 ENSMUSG00000061388 30555.67 10.172377 0.5858320 17.36398
6 ENSMUSG00000024903 60039.52 8.772293 0.5388825 16.27868
pvalue padj
1 1.612341e-116 2.785481e-112
2 1.268851e-86 1.096033e-82
3 1.748099e-73 1.006672e-69
4 3.350312e-72 1.447000e-68
5 1.546170e-67 5.342325e-64
6 1.398617e-59 4.027083e-56In DESeq2, the function plotMA shows the log2 fold changes attributable to a given variable over the mean of normalized counts for all the samples in the DESeqDataSet. Points will be colored red if the adjusted p value is less than 0.1. Points which fall out of the window are plotted as open triangles pointing either up or down.
The log2 fold change for a particular comparison is plotted on the y-axis and the average of the counts normalized by size factor is shown on the x-axis (“M” for minus, because a log ratio is equal to log minus log, and “A” for average). Each gene is represented with a dot. Genes with an adjusted p value below a threshold (here 0.1, the default) are shown in red.
plotMA(results(de.mf,contrast=c("Status","lactation","virgin")))
Note You may see an error message when trying to make the above MA plot. This could be because both limma and DESeq2 have a function called plotMA, and R can sometimes pick the wrong function. To explictly use the DESeq2 function you can use:-
DESeq2::plotMA(results(de.mf,contrast=c("Status","lactation","virgin")))
MA-plots often display a fanning-effect at the left-hand side (genes with low numbers of counts) due to the high variability of the measurements for these genes. For more informative visualization and more accurate ranking of genes by effect size (the log fold change may sometimes be referred to as an effect size), the DESeq2 authors recommend “shrinking” the log fold-changes which is available in DESeq2’s lfcShrink function. This results in more stable fold change values. The p-values are unaffected.
res_LvsV <- lfcShrink(de.mf,contrast=c("Status","lactation","virgin"))
DESeq2::plotMA(res_LvsV)
We will re-define our results object to use these new fold-changes.
results.ordered <- as.data.frame(res_LvsV) %>%
rownames_to_column("ENSEMBL") %>%
arrange(padj)
head(results.ordered) ENSEMBL baseMean log2FoldChange lfcSE stat
1 ENSMUSG00000000381 398943.44 8.722325 0.3815879 22.94608
2 ENSMUSG00000061937 1134857.22 7.351405 0.3761569 19.72683
3 ENSMUSG00000026417 11320.43 6.078220 0.3332017 18.13304
4 ENSMUSG00000022491 514850.04 7.765956 0.4539439 17.96995
5 ENSMUSG00000061388 30555.67 8.282907 0.4740455 17.36398
6 ENSMUSG00000024903 60039.52 7.312158 0.4522371 16.27868
pvalue padj
1 1.612341e-116 2.785481e-112
2 1.268851e-86 1.096033e-82
3 1.748099e-73 1.006672e-69
4 3.350312e-72 1.447000e-68
5 1.546170e-67 5.342325e-64
6 1.398617e-59 4.027083e-56Another common plot for displaying the results of a differential expression analysis is a volcano plot
library(ggplot2)
results.ordered %>%
ggplot(aes(x = log2FoldChange, y = -log10(padj))) + geom_point()
It can also be useful to examine the counts of reads for a single gene across the groups. A simple function for making this plot is plotCounts, which normalizes counts by sequencing depth and adds a pseudocount of 1/2 to allow for log scale plotting. The counts are grouped by the variables in intgroup, where more than one variable can be specified. Here we specify the gene which had the smallest p value from the results table created above. You can select the gene to plot by rowname or by numeric index:-
plotCounts(dds, "ENSMUSG00000000381",intgroup = c("Status"))
If we want greater control over how to visualise the data, we can use the plotCounts function to return the count data, but not actually produce the plot:-
plotCounts(dds, "ENSMUSG00000000381",intgroup = c("Status"),returnData=TRUE) count Status
SRR1552444 1.569482e+03 virgin
SRR1552445 1.671109e+03 virgin
SRR1552446 9.802440e+04 pregnancy
SRR1552447 2.205304e+05 pregnancy
SRR1552448 2.292659e+06 lactation
SRR1552449 2.137836e+06 lactation
SRR1552450 2.427758e+01 virgin
SRR1552451 3.357483e+01 virgin
SRR1552452 3.410069e+03 pregnancy
SRR1552453 4.660197e+03 pregnancy
SRR1552454 1.368587e+04 lactation
SRR1552455 1.322326e+04 lactationChallenge 1
- Use the option
returnData=TRUEto get a data frame containing the counts ofENSMUSG00000000381in the different development stages. Visualise these data usingggplot2(see plot A below).- Repeat the volcano plot from above, but use a different colour to indicate which genes are significant with an adjusted p-value less than 0.05. See plot B below
- (Optional) The argument
intgroup=can be used to retrieve and plot data from multiple variables of interest in the data. Use the valueintgroup=c("Status","CellType")and compare the counts between different cell types and status. See plot C below. HINT: To get the counts on the same scale as displayed by the plotCounts function you will need to add+scale_y_log10in your ggplot2 code
However, it is hard to assess the biological significance of such a gene without more information about . To perform such a task we need to map between the identifiers we have in the DESeq2 output and more familiar names.
Adding annotation to the DESeq2 results
There are a number of ways to add annotation, but we will demonstrate how to do this using the org.Mm.eg.db package. This package is one of several organism-level packages which are re-built every 6 months. These packages are listed on the annotation section of the Bioconductor, and are installed in the same way as regular Bioconductor packages. An alternative approach is to use biomaRt, an interface to the BioMart resource. BioMart is much more comprehensive, but the organism packages fit better into the Bioconductor workflow.
### Only execute when you need to install the package
install.packages("BiocManager")
BiocManager::install("org.Mm.eg.db")
# For Human
BiocManager::install("org.Hs.eg.db")The packages are larger in size that Bioconductor software pacakges, but essentially they are databases that can be used to make offline queries.
library(org.Mm.eg.db)First we need to decide what information we want. In order to see what we can extract we can run the columns function on the annotation database.
columns(org.Mm.eg.db) [1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT"
[5] "ENSEMBLTRANS" "ENTREZID" "ENZYME" "EVIDENCE"
[9] "EVIDENCEALL" "GENENAME" "GO" "GOALL"
[13] "IPI" "MGI" "ONTOLOGY" "ONTOLOGYALL"
[17] "PATH" "PFAM" "PMID" "PROSITE"
[21] "REFSEQ" "SYMBOL" "UNIGENE" "UNIPROT" We are going to filter the database by a key or set of keys in order to extract the information we want. Valid names for the key can be retrieved with the keytypes function.
keytypes(org.Mm.eg.db) [1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT"
[5] "ENSEMBLTRANS" "ENTREZID" "ENZYME" "EVIDENCE"
[9] "EVIDENCEALL" "GENENAME" "GO" "GOALL"
[13] "IPI" "MGI" "ONTOLOGY" "ONTOLOGYALL"
[17] "PATH" "PFAM" "PMID" "PROSITE"
[21] "REFSEQ" "SYMBOL" "UNIGENE" "UNIPROT" We should see ENSEMBL, which is the type of key we are going to use in this case. If we are unsure what values are acceptable for the key, we can check what keys are valid with keys
keys(org.Mm.eg.db, keytype="ENSEMBL")[1:10] [1] "ENSMUSG00000030359" "ENSMUSG00000020804" "ENSMUSG00000025375"
[4] "ENSMUSG00000015243" "ENSMUSG00000028125" "ENSMUSG00000026944"
[7] "ENSMUSG00000031333" "ENSMUSG00000024030" "ENSMUSG00000058835"
[10] "ENSMUSG00000026842"For the top gene in our analysis the call to the function would be:-
select(org.Mm.eg.db, keys="ENSMUSG00000000381",
keytype = "ENSEMBL",columns=c("SYMBOL","GENENAME")
)Unfortunately, the authors of dplyr and AnnotationDbi have both decided to use the name select in their packages. To avoid confusion, the following code is sometimes used:-
AnnotationDbi::select(org.Mm.eg.db, keys="ENSMUSG00000000381",keytype = "ENSEMBL",columns=c("SYMBOL","GENENAME")) ENSEMBL SYMBOL GENENAME
1 ENSMUSG00000000381 Wap whey acidic proteinTo annotate our results, we definitely want gene symbols and perhaps the full gene name. Let’s build up our annotation information into a new data frame using the select function.
anno <- AnnotationDbi::select(org.Mm.eg.db,keys=results.ordered$ENSEMBL,
columns=c("SYMBOL","GENENAME"),
keytype="ENSEMBL")
# Have a look at the annotation
head(anno) ENSEMBL SYMBOL
1 ENSMUSG00000000381 Wap
2 ENSMUSG00000061937 Csn1s2a
3 ENSMUSG00000026417 Pigr
4 ENSMUSG00000022491 Glycam1
5 ENSMUSG00000061388 Csn1s2b
6 ENSMUSG00000024903 Lao1
GENENAME
1 whey acidic protein
2 casein alpha s2-like A
3 polymeric immunoglobulin receptor
4 glycosylation dependent cell adhesion molecule 1
5 casein alpha s2-like B
6 L-amino acid oxidase 1However, we have a problem that the resulting data frame has more rows than our results table. This is due to the one-to-many relationships that often occur when mapping between various identifiers.
dim(anno)[1] 18071 3dim(results.ordered)[1] 17973 7Such duplicated entries can be identified using the duplicated function.
dup_ids <- anno$ENSEMBL[duplicated(anno$ENSEMBL)]
filter(anno, ENSEMBL %in% dup_ids) %>%
arrange(ENSEMBL) %>% head ENSEMBL SYMBOL
1 ENSMUSG00000000486 Sept1
2 ENSMUSG00000000486 Gm4532
3 ENSMUSG00000000562 Adora3
4 ENSMUSG00000000562 Tmigd3
5 ENSMUSG00000001768 Rin2
6 ENSMUSG00000001768 BC039771
GENENAME
1 septin 1
2 predicted gene 4532
3 adenosine A3 receptor
4 transmembrane and immunoglobulin domain containing 3
5 Ras and Rab interactor 2
6 cDNA sequence BC039771Fortunately, there are not too many so hopefully we won’t lose too much information if we discard the entries that are duplicated. The first occurence of the duplicated ID will still be included in the table.
anno <- AnnotationDbi::select(org.Mm.eg.db,keys=results.ordered$ENSEMBL,
columns=c("ENSEMBL","SYMBOL","GENENAME","ENTREZID"),
keytype="ENSEMBL") %>%
filter(!duplicated(ENSEMBL))
dim(anno)[1] 17973 4We can bind in the annotation information to the results data frame.
results.annotated <- left_join(results.ordered, anno,by="ENSEMBL")
head(results.annotated) ENSEMBL baseMean log2FoldChange lfcSE stat
1 ENSMUSG00000000381 398943.44 8.722325 0.3815879 22.94608
2 ENSMUSG00000061937 1134857.22 7.351405 0.3761569 19.72683
3 ENSMUSG00000026417 11320.43 6.078220 0.3332017 18.13304
4 ENSMUSG00000022491 514850.04 7.765956 0.4539439 17.96995
5 ENSMUSG00000061388 30555.67 8.282907 0.4740455 17.36398
6 ENSMUSG00000024903 60039.52 7.312158 0.4522371 16.27868
pvalue padj SYMBOL
1 1.612341e-116 2.785481e-112 Wap
2 1.268851e-86 1.096033e-82 Csn1s2a
3 1.748099e-73 1.006672e-69 Pigr
4 3.350312e-72 1.447000e-68 Glycam1
5 1.546170e-67 5.342325e-64 Csn1s2b
6 1.398617e-59 4.027083e-56 Lao1
GENENAME ENTREZID
1 whey acidic protein 22373
2 casein alpha s2-like A 12993
3 polymeric immunoglobulin receptor 18703
4 glycosylation dependent cell adhesion molecule 1 14663
5 casein alpha s2-like B 12992
6 L-amino acid oxidase 1 100470We can save the results table using the write.csv function, which writes the results out to a csv file that you can open in excel.
write.csv(results.annotated,file="virgin_vs_lactation_DESeq_annotated.csv",row.names=FALSE)We have already seen the use of a heatmap as a quality assessment tool to visualise the relationship between samples in an experiment. Another common use-case for such a plot is to visualise the results of a differential expression analysis.
Here we will take the top 10 genes from the differential expression analysis and produce a heatmap with the pheatmap package. The default colour palette goes from low expression in blue to high expression in red, which is a good alternative to the traditional red/green heatmaps which are not suitable for those with forms of colour-blindness.
The counts we are visualising are the variance-stablised counts, which are more appropriate for visualisation.
library(pheatmap)
top_genes <- results.annotated$ENSEMBL[1:10]
vsd <- vst(dds)
pheatmap(assay(vsd)[top_genes,])
The heatmap is more informative if we add colours underneath the sample dendrogram to indicate which sample group each sample belongs to. This we can do by creating a data frame containing metadata for each of the samples in our dataset. With the DESeq2 workflow we have already created such a data frame. We have to make sure the the rownames of the data frame are the same as the column names of the counts matrix.
sampleInfo <- as.data.frame(colData(dds)[,c("Status","CellType")])
pheatmap(assay(vsd)[top_genes,],
annotation_col = sampleInfo)
Any plot we create in RStudio can be saved as a png or pdf file. We use the png or pdf function to create a file for the plot to be saved into and run the rest of the code as normal. The plot does not get displayed in RStudio, but printed to the specified file.
png("heatmap_top10_genes.png",width=800,height=800)
pheatmap(assay(vsd)[top_genes,],
annotation_col = sampleInfo)
# dev.off()Challenge 2
- Repeat the same heatmap as above, but for the top 100 most differentially-expressed genes between pregnant and luminal
- Change the plot so that gene names are displayed rather than Ensembl IDs
- Save the plot to a pdf file HINT: check the help for the
pheatmapfunction to see how column and row labels can be changed
Accessing the sample or gene clusters
The heatmap displays relationships between samples and genes in our study as a useful visualisation. In this example we can easily identify which samples are most similar based on their expression patterns. However, for larger dataset this may be more problematic. We can extract data the sample relationships about if we manually perform the clustering steps used by pheatmap. First is to cluster the samples with the default distance matrix and clustering algorithms.
mat <- assay(vsd)[top_genes,]
## Calculate the distance matrix between samples
d_samples <- dist(t(mat))
plot(hclust(d_samples))
rect.hclust(hclust(d_samples),k=2)
We can then “cut” the dendrogram to give a set number of clusters. Each sample has been assigned a label of 1 or 2 depending on which cluster it belongs to.
clusters <- cutree(hclust(d_samples),k = 2)
clustersSRR1552444 SRR1552445 SRR1552446 SRR1552447 SRR1552448 SRR1552449
1 1 2 2 2 2
SRR1552450 SRR1552451 SRR1552452 SRR1552453 SRR1552454 SRR1552455
1 1 1 1 1 1 The groupings could then be tabulated against with the sample metadata to see if particular biological groups are associated with the new clusters we have identified.
table(clusters, colData(dds)$Status)
clusters lactation pregnancy virgin
1 2 2 4
2 2 2 0Adding gene names to a volcano plot
Now that we have an annotated table of results, we can add the gene names to some of the other plots we have created. This should be straightforward as ggplot2 has a label aesthetic that can be mapped to columns in a data frame. The geom_text plot will then display the labels. However, the following plot is a bit crowded.
## Not a good idea to run this!!
results.annotated %>%
ggplot(aes(x = log2FoldChange, y = -log10(padj), label=SYMBOL)) + geom_point() + geom_text()
The problem here is that ggplot2 is trying to label every point with a name; not quite what we want. The trick is to create a label that is blank for most genes and only labels the points we are interested in. The ifelse function in R is a convenient way to set the entries in a vector based on a logical expression. In this case, make the values in Label the same as the gene symbol if the gene is in our list of “top genes”. Otherwise, points get labeled with a blank string "".
For clarity, we also make the points slightly transparent and use a different colour for the text.
N <- 10
top_genes <- results.annotated$ENSEMBL[1:N]
results.annotated %>%
mutate(Label = ifelse(ENSEMBL %in% top_genes, SYMBOL, "")) %>%
ggplot(aes(x = log2FoldChange, y = -log10(padj), label=Label)) + geom_point(alpha=0.4) + geom_text(col="blue")
Finally, a slightly better positioning of text is given by the ggrepel package.
if(!require(ggrepel)) install.packages("ggrepel")
results.annotated %>%
mutate(Label = ifelse(ENSEMBL %in% top_genes, SYMBOL, "")) %>%
ggplot(aes(x = log2FoldChange, y = -log10(padj), label=Label)) + geom_point(alpha=0.4) + geom_text_repel(col="blue")
Annotation with the biomaRt resource
The Bioconductor package have the convenience of being able to make queries offline. However, they are only available for certain organisms. If your organism does not have an org.XX.eg.db package listed on the Bioconductor annotation page (http://bioconductor.org/packages/release/BiocViews.html#___AnnotationData), an alternative is to use biomaRt which provides an interface to the popular biomart annotation resource.
The first step is to find the name of a database that you want to connect to
library(biomaRt)
listMarts() biomart version
1 ENSEMBL_MART_ENSEMBL Ensembl Genes 99
2 ENSEMBL_MART_MOUSE Mouse strains 99
3 ENSEMBL_MART_SNP Ensembl Variation 99
4 ENSEMBL_MART_FUNCGEN Ensembl Regulation 99ensembl=useMart("ENSEMBL_MART_ENSEMBL")
# list the available datasets (species). Replace mouse with the name of your organism
listDatasets(ensembl) %>% filter(grepl("Mouse",description)) dataset description version
1 mmurinus_gene_ensembl Mouse Lemur genes (Mmur_3.0) Mmur_3.0
2 mmusculus_gene_ensembl Mouse genes (GRCm38.p6) GRCm38.p6ensembl = useDataset("mmusculus_gene_ensembl", mart=ensembl)Queries to biomaRt are constructed in a similar way to the queries we performed with the org.Mm.eg.db package. Instead of keys we have filters, and instead of columns we have attributes. The list of acceptable values is much more comprehensive that for the org.Mm.eg.db package.
listFilters(ensembl) %>%
filter(grepl("ensembl",name)) name
1 with_clone_based_ensembl_gene
2 with_clone_based_ensembl_transcript
3 ensembl_gene_id
4 ensembl_gene_id_version
5 ensembl_transcript_id
6 ensembl_transcript_id_version
7 ensembl_peptide_id
8 ensembl_peptide_id_version
9 ensembl_exon_id
10 clone_based_ensembl_gene
11 clone_based_ensembl_transcript
description
1 With Clone-based (Ensembl) gene ID(s)
2 With Clone-based (Ensembl) transcript ID(s)
3 Gene stable ID(s) [e.g. ENSMUSG00000000001]
4 Gene stable ID(s) with version [e.g. ENSMUSG00000000001.4]
5 Transcript stable ID(s) [e.g. ENSMUST00000000001]
6 Transcript stable ID(s) with version [e.g. ENSMUST00000000001.4]
7 Protein stable ID(s) [e.g. ENSMUSP00000000001]
8 Protein stable ID(s) with version [e.g. ENSMUSP00000000001.4]
9 Exon ID(s) [e.g. ENSMUSE00000097910]
10 Clone-based (Ensembl) gene ID(s) [e.g. AC015535.1]
11 Clone-based (Ensembl) transcript ID(s) [e.g. AC015535.1-201]listAttributes(ensembl) %>%
filter(grepl("gene",name))An advantage over the org.. packages is that positional information can be retrieved
attributeNames <- c('ensembl_gene_id', 'entrezgene_id', 'external_gene_name')
getBM(attributes = attributeNames,
filters = "ensembl_gene_id",
values=top_genes,
mart=ensembl) ensembl_gene_id entrezgene_id external_gene_name
1 ENSMUSG00000000381 22373 Wap
2 ENSMUSG00000022491 14663 Glycam1
3 ENSMUSG00000024331 13506 Dsc2
4 ENSMUSG00000024903 100470 Lao1
5 ENSMUSG00000026417 18703 Pigr
6 ENSMUSG00000040260 76441 Daam2
7 ENSMUSG00000061388 12992 Csn1s2b
8 ENSMUSG00000061937 12993 Csn1s2a
9 ENSMUSG00000063275 30963 Hacd1
10 ENSMUSG00000070702 12990 Csn1s1Challenge 3
- Use biomaRt to create an data frame containing the entrezgene, gene symbol and genomic coordinates (chromosome, start, end) for the Ensembl IDs in the DESeq2 results
- Remove duplicates entries from the new data frame
- Join the biomaRt annotation to the DESeq2 results to produce a data frame with differential expression results and annotation
- Write the joined data frame to a csv file
Obtaining gene models with Bioconductor
Using biomaRt as above allows us to retrieve the genomic coordinates of a given gene. If we want more-comprehensive information about the structue of a gene (and indeed all genes in the transcriptome), we can use one of several pre-built transcript database packages.
Retrieving the coordinates for a particular gene (or set of genes) uses the same select function as for the org.Mm.eg.db package, but checking for different columns and keytypes.
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
columns(txdb) [1] "CDSCHROM" "CDSEND" "CDSID" "CDSNAME" "CDSPHASE"
[6] "CDSSTART" "CDSSTRAND" "EXONCHROM" "EXONEND" "EXONID"
[11] "EXONNAME" "EXONRANK" "EXONSTART" "EXONSTRAND" "GENEID"
[16] "TXCHROM" "TXEND" "TXID" "TXNAME" "TXSTART"
[21] "TXSTRAND" "TXTYPE" AnnotationDbi::select(txdb, columns = c("EXONID","EXONSTART","EXONCHROM"),
keys="22373",
keytype = "GENEID") GENEID EXONID EXONCHROM EXONSTART
1 22373 163805 chr11 6638535
2 22373 163804 chr11 6637190
3 22373 163803 chr11 6636710
4 22373 163802 chr11 6635483Alternatively we can grab the coordinates of all genes in a single object and then subset accordingly. This allows us to perform all kinds of subset operation using the GenomicFeatures framework in Bioconductor.
exons <- exonsBy(txdb,"gene")
exonsGRangesList object of length 24421:
$100009600
GRanges object with 7 ranges and 2 metadata columns:
seqnames ranges strand | exon_id exon_name
<Rle> <IRanges> <Rle> | <integer> <character>
[1] chr9 21062393-21062717 - | 134805 <NA>
[2] chr9 21062894-21062987 - | 134806 <NA>
[3] chr9 21063314-21063396 - | 134807 <NA>
[4] chr9 21066024-21066377 - | 134808 <NA>
[5] chr9 21066940-21067925 - | 134809 <NA>
[6] chr9 21068030-21068117 - | 134810 <NA>
[7] chr9 21073075-21075496 - | 134812 <NA>
$100009609
GRanges object with 6 ranges and 2 metadata columns:
seqnames ranges strand | exon_id exon_name
[1] chr7 84940169-84941088 - | 110194 <NA>
[2] chr7 84943141-84943264 - | 110195 <NA>
[3] chr7 84943504-84943722 - | 110196 <NA>
[4] chr7 84946200-84947000 - | 110197 <NA>
[5] chr7 84947372-84947651 - | 110198 <NA>
[6] chr7 84963816-84964009 - | 110199 <NA>
$100009614
GRanges object with 1 range and 2 metadata columns:
seqnames ranges strand | exon_id exon_name
[1] chr10 77711446-77712009 + | 144278 <NA>
...
<24418 more elements>
-------
seqinfo: 66 sequences (1 circular) from mm10 genomeexons[["22373"]]GRanges object with 4 ranges and 2 metadata columns:
seqnames ranges strand | exon_id exon_name
<Rle> <IRanges> <Rle> | <integer> <character>
[1] chr11 6635483-6635629 - | 163802 <NA>
[2] chr11 6636710-6636874 - | 163803 <NA>
[3] chr11 6637190-6637324 - | 163804 <NA>
[4] chr11 6638535-6638649 - | 163805 <NA>
-------
seqinfo: 66 sequences (1 circular) from mm10 genomeInteractive graphs and tables
It is often useful to be able to explore our results in an interactive manner; searching for our favourite genes of interest and plotting on-the-fly whether they are statistically significant in our dataset or not.
Such a visualisation is possible with the Glimma Bioconductor package.
It takes our DESeq2 results object, annotation table and normalized counts, and produces a HTML page including a sortable results table, MA-plot and scatter plot. Particular genes can be searched among the table and their expression patterns can be displayed. Alternatively we can click on particular point in the plot and display their stats.
results <- results(de.mf,contrast=c("Status","lactation","virgin"))
## Repeat the annotation, as the previous annotation table was created using an ordered results table
anno <- AnnotationDbi::select(org.Mm.eg.db,keys=rownames(results),
columns=c("SYMBOL","GENENAME"),
keytype="ENSEMBL") %>%
filter(!duplicated(ENSEMBL))## Make sure we have normalised counts before proceeding
dds <- estimateSizeFactors(dds)## Load the Glimma package and create the report
library(Glimma)
glMDPlot(results,
anno,
groups = colData(dds)$Status,
counts = counts(dds,normalized=TRUE),
transform = TRUE,
side.main = "ENSEMBL")