Contents

Included is an example of how to download and analyze expression data from SRA study SRP032798. The data come from human breast cancer samples, and we compare the transcriptomes of TNBC samples (triple negative breast cancer) and HER2-positive breast cancer samples (breast cancer type that tests positive for a protein called human epidermal growth factor receptor 2). Code here demonstrates how to carry out differential expression analyses on gene, exon, junction, and differential expressed region (DER) levels within a single study using limma and voom. We test for concordance among the results of each analysis and demonstrate how to carry out gene ontology analysis using topGO to characterize top hits from differential expression analyses.

We first load the required packages.

## load libraries
library('recount')
library('SummarizedExperiment')
library('limma')
library('edgeR')
library('qvalue')
library('topGO')
library('matrixStats')
library('RSkittleBrewer')
library('derfinder')

1 Gene level analysis

We first download the project of interest (SRP032798), obtaining expression data for the study of interest. We obtain summaries of the number of samples and genes included using colData() and rowData(), respectively.

## Find the project of interest (SRP032789), e.g. with parts of the abstract
project_info <- abstract_search('To define the digital transcriptome of three breast cancer')

## Explore information
project_info
##     number_samples species
## 865             20   human
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               abstract
## 865 Goal: To define the digital transcriptome of three breast cancer subtypes (TNBC, Non-TNBC, and HER2-positive) using RNA-sequencing technology. To elucidate differentially expressed known and novel transcripts, alternatively spliced genes and differential isoforms and lastly expressed variants in our dataset. Method: Dr. Suzanne Fuqua (Baylor College of Medicine) provided the human breast cancer tissue RNA samples. All of the human samples were used in accordance with the IRB procedures of Baylor College of Medicine. The breast tumour types, TNBC, Non-TNBC and HER2-positive, were classified on the basis of immunohistochemical and RT-qPCR classification. Results: Comparative transcriptomic analyses elucidated differentially expressed transcripts between the three breast cancer groups, identifying several new modulators of breast cancer. We discovered subtype specific differentially spliced genes and splice isoforms not previously recognized in human transcriptome. Further, we showed that exon skip and intron retention are predominant splice events in breast cancer. In addition, we found that differential expression of primary transcripts and promoter switching are significantly deregulated in breast cancer compared to normal breast.  We also report novel expressed variants, allelic prevalence and abundance, and coexpression with other variation, and splicing signatures. Additionally we describe novel SNPs and INDELs in cancer relevant genes with no prior reported association of point mutations with cancer Overall design: mRNA profiles of 17 breast tumor samples of three different subtypes (TNBC, non-TNBC and HER2-positive) and normal human breast organoids (epithelium) samples (NBS) were sequenced using Illumina HiSeq.
##       project
## 865 SRP032789
## Browse the project at SRA
browse_study(project_info$project)

## Download the gene level RangedSummarizedExperiment data
if(!file.exists(file.path('SRP032789', 'rse_gene.Rdata'))) {
    download_study(project_info$project)
}

