1 Analysis
To illustrate the ease of combining data from multiple projects included in recount as part of a cross-study meta-analysis, we present here a cross-tissue differential expression (DE) analysis comparing gene expression between colon and whole blood. We perform an initial analysis on samples from a variety of studies found in the recount project and then compare the results we obtain to those from analysis of GTEx data in the same tissues.
1.1 Load libraries we will need
library('dplyr')
library('recount')
library('magrittr')
library('limma')
library('edgeR')
library('ffpe')
library('RSkittleBrewer')
library('SummarizedExperiment')
library('devtools')
trop <- RSkittleBrewer::RSkittleBrewer('tropical')1.2 Get data sets
As an initial analysis, colon samples labeled as controls were taken from studies SRP029880 (a study of colorectal cancer (Kim, Kim, Kim, Roh, et al., 2014), n=19) and SRP042228 (a study of Crohn’s disease (Denson, 2014), n=41). Whole blood samples labeled as controls were taken from SRP059039 (a study of virus-caused diarrhea, unpublished, n=24), SRP059172 (a study of blood biomarkers for brucellosis, unpublished, n=47) and SRP062966 (a study of lupus, unpublished, n=18).
We start by downloading .Rdata files containing SummarizedExperiment objects for each of these studies, which contain RNA-seq counts at the gene level for each individual in the study, as well as accompanying sample-specific phenotype data.
colon_proj <- c('SRP029880', 'SRP042228')
if(any(!file.exists(file.path(colon_proj, 'rse_gene.Rdata')))) {
sapply(colon_proj, download_study)
}
blood_proj <- c('SRP059039', 'SRP059172', 'SRP062966')
if(any(!file.exists(file.path(blood_proj, 'rse_gene.Rdata')))) {
sapply(blood_proj, download_study)
}
proj <- c(colon_proj,blood_proj)1.3 Load the data
Once the data has been downloaded, we load all five data sets into R and calculate the number of genes (number of rows) and samples (number of columns) for each data set.
dat <- lapply(proj, function(x) {
load(file.path(x, 'rse_gene.Rdata'))
return(rse_gene)
})
proj## [1] "SRP029880" "SRP042228" "SRP059039" "SRP059172" "SRP062966"sapply(dat, dim)## [,1] [,2] [,3] [,4] [,5]
## [1,] 58037 58037 58037 58037 58037
## [2,] 54 314 205 169 1171.4 Load the metadata file
Since the data contained in the SummarizedExperiment objects for each study are somewhat limited, we also load a data frame containing metadata for each sample, including Gene Expression Omnibus (GEO) sample identifiers. We will use this in order to identify and extract the control subjects from each study for analysis.
metadata <- all_metadata('sra')## 2017-05-15 13:10:24 downloading the metadata to /var/folders/cx/n9s558kx6fb7jf5z_pgszgb80000gn/T//RtmpIlquGX/metadata_clean_sra.Rdata1.5 Get additional GEI data for these samples
Now we go through and collect GEO information for the samples. We label them with their respective tissue and identify which samples are supposed to be normal controls, with labels varying from data set to data set. We also label the data sets with their tissue of origin.
if(!file.exists('charvec.Rdata')) {
charvec <- vector('list', 5)
dir.create('geoinfo', showWarnings = FALSE)
for(i in 1:5){
index <- match(colData(dat[[i]])$run, metadata$run)
colData(dat[[i]])$geo <- metadata$geo_accession[index]
info <- sapply(colData(dat[[i]])$geo, geo_info, destdir = 'geoinfo')
charvec[[i]] <- sapply(info, geo_characteristics)
}
save(charvec, file = 'charvec.Rdata')
} else {
load('charvec.Rdata')
}
## first data set - normals called 'normal-looking surrounding colonic epithelium'
colData(dat[[1]])$normal <- grepl('normal', unlist(charvec[1])[(1:54) * 2 - 1])
colData(dat[[1]])$tissue <- 'colon'
## second data set - normals called 'not ibd'
colData(dat[[2]])$normal <- grepl('not ibd', tolower(unlist(charvec[[2]][5, ])))
colData(dat[[2]])$tissue <- 'colon'
## third data set - normals called 'Control'
colData(dat[[3]])$normal <- grepl('Control', unlist(charvec[[3]][2, ]))
colData(dat[[3]])$tissue <- 'blood'
## fourth data set - normals called 'Control'
colData(dat[[4]])$normal <- grepl('Control', unlist(charvec[[4]][1, ]))
colData(dat[[4]])$tissue <- 'blood'
## fifth data set - normals called 'healthy'
colData(dat[[5]])$normal <- grepl('healthy', unlist(charvec[[5]][1, ]))
colData(dat[[5]])$tissue <- 'blood'1.6 Merge the data sets
We merge the data sets into one ranged SummarizedExperiment.
mdat <- do.call(cbind, dat)1.7 Subset to normal control individuals
Using the label determined above, we find out how many samples are normal in each study and subset to just the normal samples for further analysis.
table(colData(mdat)$normal, colData(mdat)$project)##
## SRP029880 SRP042228 SRP059039 SRP059172 SRP062966
## FALSE 35 273 181 122 99
## TRUE 19 41 24 47 18ndat <- mdat[, colData(mdat)$normal]1.8 Analysis comparing blood to colon
Here we do a differential expression analysis comparing blood to colon using voom (Law, Chen, Shi, and Smyth, 2014). We start by normalizing the counts and then consider only the genes that have an average normalized count of at least 5 across the data set.
ndat <- scale_counts(ndat)
ndat_counts <- assays(ndat)[[1]]
keep <- rowMeans(ndat_counts) > 5
ndat_counts = ndat_counts[keep, ]
design <- model.matrix(~colData(ndat)$tissue)
dge <- DGEList(counts = ndat_counts)
dge <- calcNormFactors(dge)
v <- voom(dge, design, plot=FALSE)
fit <- lmFit(v, design)
fit <- eBayes(fit)
topTable(fit)## Removing intercept from test coefficients## logFC AveExpr t P.Value
## ENSG00000164318.17 7.529724 -3.06463060 80.36671 2.831746e-126
## ENSG00000006611.15 13.371932 -0.75728329 75.76053 1.768905e-122
## ENSG00000127324.8 15.539959 0.08913163 75.61819 2.335941e-122
## ENSG00000126785.12 8.823754 -2.29940093 74.33356 2.938672e-121
## ENSG00000237187.8 7.355825 -3.08459539 74.51823 2.036866e-121
## ENSG00000101825.7 11.228958 -1.57276438 71.87321 4.221307e-119
## ENSG00000151617.15 7.904513 -2.91782300 70.75199 4.283362e-118
## ENSG00000173702.7 15.051876 0.25301464 70.31748 1.061250e-117
## ENSG00000118526.6 10.537333 -1.83523586 68.77996 2.746173e-116
## ENSG00000183580.9 8.579915 -2.49579830 68.87721 2.230948e-116
## adj.P.Val B
## ENSG00000164318.17 7.849883e-122 272.6540
## ENSG00000006611.15 2.158488e-118 265.7404
## ENSG00000127324.8 2.158488e-118 265.4671
## ENSG00000126785.12 1.629258e-117 262.8825
## ENSG00000237187.8 1.411599e-117 262.4106
## ENSG00000101825.7 1.950314e-115 258.4736
## ENSG00000151617.15 1.696273e-114 255.6127
## ENSG00000173702.7 3.677364e-114 255.5598
## ENSG00000118526.6 7.612666e-113 252.3701
## ENSG00000183580.9 6.871567e-113 252.28461.8.1 GTEx analysis: comparing same tissues
In order to have something to compare to the results generated on the data sets selected above, we do a differential expression analysis comparing blood to colon samples from GTEx, project SRP012682. We carry out the same steps as above for the analysis, but in this case we also control for batch as included in the GTEx metadata.
## Download the GTEx data
if(!file.exists(file.path('SRP012682', 'rse_gene.Rdata'))) {
download_study('SRP012682')
}
load(file.path('SRP012682', 'rse_gene.Rdata'))
## Download GTEx metadata
gtex_metadata <- all_metadata('gtex')## 2017-05-15 13:11:18 downloading the metadata to /var/folders/cx/n9s558kx6fb7jf5z_pgszgb80000gn/T//RtmpIlquGX/metadata_clean_gtex.Rdata## Subset to samples from tissues of interest
gtex_blood <- rse_gene[, subset(gtex_metadata, smtsd == 'Whole Blood')$run]
colData(gtex_blood)$tissue <- 'wholeblood'
gtex_colon <- rse_gene[, subset(gtex_metadata, smts == 'Colon')$run]
colData(gtex_colon)$tissue <- 'colon'
## Combine tissue subsets and include batch label
gtex_both <- do.call(cbind, list(gtex_blood, gtex_colon))
colData(gtex_both)$batch <- gtex_metadata[match(colData(gtex_both)$run,
gtex_metadata$run), 'smgebtch']
## Scale and filter data to the same genes kept in initial analysis
gtex_both <- scale_counts(gtex_both)
gtex_both_counts <- assays(gtex_both)[[1]]
gtex_both_counts <- gtex_both_counts[keep, ]
## Carry out differential expression analysis
design_gtex <- model.matrix(~colData(gtex_both)$tissue +
colData(gtex_both)$batch)
dge_gtex <- DGEList(counts = gtex_both_counts)
dge_gtex <- calcNormFactors(dge_gtex)
v_gtex <- voom(dge_gtex, design_gtex, plot=FALSE)
fit_gtex <- lmFit(v_gtex, design_gtex)
fit_gtex <- eBayes(fit_gtex)
topTable(fit_gtex, coef = 2)## logFC AveExpr t P.Value adj.P.Val B
## ENSG00000072786.12 4.595693 6.484796 87.13601 0 0 852.9246
## ENSG00000150760.12 -5.426947 2.636342 -85.76042 0 0 842.4991
## ENSG00000104894.11 5.428923 6.456112 85.31817 0 0 839.5314
## ENSG00000198752.10 -3.524774 4.592607 -84.15769 0 0 830.8078
## ENSG00000141480.17 5.798289 7.545908 83.97986 0 0 829.5412
## ENSG00000083857.13 -7.000981 2.562796 -83.19560 0 0 823.2757
## ENSG00000136167.13 7.773033 8.014350 83.15074 0 0 823.2608
## ENSG00000142347.16 7.257857 6.974144 81.95287 0 0 814.1077
## ENSG00000116701.14 8.176832 6.252451 80.64959 0 0 803.9663
## ENSG00000151651.15 7.258016 5.728969 80.47561 0 0 802.62141.8.2 GTEx analysis: comparing different tissues
To have a way of gauguing how much overlap we expect to see between the results from the five combined data sets and those from GTEx data, we also carry out a differential analysis using GTEx data comparing blood to lung. We expect to see much less concordance between our initial results and the GTEx blood-lung comparison than we see between our initial results and the GTEx blood-colon comparison.
## Subset lung data
gtex_lung <- rse_gene[, subset(gtex_metadata, smts=='Lung')$run]
colData(gtex_lung)$tissue <- 'lung'
## Combine lung and blood data and label with batch
gtex_both_lung <- do.call(cbind, list(gtex_blood, gtex_lung))
colData(gtex_both_lung)$batch <- gtex_metadata[
match(
colData(gtex_both_lung)$run,
gtex_metadata$run
), 'smgebtch']
## Scale and filter to correct gene set
gtex_both_lung <- scale_counts(gtex_both_lung)
gtex_both_lung_counts <- assays(gtex_both_lung)[[1]]
gtex_both_lung_counts <- gtex_both_lung_counts[keep,]
## Carry out differential expression analysis
design_gtex_lung <- model.matrix(~colData(gtex_both_lung)$tissue +
colData(gtex_both_lung)$batch)
dge_gtex_lung <- DGEList(counts = gtex_both_lung_counts)
dge_gtex_lung <- calcNormFactors(dge_gtex_lung)
v_gtex_lung <- voom(dge_gtex_lung, design_gtex_lung, plot = FALSE)
fit_gtex_lung <- lmFit(v_gtex_lung, design_gtex_lung)
fit_gtex_lung <- eBayes(fit_gtex_lung)
topTable(fit_gtex_lung, coef = 2)## logFC AveExpr t P.Value adj.P.Val
## ENSG00000069122.18 -7.834493 3.1073837 -109.38749 0 0
## ENSG00000118526.6 -8.943232 1.4057694 -104.47541 0 0
## ENSG00000175899.14 -8.082122 7.0747207 -102.37996 0 0
## ENSG00000129473.9 -2.823777 4.0137367 -100.97719 0 0
## ENSG00000121068.13 -8.283958 2.6788408 -100.14008 0 0
## ENSG00000128567.16 -6.259427 3.7794436 -100.06553 0 0
## ENSG00000112655.15 -5.255477 2.6448203 -99.02082 0 0
## ENSG00000091136.13 -6.245406 4.6398510 -98.42636 0 0
## ENSG00000225383.7 -8.370107 -0.8413745 -96.93016 0 0
## ENSG00000227954.6 -8.626676 0.6834325 -96.88464 0 0
## B
## ENSG00000069122.18 1011.7457
## ENSG00000118526.6 981.1413
## ENSG00000175899.14 968.0088
## ENSG00000129473.9 958.9134
## ENSG00000121068.13 953.1886
## ENSG00000128567.16 952.8221
## ENSG00000112655.15 945.8699
## ENSG00000091136.13 941.9876
## ENSG00000225383.7 931.7232
## ENSG00000227954.6 931.42071.9 Compare concordance at the top (CAT) plots
To look for overlap among the most differentially expressed genes across the three analyses, we construct concordance at the top (CAT) plots. We make three comparisons: our five-study results compared to GTEx blood-colon results (same tissue comparison), our five-study results compared to GTEx blood-lung results (different tissue comparison) and our five-study results compared to a ranked set of coefficients for the batch variable in the GTEx blood-lung results (which should not show any biologically meaningful correspondence).
## Comparison of initial analysis to GTEx, same tissues
cat_sra_gtex <- CATplot(
-rank(fit$coefficients[, 2]),
-rank(-fit_gtex$coefficients[, 2]), maxrank = 1000, ylim = c(0,1), make.plot=FALSE)
## Comparison of initial analysis to GTEx, different tissues
cat_sra_gtex_lung = CATplot(
-rank(fit$coefficients[, 2]),
-rank(-fit_gtex_lung$coefficients[, 2]), maxrank = 1000, ylim = c(0,1), make.plot=FALSE)
## Comparison of initial analysis to GTEx, coefficient of batch in model
cat_sra_gtex_batch = CATplot(
-rank(fit$coefficients[, 2]),
-rank(-fit_gtex_lung$coefficients[, 3]), maxrank = 1000, ylim = c(0,1), make.plot=FALSE)
plot(cat_sra_gtex, type = 'l', col = trop[1], lwd = 3)
lines(cat_sra_gtex_lung, type = 'l', col = trop[2], lwd = 3)
lines(cat_sra_gtex_batch, type = 'l', col = trop[3], lwd = 3)
legend(0, 0.5, legend=c('Same Tissue', 'Different Tissues',
'Tissue vs. Batch'), col = trop[1:3], lwd = 3)
From the figure above, we can see that the best concordance is for the comparison to GTEx where we analyzed the same pair of tissues. The GTEx analysis with lung looks considerably less overlapping, with the overlap with the batch coefficient well below the other two.