RNAseq in Bioconductor exercises (part 4)

In todays session we will work with some of the RNAseq data of adult mouse tissues from Bing Ren’s lab, Liver and Heart.

• More information on liver data can be found here

• More information on heart data can be found here

• More information on Kidney data can be found here

• Precounted RNAseq reads in genes for these tissues can be found as an R data object in data/gC_TissueFull.RData.

Exercises

1. Differential Transcript Usage

We have run the DEXseq workflow to compare Heart and Liver. Load the DEXseqResults object called dxr1 (found in data/dxr_HeartVsLiver.RData) and make a plot of the differential splicing for the Atp2a2 gene.

library(DEXSeq)
load("data/dxr_HeartVsLiver.RData")

dxr1DF <- as.data.frame(dxr1)
dxr1DF <- dxr1DF[order(dxr1DF$pvalue),]
library(org.Mm.eg.db)
eToSym <- AnnotationDbi::select(org.Mm.eg.db,
                 keys = unique(dxr1DF[,1]),
                 keytype = "ENTREZID",
                 columns="SYMBOL")

## 'select()' returned 1:1 mapping between keys and columns

annotatedRes <- merge(eToSym,dxr1DF,
                      by.x=1,
                      by.y=1,all=FALSE)
annotatedRes <- annotatedRes[order(annotatedRes$pvalue),]
annotatedRes[1:3,]

##       ENTREZID SYMBOL featureID exonBaseMean  dispersion     stat        pvalue
## 10613    22003   Tpm1      E001     117.4021 0.009979611 790.0469 7.871657e-174
## 3297     12870     Cp      E026      72.2033 0.003310883 715.7727 1.111406e-157
## 2345     11938 Atp2a2      E002     651.1539 0.014388628 389.7594  9.338217e-87
##                padj    Heart    Liver log2fold_Liver_Heart genomicData.seqnames
## 10613 1.522772e-169 1.030876 3.134205             7.127092                 chr9
## 3297  1.075007e-153 2.429502 1.079348            -4.605220                 chr3
## 2345   6.021594e-83 2.146096 3.454840             4.357391                 chr5
##       genomicData.start genomicData.end genomicData.width genomicData.strand
## 10613          67022593        67023441               849                  -
## 3297           20008540        20009750              1211                  +
## 2345          122455893       122457302              1410                  -
##       countData.Sorted_Heart_1.bam countData.Sorted_Heart_2.bam
## 10613                          174                          139
## 3297                           156                          112
## 2345                           985                          829
##       countData.Sorted_Liver_1.bam countData.Sorted_Liver_2.bam  transcripts
## 10613                           80                           97 uc009qfo....
## 3297                            25                           21   uc008osc.2
## 2345                           505                          403   uc008zlj.2

plotDEXSeq(dxr1,11938,fitExpToVar = "tissue",displayTranscripts = TRUE)

2. Differential Transcript Usage

Load the RangedSummarizedExperiment object called tissueExonCounts (found in data/RSE_HeartAndLiver.RData). Perform a differential exon analysis. Make a plot of the gene containing the most significantly differentially used exon.

load(file="data/RSE_HeartAndLiver.RData")

ddxTissue <- DEXSeqDataSetFromSE(tissueExonCounts,
                      design= ~ sample + exon + tissue:exon)

## converting counts to integer mode

ddxTissue <- DEXSeq(ddxTissue,fitExpToVar = "tissue")

## -- note: fitType='parametric', but the dispersion trend was not well captured by the
##    function: y = a/x + b, and a local regression fit was automatically substituted.
##    specify fitType='local' or 'mean' to avoid this message next time.

ddxTissue <- ddxTissue[order(ddxTissue$pvalue),]
plotDEXSeq(ddxTissue,16561,fitExpToVar = "tissue",displayTranscripts = TRUE)