## Load the data
load(file.path(project_info$project, 'rse_gene.Rdata'))
rse_gene
## class: RangedSummarizedExperiment 
## dim: 58037 20 
## metadata(0):
## assays(1): counts
## rownames(58037): ENSG00000000003.14 ENSG00000000005.5 ...
##   ENSG00000283698.1 ENSG00000283699.1
## rowData names(3): gene_id bp_length symbol
## colnames(20): SRR1027171 SRR1027173 ... SRR1027190 SRR1027172
## colData names(21): project sample ... title characteristics
## This is the phenotype data provided by the recount project
colData(rse_gene)
## DataFrame with 20 rows and 21 columns
##                project      sample  experiment         run
##            <character> <character> <character> <character>
## SRR1027171   SRP032789   SRS500214   SRX374850  SRR1027171
## SRR1027173   SRP032789   SRS500216   SRX374852  SRR1027173
## SRR1027174   SRP032789   SRS500217   SRX374853  SRR1027174
## SRR1027175   SRP032789   SRS500218   SRX374854  SRR1027175
## SRR1027176   SRP032789   SRS500219   SRX374855  SRR1027176
## ...                ...         ...         ...         ...
## SRR1027187   SRP032789   SRS500230   SRX374866  SRR1027187
## SRR1027188   SRP032789   SRS500231   SRX374867  SRR1027188
## SRR1027189   SRP032789   SRS500232   SRX374868  SRR1027189
## SRR1027190   SRP032789   SRS500233   SRX374869  SRR1027190
## SRR1027172   SRP032789   SRS500215   SRX374851  SRR1027172
##            read_count_as_reported_by_sra reads_downloaded
##                                <integer>        <integer>
## SRR1027171                      88869444         88869444
## SRR1027173                     107812596        107812596
## SRR1027174                      98563260         98563260
## SRR1027175                      91327892         91327892
## SRR1027176                      96513572         96513572
## ...                                  ...              ...
## SRR1027187                      75260678         75260678
## SRR1027188                      65709192         65709192
## SRR1027189                      65801392         65801392
## SRR1027190                      74356276         74356276
## SRR1027172                      80986440         58902122
##            proportion_of_reads_reported_by_sra_downloaded paired_end
##                                                 <numeric>  <logical>
## SRR1027171                                              1       TRUE
## SRR1027173                                              1       TRUE
## SRR1027174                                              1       TRUE
## SRR1027175                                              1       TRUE
## SRR1027176                                              1       TRUE
## ...                                                   ...        ...
## SRR1027187                                      1.0000000       TRUE
## SRR1027188                                      1.0000000       TRUE
## SRR1027189                                      1.0000000       TRUE
## SRR1027190                                      1.0000000       TRUE
## SRR1027172                                      0.7273084       TRUE
##            sra_misreported_paired_end mapped_read_count        auc
##                             <logical>         <integer>  <numeric>
## SRR1027171                      FALSE          86949307 5082692127
## SRR1027173                      FALSE         104337779 6077034329
## SRR1027174                      FALSE          95271238 5504462845
## SRR1027175                      FALSE          88820239 5150234117
## SRR1027176                      FALSE          93464650 5416681912
## ...                               ...               ...        ...
## SRR1027187                      FALSE          64697612 3567078255
## SRR1027188                      FALSE          65278500 4856453823
## SRR1027189                      FALSE          65328289 4858587600
## SRR1027190                      FALSE          73911898 5501089036
## SRR1027172                      FALSE          57523391 3351013968
##            sharq_beta_tissue sharq_beta_cell_type
##                  <character>          <character>
## SRR1027171            breast                  esc
## SRR1027173            breast                  esc
## SRR1027174            breast                  esc
## SRR1027175            breast                  esc
## SRR1027176            breast                  esc
## ...                      ...                  ...
## SRR1027187            breast                  esc
## SRR1027188            breast                  esc
## SRR1027189            breast                  esc
## SRR1027190            breast                  esc
## SRR1027172            breast                  esc
##            biosample_submission_date biosample_publication_date
##                          <character>                <character>
## SRR1027171   2013-11-07T12:40:22.203    2013-11-08T01:11:17.160
## SRR1027173   2013-11-07T12:40:32.283    2013-11-08T01:11:14.827
## SRR1027174   2013-11-07T12:40:28.283    2013-11-08T01:11:52.283
## SRR1027175   2013-11-07T12:40:34.343    2013-11-08T01:11:15.963
## SRR1027176   2013-11-07T12:40:36.303    2013-11-08T01:11:46.430
## ...                              ...                        ...
## SRR1027187   2013-11-07T12:40:56.180    2013-11-08T01:11:29.587
## SRR1027188   2013-11-07T12:40:58.170    2013-11-08T01:12:06.660
## SRR1027189   2013-11-07T12:40:20.227    2013-11-08T01:11:33.080
## SRR1027190   2013-11-07T12:40:18.090    2013-11-08T01:12:11.320
## SRR1027172   2013-11-07T12:40:26.217    2013-11-08T01:11:45.250
##              biosample_update_date avg_read_length geo_accession
##                        <character>       <integer>   <character>
## SRR1027171 2014-03-07T16:09:38.542             120    GSM1261016
## SRR1027173 2014-03-07T16:09:38.698             120    GSM1261018
## SRR1027174 2014-03-07T16:09:38.637             120    GSM1261019
## SRR1027175 2014-03-07T16:09:38.731             120    GSM1261020
## SRR1027176 2014-03-07T16:09:38.768             120    GSM1261021
## ...                            ...             ...           ...
## SRR1027187 2014-03-07T16:09:39.093             120    GSM1261032
## SRR1027188 2014-03-07T16:09:39.130             150    GSM1261033
## SRR1027189 2014-03-07T16:09:38.498             150    GSM1261034
## SRR1027190 2014-03-07T16:09:38.469             150    GSM1261035
## SRR1027172 2014-03-07T16:09:38.604              87    GSM1261017
##              bigwig_file       title
##              <character> <character>
## SRR1027171 SRR1027171.bw       TNBC1
## SRR1027173 SRR1027173.bw       TNBC3
## SRR1027174 SRR1027174.bw       TNBC4
## SRR1027175 SRR1027175.bw       TNBC5
## SRR1027176 SRR1027176.bw       TNBC6
## ...                  ...         ...
## SRR1027187 SRR1027187.bw      HER2-5
## SRR1027188 SRR1027188.bw        NBS1
## SRR1027189 SRR1027189.bw        NBS2
## SRR1027190 SRR1027190.bw        NBS3
## SRR1027172 SRR1027172.bw       TNBC2
##                                   characteristics
##                                   <CharacterList>
## SRR1027171          tumor type: TNBC Breast Tumor
## SRR1027173          tumor type: TNBC Breast Tumor
## SRR1027174          tumor type: TNBC Breast Tumor
## SRR1027175          tumor type: TNBC Breast Tumor
## SRR1027176          tumor type: TNBC Breast Tumor
## ...                                           ...
## SRR1027187 tumor type: HER2 Positive Breast Tumor
## SRR1027188    tumor type: Normal Breast Organoids
## SRR1027189    tumor type: Normal Breast Organoids
## SRR1027190    tumor type: Normal Breast Organoids
## SRR1027172          tumor type: TNBC Breast Tumor
## At the gene level, the row data includes the names of the genes and
## the sum of the reduced exons widths, which can be used for taking into
## account the gene length.
rowData(rse_gene)
## DataFrame with 58037 rows and 3 columns
##                  gene_id bp_length          symbol
##              <character> <integer> <CharacterList>
## 1     ENSG00000000003.14      4535          TSPAN6
## 2      ENSG00000000005.5      1610            TNMD
## 3     ENSG00000000419.12      1207            DPM1
## 4     ENSG00000000457.13      6883           SCYL3
## 5     ENSG00000000460.16      5967        C1orf112
## ...                  ...       ...             ...
## 58033  ENSG00000283695.1        61              NA
## 58034  ENSG00000283696.1       997              NA
## 58035  ENSG00000283697.1      1184    LOC101928917
## 58036  ENSG00000283698.1       940              NA
## 58037  ENSG00000283699.1        60         MIR4481

