Here we provide the code which was used to contruct the RangedSummarizedExperiment object of the airway experiment data package. The experiment citation is:
Himes BE, Jiang X, Wagner P, Hu R, Wang Q, Klanderman B, Whitaker RM, Duan Q, Lasky-Su J, Nikolos C, Jester W, Johnson M, Panettieri R Jr, Tantisira KG, Weiss ST, Lu Q. “RNA-Seq Transcriptome Profiling Identifies CRISPLD2 as a Glucocorticoid Responsive Gene that Modulates Cytokine Function in Airway Smooth Muscle Cells.” PLoS One. 2014 Jun 13;9(6):e99625. PMID: 24926665. GEO: GSE52778.
From the abstract, a brief description of the RNA-Seq experiment on airway smooth muscle (ASM) cell lines: “Using RNA-Seq, a high-throughput sequencing method, we characterized transcriptomic changes in four primary human ASM cell lines that were treated with dexamethasone - a potent synthetic glucocorticoid (1 micromolar for 18 hours).”
The following code chunk obtains the sample information from the series matrix file downloaded from GEO. The columns are then parsed and new columns with shorter names and factor levels are added.
suppressPackageStartupMessages( library( "GEOquery" ) )
suppressPackageStartupMessages( library( "airway" ) )
dir <- system.file("extdata",package="airway")
geofile <- file.path(dir, "GSE52778_series_matrix.txt")
gse <- getGEO(filename=geofile)
## Parsed with column specification:
## cols(
## ID_REF = col_character(),
## GSM1275862 = col_character(),
## GSM1275863 = col_character(),
## GSM1275864 = col_character(),
## GSM1275865 = col_character(),
## GSM1275866 = col_character(),
## GSM1275867 = col_character(),
## GSM1275868 = col_character(),
## GSM1275869 = col_character(),
## GSM1275870 = col_character(),
## GSM1275871 = col_character(),
## GSM1275872 = col_character(),
## GSM1275873 = col_character(),
## GSM1275874 = col_character(),
## GSM1275875 = col_character(),
## GSM1275876 = col_character(),
## GSM1275877 = col_character()
## )
## File stored at:
## /tmp/Rtmp97ML1W/GPL11154.soft
pdata <- pData(gse)[,grepl("ch1",names(pData(gse)))]
names(pdata) <- c("treatment","tissue","ercc_mix","cell","celltype")
pdataclean <- data.frame(treatment=sub("treatment: (.*)","\\1",pdata$treatment),
cell=sub("cell line: (.*)","\\1",pdata$cell),
row.names=rownames(pdata))
pdataclean$dex <- ifelse(grepl("Dex",pdataclean$treatment),"trt","untrt")
pdataclean$albut <- ifelse(grepl("Albut",pdataclean$treatment),"trt","untrt")
pdataclean$SampleName <- rownames(pdataclean)
pdataclean$treatment <- NULL
The information which connects the sample information from GEO with the SRA run id is downloaded from SRA using the Send to: File button.
srafile <- file.path(dir, "SraRunInfo_SRP033351.csv")
srp <- read.csv(srafile)
srpsmall <- srp[,c("Run","avgLength","Experiment","Sample","BioSample","SampleName")]
These two data.frames are merged and then we subset to only the samples not treated with albuterol (these samples were not included in the analysis of the publication).
coldata <- merge(pdataclean, srpsmall, by="SampleName")
rownames(coldata) <- coldata$Run
coldata <- coldata[coldata$albut == "untrt",]
coldata$albut <- NULL
coldata
## SampleName cell dex
## SRR1039508 GSM1275862 tissue: human airway smooth muscle cells untrt
## SRR1039509 GSM1275863 tissue: human airway smooth muscle cells untrt
## SRR1039510 GSM1275864 tissue: human airway smooth muscle cells untrt
## SRR1039511 GSM1275865 tissue: human airway smooth muscle cells untrt
## SRR1039512 GSM1275866 tissue: human airway smooth muscle cells untrt
## SRR1039513 GSM1275867 tissue: human airway smooth muscle cells untrt
## SRR1039514 GSM1275868 tissue: human airway smooth muscle cells untrt
## SRR1039515 GSM1275869 tissue: human airway smooth muscle cells untrt
## SRR1039516 GSM1275870 tissue: human airway smooth muscle cells untrt
## SRR1039517 GSM1275871 tissue: human airway smooth muscle cells untrt
## SRR1039518 GSM1275872 tissue: human airway smooth muscle cells untrt
## SRR1039519 GSM1275873 tissue: human airway smooth muscle cells untrt
## SRR1039520 GSM1275874 tissue: human airway smooth muscle cells untrt
## SRR1039521 GSM1275875 tissue: human airway smooth muscle cells untrt
## SRR1039522 GSM1275876 tissue: human airway smooth muscle cells untrt
## SRR1039523 GSM1275877 tissue: human airway smooth muscle cells untrt
## Run avgLength Experiment Sample BioSample
## SRR1039508 SRR1039508 126 SRX384345 SRS508568 SAMN02422669
## SRR1039509 SRR1039509 126 SRX384346 SRS508567 SAMN02422675
## SRR1039510 SRR1039510 126 SRX384347 SRS508570 SAMN02422668
## SRR1039511 SRR1039511 126 SRX384348 SRS508569 SAMN02422667
## SRR1039512 SRR1039512 126 SRX384349 SRS508571 SAMN02422678
## SRR1039513 SRR1039513 87 SRX384350 SRS508572 SAMN02422670
## SRR1039514 SRR1039514 126 SRX384351 SRS508574 SAMN02422681
## SRR1039515 SRR1039515 114 SRX384352 SRS508573 SAMN02422671
## SRR1039516 SRR1039516 120 SRX384353 SRS508575 SAMN02422682
## SRR1039517 SRR1039517 126 SRX384354 SRS508576 SAMN02422673
## SRR1039518 SRR1039518 126 SRX384355 SRS508578 SAMN02422679
## SRR1039519 SRR1039519 107 SRX384356 SRS508577 SAMN02422672
## SRR1039520 SRR1039520 101 SRX384357 SRS508579 SAMN02422683
## SRR1039521 SRR1039521 98 SRX384358 SRS508580 SAMN02422677
## SRR1039522 SRR1039522 125 SRX384359 SRS508582 SAMN02422680
## SRR1039523 SRR1039523 126 SRX384360 SRS508581 SAMN02422674
Finally, the sample table was saved to a CSV file for future
reference. This file is included in the inst/extdata directory of
this package.
write.csv(coldata, file="sample_table.csv")
A file containing the SRA run numbers was created: files. This
file was used to download the sequenced reads from the SRA using
wget. The following command was used to extract the FASTQ file from
the .sra files, using the
SRA Toolkit
cat files | parallel -j 7 fastq-dump --split-files {}.sra
The reads were aligned using the STAR read aligner to GRCh37 using the annotations from Ensembl release 75.
for f in `cat files`; do STAR --genomeDir ../STAR/ENSEMBL.homo_sapiens.release-75 \
--readFilesIn fastq/$f\_1.fastq fastq/$f\_2.fastq \
--runThreadN 12 --outFileNamePrefix aligned/$f.; done
SAMtools was used to generate BAM files.
cat files | parallel -j 7 samtools view -bS aligned/{}.Aligned.out.sam -o aligned/{}.bam
A transcript database for the homo sapiens Ensembl genes was obtained from Biomart.
library( "GenomicFeatures" )
txdb <- makeTranscriptDbFromBiomart( biomart="ensembl", dataset="hsapiens_gene_ensembl")
exonsByGene <- exonsBy( txdb, by="gene" )
The BAM files were specified using the SRR id from the SRA. A yield
size of 2 million reads was used to cap the memory used during
read counting.
sampleTable <- read.csv( "sample_table.csv", row.names=1 )
fls <- file.path("aligned",rownames(sampleTable), ".bam")
library( "Rsamtools" )
bamLst <- BamFileList( fls, yieldSize=2000000 )
The following summarizeOverlaps call distributed the 8 paired-end
BAM files to 8 workers. This used a maximum of 16 Gb per worker and
the time elapsed was 50 minutes.
library( "BiocParallel" )
register( MulticoreParam( workers=8 ) )
library( "GenomicAlignments" )
airway <- summarizeOverlaps( features=exonsByGene, reads=bamLst,
mode="Union", singleEnd=FALSE,
ignore.strand=TRUE, fragments=TRUE )
The sample information was then added as column data.
colData(airway) <- DataFrame( sampleTable )
Finally, we attached the MIAME information using the Pubmed ID.
library( "annotate" )
miame <- list(pmid2MIAME("24926665"))
miame[[1]]@url <- "http://www.ncbi.nlm.nih.gov/pubmed/24926665"
# because R's CHECK doesn't like non-ASCII characters in data objects
# or in vignettes. the actual char was used in the first argument
miame[[1]]@abstract <- gsub("micro","micro",abstract(miame[[1]]))
miame[[1]]@abstract <- gsub("beta","beta",abstract(miame[[1]]))
metadata(airway) <- miame
save(airway, file="airway.RData")
Below we print out some basic summary statistics on the airway
object which is provided by this experiment data package.
library("airway")
data(airway)
airway
## class: RangedSummarizedExperiment
## dim: 64102 8
## metadata(1): ''
## assays(1): counts
## rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99
## rowData names(0):
## colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
## colData names(9): SampleName cell ... Sample BioSample
as.data.frame(colData(airway))
## SampleName cell dex albut Run avgLength Experiment
## SRR1039508 GSM1275862 N61311 untrt untrt SRR1039508 126 SRX384345
## SRR1039509 GSM1275863 N61311 trt untrt SRR1039509 126 SRX384346
## SRR1039512 GSM1275866 N052611 untrt untrt SRR1039512 126 SRX384349
## SRR1039513 GSM1275867 N052611 trt untrt SRR1039513 87 SRX384350
## SRR1039516 GSM1275870 N080611 untrt untrt SRR1039516 120 SRX384353
## SRR1039517 GSM1275871 N080611 trt untrt SRR1039517 126 SRX384354
## SRR1039520 GSM1275874 N061011 untrt untrt SRR1039520 101 SRX384357
## SRR1039521 GSM1275875 N061011 trt untrt SRR1039521 98 SRX384358
## Sample BioSample
## SRR1039508 SRS508568 SAMN02422669
## SRR1039509 SRS508567 SAMN02422675
## SRR1039512 SRS508571 SAMN02422678
## SRR1039513 SRS508572 SAMN02422670
## SRR1039516 SRS508575 SAMN02422682
## SRR1039517 SRS508576 SAMN02422673
## SRR1039520 SRS508579 SAMN02422683
## SRR1039521 SRS508580 SAMN02422677
summary(colSums(assay(airway))/1e6)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 15.16 19.05 20.90 21.94 24.67 30.82
metadata(rowRanges(airway))
## $genomeInfo
## $genomeInfo$`Db type`
## [1] "TranscriptDb"
##
## $genomeInfo$`Supporting package`
## [1] "GenomicFeatures"
##
## $genomeInfo$`Data source`
## [1] "BioMart"
##
## $genomeInfo$Organism
## [1] "Homo sapiens"
##
## $genomeInfo$`Resource URL`
## [1] "www.biomart.org:80"
##
## $genomeInfo$`BioMart database`
## [1] "ensembl"
##
## $genomeInfo$`BioMart database version`
## [1] "ENSEMBL GENES 75 (SANGER UK)"
##
## $genomeInfo$`BioMart dataset`
## [1] "hsapiens_gene_ensembl"
##
## $genomeInfo$`BioMart dataset description`
## [1] "Homo sapiens genes (GRCh37.p13)"
##
## $genomeInfo$`BioMart dataset version`
## [1] "GRCh37.p13"
##
## $genomeInfo$`Full dataset`
## [1] "yes"
##
## $genomeInfo$`miRBase build ID`
## [1] NA
##
## $genomeInfo$transcript_nrow
## [1] "215647"
##
## $genomeInfo$exon_nrow
## [1] "745593"
##
## $genomeInfo$cds_nrow
## [1] "537555"
##
## $genomeInfo$`Db created by`
## [1] "GenomicFeatures package from Bioconductor"
##
## $genomeInfo$`Creation time`
## [1] "2014-07-10 14:55:55 -0400 (Thu, 10 Jul 2014)"
##
## $genomeInfo$`GenomicFeatures version at creation time`
## [1] "1.17.9"
##
## $genomeInfo$`RSQLite version at creation time`
## [1] "0.11.4"
##
## $genomeInfo$DBSCHEMAVERSION
## [1] "1.0"
sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.2 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.9-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.9-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] airway_1.4.0 SummarizedExperiment_1.14.0
## [3] DelayedArray_0.10.0 BiocParallel_1.18.0
## [5] matrixStats_0.54.0 GenomicRanges_1.36.0
## [7] GenomeInfoDb_1.20.0 IRanges_2.18.0
## [9] S4Vectors_0.22.0 GEOquery_2.52.0
## [11] Biobase_2.44.0 BiocGenerics_0.30.0
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.1 pillar_1.3.1 compiler_3.6.0
## [4] XVector_0.24.0 bitops_1.0-6 tools_3.6.0
## [7] zlibbioc_1.30.0 lattice_0.20-38 evaluate_0.13
## [10] tibble_2.1.1 pkgconfig_2.0.2 rlang_0.3.4
## [13] Matrix_1.2-17 xfun_0.6 GenomeInfoDbData_1.2.1
## [16] dplyr_0.8.0.1 stringr_1.4.0 xml2_1.2.0
## [19] knitr_1.22 hms_0.4.2 grid_3.6.0
## [22] tidyselect_0.2.5 glue_1.3.1 R6_2.4.0
## [25] limma_3.40.0 tidyr_0.8.3 readr_1.3.1
## [28] purrr_0.3.2 magrittr_1.5 assertthat_0.2.1
## [31] stringi_1.4.3 RCurl_1.95-4.12 crayon_1.3.4