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)
## Warning: attributes are not identical across measure variables;
## they will be dropped
## 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/RtmprTEP9B/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
## SRR1039508 GSM1275862
## SRR1039509 GSM1275863
## SRR1039510 GSM1275864
## SRR1039511 GSM1275865
## SRR1039512 GSM1275866
## SRR1039513 GSM1275867
## SRR1039514 GSM1275868
## SRR1039515 GSM1275869
## SRR1039516 GSM1275870
## SRR1039517 GSM1275871
## SRR1039518 GSM1275872
## SRR1039519 GSM1275873
## SRR1039520 GSM1275874
## SRR1039521 GSM1275875
## SRR1039522 GSM1275876
## SRR1039523 GSM1275877
## cell
## SRR1039508 Primary ASM cells were isolated from four white aborted lung transplant donors with no chronic illness. Passages 4 to 7 ASM cells maintained in Ham's F12 medium supplemented with 10% FBS were used in all experiments.
## SRR1039509 Primary ASM cells were isolated from four white aborted lung transplant donors with no chronic illness. Passages 4 to 7 ASM cells maintained in Ham's F12 medium supplemented with 10% FBS were used in all experiments.
## SRR1039510 Primary ASM cells were isolated from four white aborted lung transplant donors with no chronic illness. Passages 4 to 7 ASM cells maintained in Ham's F12 medium supplemented with 10% FBS were used in all experiments.
## SRR1039511 Primary ASM cells were isolated from four white aborted lung transplant donors with no chronic illness. Passages 4 to 7 ASM cells maintained in Ham's F12 medium supplemented with 10% FBS were used in all experiments.
## SRR1039512 Primary ASM cells were isolated from four white aborted lung transplant donors with no chronic illness. Passages 4 to 7 ASM cells maintained in Ham's F12 medium supplemented with 10% FBS were used in all experiments.
## SRR1039513 Primary ASM cells were isolated from four white aborted lung transplant donors with no chronic illness. Passages 4 to 7 ASM cells maintained in Ham's F12 medium supplemented with 10% FBS were used in all experiments.
## SRR1039514 Primary ASM cells were isolated from four white aborted lung transplant donors with no chronic illness. Passages 4 to 7 ASM cells maintained in Ham's F12 medium supplemented with 10% FBS were used in all experiments.
## SRR1039515 Primary ASM cells were isolated from four white aborted lung transplant donors with no chronic illness. Passages 4 to 7 ASM cells maintained in Ham's F12 medium supplemented with 10% FBS were used in all experiments.
## SRR1039516 Primary ASM cells were isolated from four white aborted lung transplant donors with no chronic illness. Passages 4 to 7 ASM cells maintained in Ham's F12 medium supplemented with 10% FBS were used in all experiments.
## SRR1039517 Primary ASM cells were isolated from four white aborted lung transplant donors with no chronic illness. Passages 4 to 7 ASM cells maintained in Ham's F12 medium supplemented with 10% FBS were used in all experiments.
## SRR1039518 Primary ASM cells were isolated from four white aborted lung transplant donors with no chronic illness. Passages 4 to 7 ASM cells maintained in Ham's F12 medium supplemented with 10% FBS were used in all experiments.
## SRR1039519 Primary ASM cells were isolated from four white aborted lung transplant donors with no chronic illness. Passages 4 to 7 ASM cells maintained in Ham's F12 medium supplemented with 10% FBS were used in all experiments.
## SRR1039520 Primary ASM cells were isolated from four white aborted lung transplant donors with no chronic illness. Passages 4 to 7 ASM cells maintained in Ham's F12 medium supplemented with 10% FBS were used in all experiments.
## SRR1039521 Primary ASM cells were isolated from four white aborted lung transplant donors with no chronic illness. Passages 4 to 7 ASM cells maintained in Ham's F12 medium supplemented with 10% FBS were used in all experiments.
## SRR1039522 Primary ASM cells were isolated from four white aborted lung transplant donors with no chronic illness. Passages 4 to 7 ASM cells maintained in Ham's F12 medium supplemented with 10% FBS were used in all experiments.
## SRR1039523 Primary ASM cells were isolated from four white aborted lung transplant donors with no chronic illness. Passages 4 to 7 ASM cells maintained in Ham's F12 medium supplemented with 10% FBS were used in all experiments.
## dex Run avgLength Experiment Sample BioSample
## SRR1039508 untrt SRR1039508 126 SRX384345 SRS508568 SAMN02422669
## SRR1039509 untrt SRR1039509 126 SRX384346 SRS508567 SAMN02422675
## SRR1039510 untrt SRR1039510 126 SRX384347 SRS508570 SAMN02422668
## SRR1039511 untrt SRR1039511 126 SRX384348 SRS508569 SAMN02422667
## SRR1039512 untrt SRR1039512 126 SRX384349 SRS508571 SAMN02422678
## SRR1039513 untrt SRR1039513 87 SRX384350 SRS508572 SAMN02422670
## SRR1039514 untrt SRR1039514 126 SRX384351 SRS508574 SAMN02422681
## SRR1039515 untrt SRR1039515 114 SRX384352 SRS508573 SAMN02422671
## SRR1039516 untrt SRR1039516 120 SRX384353 SRS508575 SAMN02422682
## SRR1039517 untrt SRR1039517 126 SRX384354 SRS508576 SAMN02422673
## SRR1039518 untrt SRR1039518 126 SRX384355 SRS508578 SAMN02422679
## SRR1039519 untrt SRR1039519 107 SRX384356 SRS508577 SAMN02422672
## SRR1039520 untrt SRR1039520 101 SRX384357 SRS508579 SAMN02422683
## SRR1039521 untrt SRR1039521 98 SRX384358 SRS508580 SAMN02422677
## SRR1039522 untrt SRR1039522 125 SRX384359 SRS508582 SAMN02422680
## SRR1039523 untrt 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.4.2 (2017-09-28)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.6-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.6-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] bindrcpp_0.2 airway_0.112.0
## [3] SummarizedExperiment_1.8.0 DelayedArray_0.4.0
## [5] matrixStats_0.52.2 GenomicRanges_1.30.0
## [7] GenomeInfoDb_1.14.0 IRanges_2.12.0
## [9] S4Vectors_0.16.0 GEOquery_2.46.0
## [11] Biobase_2.38.0 BiocGenerics_0.24.0
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.13 compiler_3.4.2
## [3] bindr_0.1 XVector_0.18.0
## [5] bitops_1.0-6 tools_3.4.2
## [7] zlibbioc_1.24.0 lattice_0.20-35
## [9] evaluate_0.10.1 tibble_1.3.4
## [11] pkgconfig_2.0.1 rlang_0.1.2
## [13] Matrix_1.2-11 GenomeInfoDbData_0.99.1
## [15] dplyr_0.7.4 httr_1.3.1
## [17] stringr_1.2.0 xml2_1.1.1
## [19] knitr_1.17 hms_0.3
## [21] tidyselect_0.2.2 grid_3.4.2
## [23] glue_1.2.0 R6_2.2.2
## [25] XML_3.98-1.9 tidyr_0.7.2
## [27] purrr_0.2.4 readr_1.1.1
## [29] magrittr_1.5 assertthat_0.2.0
## [31] stringi_1.1.5 RCurl_1.95-4.8