Downloaded count data are first scaled to take into account differing coverage between samples. Phenotype data (pheno) are obtained and ordered to match the sample order of the gene expression data (rse_gene). Only those samples that are HER2-positive or TNBC are included for analysis. Prior to differential gene expression analysis, count data are obtained in matrix format and then filtered to only include those genes with greater than five average normalized counts across all samples.

## Scale counts by taking into account the total coverage per sample
rse <- scale_counts(rse_gene)

## Download additional phenotype data from 
## http://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP032789
pheno <- read.table('SraRunTable_SRP032789.txt', sep = '\t', 
    header = TRUE,
    stringsAsFactors = FALSE)

## Obtain correct order for pheno data
pheno <- pheno[match(rse$run, pheno$Run_s), ]
identical(pheno$Run_s, rse$run)
## [1] TRUE
head(cbind(pheno$Run_s, rse$run))
##      [,1]         [,2]        
## [1,] "SRR1027171" "SRR1027171"
## [2,] "SRR1027173" "SRR1027173"
## [3,] "SRR1027174" "SRR1027174"
## [4,] "SRR1027175" "SRR1027175"
## [5,] "SRR1027176" "SRR1027176"
## [6,] "SRR1027177" "SRR1027177"
## Obtain grouping information
colData(rse)$group <- pheno$tumor_type_s
table(colData(rse)$group)   
## 
## HER2 Positive Breast Tumor      Non-TNBC Breast Tumor 
##                          5                          6 
##    Normal Breast Organoids          TNBC Breast Tumor 
##                          3                          6
## Subset data to HER2 and TNBC types
rse <- rse[, rse$group %in% c('HER2 Positive Breast Tumor',
                              'TNBC Breast Tumor')]

## Save filtered rse object
rse_gene_filt <- rse 

## Obtain count matrix
counts <- assays(rse_gene_filt)$counts

## Filter count matrix
filter <- apply(counts, 1, function(x) mean(x) > 5)
counts <- counts[filter, ]
dim(counts)
## [1] 26742    11
## Save for gene, exon and junction comparisons
counts_gene <- counts
counts_gene[1:5, 1:5]
##                    SRR1027171 SRR1027173 SRR1027174 SRR1027175 SRR1027176
## ENSG00000000003.14        293        242        321        513        181
## ENSG00000000005.5           6          6          3          2        119
## ENSG00000000419.12        583        303        336        546        391
## ENSG00000000457.13        418        334        197        343        290
## ENSG00000000460.16        381        205        151        209        264

To get a better sense of the data, we plot the mean-variance relationship for each gene. Similarly, we run principal component analysis (PCA) to identify any sample outliers within the data. We assess the variance explained by each of the first 11 PCs as well as visualize the relationship of each sample in the first two PCs.

## Set colors 
trop <- RSkittleBrewer('tropical')[c(1, 2)]
cols <- as.numeric(as.factor(rse$group))

## Look at mean variance relationship
plot(rowMeans(log2(counts + 1)), rowVars(log2(counts + 1)),
     pch = 19, col = trop[2])

## Calculate PCs with svd function
expr.pca <- svd(counts - rowMeans(counts))

## Plot PCs
par(font.lab = 2, cex.lab = 1.2, font.axis = 2, cex.axis = 1.2)
plot(expr.pca$d^2 / sum(expr.pca$d^2), pch = 19, col = trop[2], cex = 1.5,
     ylab = 'Fraction of variance explained (gene level)', xlab = 'PC #',
     main = 'PCs')

## Plot PC1 vs. PC2
par(font.lab = 2, cex.lab = 1.2, font.axis = 2, cex.axis = 1.2)
plot(expr.pca$v[, 1], expr.pca$v[, 2], pch = 19, col = trop[cols], cex = 1.5,
     xlab = 'PC1', ylab = 'PC2',
     main = 'PC (gene level)')
legend('topright', pch = 19, col = trop[c(1, 2)],
       names(summary(as.factor(rse$group))), bg="white")

Having determined there are no sample outliers in these data, we carry out differential gene expression analysis. Differential gene expression between TNBC and HER2-positive samples are determined using limma and voom. Differentially expressed genes are visualized using a volcano plot to compare the effect size of the differential expression [ as measured by the \(log_2(fold-change)\) in expression ] and its significance [ \(-log_10(p-value)\) ].

## Perform differential expression analysis with limma-voom
design <- model.matrix(~ rse$group)
design
##    (Intercept) rse$groupTNBC Breast Tumor
## 1            1                          1
## 2            1                          1
## 3            1                          1
## 4            1                          1
## 5            1                          1
## 6            1                          0
## 7            1                          0
## 8            1                          0
## 9            1                          0
## 10           1                          0
## 11           1                          1
## attr(,"assign")
## [1] 0 1
## attr(,"contrasts")
## attr(,"contrasts")$`rse$group`
## [1] "contr.treatment"
dge <- DGEList(counts = counts)
dge <- calcNormFactors(dge)
v <- voom(dge, design, plot = TRUE)

fit <- lmFit(v, design)
fit <- eBayes(fit)
log2FC <- fit$coefficients[, 2]
p.mod <- fit$p.value[, 2]
q.mod <- qvalue(p.mod)$q
res_gene <- data.frame(log2FC, p.mod, q.mod)
rownames(res_gene) <- rownames(counts)

## Determine the number of genes differentially expressed at q<0.05
sum(res_gene$q.mod < 0.05)
## [1] 2974
table(res_gene$log2FC[res_gene$q.mod < 0.05] > 0 )
## 
## FALSE  TRUE 
##  1612  1362
## Histogram of p-values
par(font.lab = 2, cex.lab = 1.2, font.axis = 2, cex.axis = 1.2)
hist(p.mod, col = trop[2], xlab = 'p-value',
     main = 'Histogramm of p-values', breaks = 100)

## Volcano plot
par(font.lab = 2, cex.lab = 1.2, font.axis = 2, cex.axis = 1.2)
rx2 <- c(-1, 1) * 1.1 * max(abs(log2FC))
ry2 <- c(-0.1, max(-log10(p.mod))) * 1.1
plot(log2FC, -log10(p.mod), 
     pch = 19, xlim = rx2, ylim = ry2, col = trop[2],
     xlab = bquote(paste(log[2], ' (fold change)')), 
     ylab = bquote(paste(-log[10], ' (p-value)')))
abline(v = seq(-10, 10, 1), col = 'lightgray', lty = 'dotted')
abline(h = seq(0, 23, 1), col = 'lightgray', lty = 'dotted')
points(log2FC, -log10(p.mod), pch = 19, col = trop[2])
title('Volcano plot: TNBC vs. HER2+ in SRP032789 (gene level)')

2 Gene set enrichment analysis

To get a better understanding of those genes showing differential gene expression, we utilize topGO, a gene set analysis library. Genes included in this analysis are those reaching a q-value cutoff less than 0.05.

names(q.mod) <- rownames(counts)
interesting <- function(x) x < 0.05

After determining which genes to include for analysis, topGO objects are generated and the enrichment tests are run. The Kolomogorov-Smirnov (ks) test is used to test for distributional differences. Here, we ask whether each GO group is “enriched” for differentially expressed (q.mod < 0.05) genes. Equivalently, we are testing whether the p-value distributions are the same for genes in and outside of each gene ontology. We run tests on the “biological processes” ontology.

toens <- function(x) {
    res <- x
    names(res) <- gsub('\\..*', '', names(x))
    return(res)
}

topgoobjBP <- new('topGOdata',
    description = 'biological process',
    ontology = 'BP', allGenes = toens(q.mod), geneSelectionFun = interesting,
    annotationFun = annFUN.org, mapping = 'org.Hs.eg.db', ID = 'ensembl')
## 
## Building most specific GOs .....
## Loading required package: org.Hs.eg.db
## 
##  ( 10869 GO terms found. )
## 
## Build GO DAG topology ..........
##  ( 14712 GO terms and 34861 relations. )
## 
## Annotating nodes ...............
##  ( 14653 genes annotated to the GO terms. )
bptest <- runTest(topgoobjBP, algorithm = 'weight01', statistic = 'ks')
## 
##           -- Weight01 Algorithm -- 
## 
##       the algorithm is scoring 14712 nontrivial nodes
##       parameters: 
##           test statistic: ks
##           score order: increasing
## 
##   Level 20:  1 nodes to be scored    (0 eliminated genes)
## 
##   Level 19:  8 nodes to be scored    (0 eliminated genes)
## 
##   Level 18:  18 nodes to be scored   (1 eliminated genes)
## 
##   Level 17:  44 nodes to be scored   (30 eliminated genes)
## 
##   Level 16:  119 nodes to be scored  (84 eliminated genes)
## 
##   Level 15:  242 nodes to be scored  (178 eliminated genes)
## 
##   Level 14:  481 nodes to be scored  (528 eliminated genes)
## 
##   Level 13:  834 nodes to be scored  (1215 eliminated genes)
## 
##   Level 12:  1212 nodes to be scored (2372 eliminated genes)
## 
##   Level 11:  1559 nodes to be scored (4473 eliminated genes)
## 
##   Level 10:  1947 nodes to be scored (6232 eliminated genes)
## 
##   Level 9:   2057 nodes to be scored (8549 eliminated genes)
## 
##   Level 8:   1949 nodes to be scored (10328 eliminated genes)
## 
##   Level 7:   1794 nodes to be scored (11631 eliminated genes)
## 
##   Level 6:   1315 nodes to be scored (12606 eliminated genes)
## 
##   Level 5:   729 nodes to be scored  (13331 eliminated genes)
## 
##   Level 4:   304 nodes to be scored  (13908 eliminated genes)
## 
##   Level 3:   77 nodes to be scored   (14160 eliminated genes)
## 
##   Level 2:   21 nodes to be scored   (14343 eliminated genes)
## 
##   Level 1:   1 nodes to be scored    (14427 eliminated genes)
bptest
## 
## Description: biological process 
## Ontology: BP 
## 'weight01' algorithm with the 'ks' test
## 14712 GO terms scored: 53 terms with p < 0.01
## Annotation data:
##     Annotated genes: 14653 
##     Significant genes: 1457 
##     Min. no. of genes annotated to a GO: 1 
##     Nontrivial nodes: 14712
bpres_gene <- GenTable(topgoobjBP, pval = bptest,
                       topNodes = length(bptest@score), numChar = 100)
head(bpres_gene, n = 10)
##         GO.ID                                      Term Annotated
## 1  GO:0016579                  protein deubiquitination       115
## 2  GO:0016569                    chromatin modification       509
## 3  GO:0030049                   muscle filament sliding        32
## 4  GO:0045494            photoreceptor cell maintenance        32
## 5  GO:0033962 cytoplasmic mRNA processing body assembly        20
## 6  GO:0050776             regulation of immune response       762
## 7  GO:0000042                protein targeting to Golgi        18
## 8  GO:0030521       androgen receptor signaling pathway        60
## 9  GO:0007050                         cell cycle arrest       229
## 10 GO:0071557              histone H3-K27 demethylation         4
##    Significant Expected    pval
## 1           18    11.43   1e-05
## 2           70    50.61 0.00019
## 3            6     3.18 0.00023
## 4            4     3.18 0.00036
## 5            4     1.99 0.00061
## 6           49    75.77 0.00097
## 7            6     1.79 0.00100
## 8           14     5.97 0.00110
## 9           32    22.77 0.00118
## 10           4     0.40 0.00144

3 Exon level analysis

As above, we are interested here in differential expression. However, rather than summarizing across genes, this analysis will look for differential expression at the exon level. In this analysis, we include all exons that map to the previous filtered genes and again carry out differential expression analysis using limma and voom.

Here, we download data from the same project as above (SRP032798); however, this time, we are interested in obtaining the exon level data.

## Find a project of interest (SRP032789)
project_info <- abstract_search('To define the digital transcriptome of three breast cancer')
project_info
##     number_samples species
## 865             20   human
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               abstract
## 865 Goal: To define the digital transcriptome of three breast cancer subtypes (TNBC, Non-TNBC, and HER2-positive) using RNA-sequencing technology. To elucidate differentially expressed known and novel transcripts, alternatively spliced genes and differential isoforms and lastly expressed variants in our dataset. Method: Dr. Suzanne Fuqua (Baylor College of Medicine) provided the human breast cancer tissue RNA samples. All of the human samples were used in accordance with the IRB procedures of Baylor College of Medicine. The breast tumour types, TNBC, Non-TNBC and HER2-positive, were classified on the basis of immunohistochemical and RT-qPCR classification. Results: Comparative transcriptomic analyses elucidated differentially expressed transcripts between the three breast cancer groups, identifying several new modulators of breast cancer. We discovered subtype specific differentially spliced genes and splice isoforms not previously recognized in human transcriptome. Further, we showed that exon skip and intron retention are predominant splice events in breast cancer. In addition, we found that differential expression of primary transcripts and promoter switching are significantly deregulated in breast cancer compared to normal breast.  We also report novel expressed variants, allelic prevalence and abundance, and coexpression with other variation, and splicing signatures. Additionally we describe novel SNPs and INDELs in cancer relevant genes with no prior reported association of point mutations with cancer Overall design: mRNA profiles of 17 breast tumor samples of three different subtypes (TNBC, non-TNBC and HER2-positive) and normal human breast organoids (epithelium) samples (NBS) were sequenced using Illumina HiSeq.
##       project
## 865 SRP032789
## Browse the project at SRA
browse_study(project_info$project)

## Download the exon level RangedSummarizedExperiment data
if(!file.exists(file.path('SRP032789', 'rse_exon.Rdata'))) {
    download_study(project_info$project, type = 'rse-exon')
}

## Load the data
load(file.path(project_info$project, 'rse_exon.Rdata'))
rse_exon
## class: RangedSummarizedExperiment 
## dim: 329092 20 
## metadata(0):
## assays(1): counts
## rownames(329092): ENSG00000000003.14 ENSG00000000003.14 ...
##   ENSG00000283698.1 ENSG00000283699.1
## rowData names(0):
## colnames(20): SRR1027171 SRR1027173 ... SRR1027190 SRR1027172
## colData names(21): project sample ... title characteristics
## This is the sample phenotype data provided by the recount project
colData(rse_exon)
## DataFrame with 20 rows and 21 columns
##                project      sample  experiment         run
##            <character> <character> <character> <character>
## SRR1027171   SRP032789   SRS500214   SRX374850  SRR1027171
## SRR1027173   SRP032789   SRS500216   SRX374852  SRR1027173
## SRR1027174   SRP032789   SRS500217   SRX374853  SRR1027174
## SRR1027175   SRP032789   SRS500218   SRX374854  SRR1027175
## SRR1027176   SRP032789   SRS500219   SRX374855  SRR1027176
## ...                ...         ...         ...         ...
## SRR1027187   SRP032789   SRS500230   SRX374866  SRR1027187
## SRR1027188   SRP032789   SRS500231   SRX374867  SRR1027188
## SRR1027189   SRP032789   SRS500232   SRX374868  SRR1027189
## SRR1027190   SRP032789   SRS500233   SRX374869  SRR1027190
## SRR1027172   SRP032789   SRS500215   SRX374851  SRR1027172
##            read_count_as_reported_by_sra reads_downloaded
##                                <integer>        <integer>
## SRR1027171                      88869444         88869444
## SRR1027173                     107812596        107812596
## SRR1027174                      98563260         98563260
## SRR1027175                      91327892         91327892
## SRR1027176                      96513572         96513572
## ...                                  ...              ...
## SRR1027187                      75260678         75260678
## SRR1027188                      65709192         65709192
## SRR1027189                      65801392         65801392
## SRR1027190                      74356276         74356276
## SRR1027172                      80986440         58902122
##            proportion_of_reads_reported_by_sra_downloaded paired_end
##                                                 <numeric>  <logical>
## SRR1027171                                              1       TRUE
## SRR1027173                                              1       TRUE
## SRR1027174                                              1       TRUE
## SRR1027175                                              1       TRUE
## SRR1027176                                              1       TRUE
## ...                                                   ...        ...
## SRR1027187                                      1.0000000       TRUE
## SRR1027188                                      1.0000000       TRUE
## SRR1027189                                      1.0000000       TRUE
## SRR1027190                                      1.0000000       TRUE
## SRR1027172                                      0.7273084       TRUE
##            sra_misreported_paired_end mapped_read_count        auc
##                             <logical>         <integer>  <numeric>
## SRR1027171                      FALSE          86949307 5082692127
## SRR1027173                      FALSE         104337779 6077034329
## SRR1027174                      FALSE          95271238 5504462845
## SRR1027175                      FALSE          88820239 5150234117
## SRR1027176                      FALSE          93464650 5416681912
## ...                               ...               ...        ...
## SRR1027187                      FALSE          64697612 3567078255
## SRR1027188                      FALSE          65278500 4856453823
## SRR1027189                      FALSE          65328289 4858587600
## SRR1027190                      FALSE          73911898 5501089036
## SRR1027172                      FALSE          57523391 3351013968
##            sharq_beta_tissue sharq_beta_cell_type
##                  <character>          <character>
## SRR1027171            breast                  esc
## SRR1027173            breast                  esc
## SRR1027174            breast                  esc
## SRR1027175            breast                  esc
## SRR1027176            breast                  esc
## ...                      ...                  ...
## SRR1027187            breast                  esc
## SRR1027188            breast                  esc
## SRR1027189            breast                  esc
## SRR1027190            breast                  esc
## SRR1027172            breast                  esc
##            biosample_submission_date biosample_publication_date
##                          <character>                <character>
## SRR1027171   2013-11-07T12:40:22.203    2013-11-08T01:11:17.160
## SRR1027173   2013-11-07T12:40:32.283    2013-11-08T01:11:14.827
## SRR1027174   2013-11-07T12:40:28.283    2013-11-08T01:11:52.283
## SRR1027175   2013-11-07T12:40:34.343    2013-11-08T01:11:15.963
## SRR1027176   2013-11-07T12:40:36.303    2013-11-08T01:11:46.430
## ...                              ...                        ...
## SRR1027187   2013-11-07T12:40:56.180    2013-11-08T01:11:29.587
## SRR1027188   2013-11-07T12:40:58.170    2013-11-08T01:12:06.660
## SRR1027189   2013-11-07T12:40:20.227    2013-11-08T01:11:33.080
## SRR1027190   2013-11-07T12:40:18.090    2013-11-08T01:12:11.320
## SRR1027172   2013-11-07T12:40:26.217    2013-11-08T01:11:45.250
##              biosample_update_date avg_read_length geo_accession
##                        <character>       <integer>   <character>
## SRR1027171 2014-03-07T16:09:38.542             120    GSM1261016
## SRR1027173 2014-03-07T16:09:38.698             120    GSM1261018
## SRR1027174 2014-03-07T16:09:38.637             120    GSM1261019
## SRR1027175 2014-03-07T16:09:38.731             120    GSM1261020
## SRR1027176 2014-03-07T16:09:38.768             120    GSM1261021
## ...                            ...             ...           ...
## SRR1027187 2014-03-07T16:09:39.093             120    GSM1261032
## SRR1027188 2014-03-07T16:09:39.130             150    GSM1261033
## SRR1027189 2014-03-07T16:09:38.498             150    GSM1261034
## SRR1027190 2014-03-07T16:09:38.469             150    GSM1261035
## SRR1027172 2014-03-07T16:09:38.604              87    GSM1261017
##              bigwig_file       title
##              <character> <character>
## SRR1027171 SRR1027171.bw       TNBC1
## SRR1027173 SRR1027173.bw       TNBC3
## SRR1027174 SRR1027174.bw       TNBC4
## SRR1027175 SRR1027175.bw       TNBC5
## SRR1027176 SRR1027176.bw       TNBC6
## ...                  ...         ...
## SRR1027187 SRR1027187.bw      HER2-5
## SRR1027188 SRR1027188.bw        NBS1
## SRR1027189 SRR1027189.bw        NBS2
## SRR1027190 SRR1027190.bw        NBS3
## SRR1027172 SRR1027172.bw       TNBC2
##                                   characteristics
##                                   <CharacterList>
## SRR1027171          tumor type: TNBC Breast Tumor
## SRR1027173          tumor type: TNBC Breast Tumor
## SRR1027174          tumor type: TNBC Breast Tumor
## SRR1027175          tumor type: TNBC Breast Tumor
## SRR1027176          tumor type: TNBC Breast Tumor
## ...                                           ...
## SRR1027187 tumor type: HER2 Positive Breast Tumor
## SRR1027188    tumor type: Normal Breast Organoids
## SRR1027189    tumor type: Normal Breast Organoids
## SRR1027190    tumor type: Normal Breast Organoids
## SRR1027172          tumor type: TNBC Breast Tumor

As above, downloaded count data are first scaled to take into account differing coverage between samples. The same phenotype data (pheno) are used and again ordered to match the sample order of the expression data (rse_exon). Only those samples that are HER2-positive or TNBC are included for analysis. Prior to differential exon expression analysis, count data are obtained in matrix format and then filtered to only include exons within genes that had been analyzed previously.

## Scale counts by taking into account the total coverage per sample
rse <- scale_counts(rse_exon)

## Download pheno data from 
## http://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP032789
pheno <- read.table('SraRunTable_SRP032789.txt', sep = '\t', 
                    header = TRUE,
                    stringsAsFactors = FALSE)

## Obtain correct order for pheno data
pheno <- pheno[match(rse$run, pheno$Run_s), ]
identical(pheno$Run_s, rse$run)
## [1] TRUE
head(cbind(pheno$Run_s, rse$run))
##      [,1]         [,2]        
## [1,] "SRR1027171" "SRR1027171"
## [2,] "SRR1027173" "SRR1027173"
## [3,] "SRR1027174" "SRR1027174"
## [4,] "SRR1027175" "SRR1027175"
## [5,] "SRR1027176" "SRR1027176"
## [6,] "SRR1027177" "SRR1027177"
## Obtain grouping information
colData(rse)$group <- pheno$tumor_type_s
table(colData(rse)$group)   
## 
## HER2 Positive Breast Tumor      Non-TNBC Breast Tumor 
##                          5                          6 
##    Normal Breast Organoids          TNBC Breast Tumor 
##                          3                          6
## Subset data to HER2 and TNBC types
rse <- rse[, rse$group %in% c('HER2 Positive Breast Tumor',
                              'TNBC Breast Tumor')]

## Save filtered rse object
rse_exon_filt <- rse
rse_exon_filt
## class: RangedSummarizedExperiment 
## dim: 329092 11 
## metadata(0):
## assays(1): counts
## rownames(329092): ENSG00000000003.14 ENSG00000000003.14 ...
##   ENSG00000283698.1 ENSG00000283699.1
## rowData names(0):
## colnames(11): SRR1027171 SRR1027173 ... SRR1027187 SRR1027172
## colData names(22): project sample ... characteristics group
## Obtain count matrix
counts <- assays(rse_exon_filt)$counts
dim(counts)
## [1] 329092     11
## Filter count matrix (keep exons that are in filtered gene counts matrix)
filter <- rownames(counts) %in% rownames(counts_gene)
counts <- counts[filter, ]
dim(counts)
## [1] 259626     11
## Save for gene, exon and junction comparisons
counts_exon <- counts
counts_exon[1:5, 1:5]
##                    SRR1027171 SRR1027173 SRR1027174 SRR1027175 SRR1027176
## ENSG00000000003.14        151        135        169        252         96
## ENSG00000000003.14         25         20         18         30         10
## ENSG00000000003.14          0          0          0          0          0
## ENSG00000000003.14         15         14         14         26          8
## ENSG00000000003.14         22         15         27         30         13

As above, to get a better sense of the data, we assess the mean-variance relationship for each exon. Similarly, we run principal component analysis (PCA) to identify any sample outliers within the data. We assess the variance explained by each of the first 11 PCs as well as visualize the relationship of each sample in the first two PCs.

## Set colors 
trop <- RSkittleBrewer('tropical')[c(1, 2)]
cols <- as.numeric(as.factor(rse$group))

## Look at mean variance relationship
plot(rowMeans(log2(counts + 1)), rowVars(log2(counts + 1)),
     pch = 19, col = trop[2])

## Calculate PCs with svd function
expr.pca <- svd(counts - rowMeans(counts))

## Plot PCs
par(font.lab = 2, cex.lab = 1.2, font.axis = 2, cex.axis = 1.2)
plot(expr.pca$d^2 / sum(expr.pca$d^2), pch = 19, col = trop[2], cex = 1.5,
     ylab = 'Fraction of variance explained', xlab = 'PC #',
     main = 'PCs (exon level)')

## Plot PC1 vs. PC2
par(font.lab = 2, cex.lab = 1.2, font.axis = 2, cex.axis = 1.2)
plot(expr.pca$v[, 1], expr.pca$v[, 2], pch = 19, col = trop[cols], cex = 1.5,
     xlab = 'PC1', ylab = 'PC2',
     main = 'PC (exon level)')
legend('topright', pch = 19, col = trop[c(1, 2)],
       names(summary(as.factor(rse$group))), bg="white")

Again, differential expression analysis is carried out using limma and voom; however, this time at the exon, rather than gene, level. Data are again visualized using a volcano plot to assess the strength [ \(log_2(fold-change)\) in expression ] and its significance [ \(-log_10(p-value)\) ].for each exon.

design <- model.matrix(~ rse$group)
design
##    (Intercept) rse$groupTNBC Breast Tumor
## 1            1                          1
## 2            1                          1
## 3            1                          1
## 4            1                          1
## 5            1                          1
## 6            1                          0
## 7            1                          0
## 8            1                          0
## 9            1                          0
## 10           1                          0
## 11           1                          1
## attr(,"assign")
## [1] 0 1
## attr(,"contrasts")
## attr(,"contrasts")$`rse$group`
## [1] "contr.treatment"
dge <- DGEList(counts = counts)
dge <- calcNormFactors(dge)
v <- voom(dge, design,plot = TRUE)