The data used in this workflow comes from an RNA-seq experiment where airway smooth muscle cells were treated with dexamethasone, a synthetic glucocorticoid steroid with anti-inflammatory effects (Himes et al. 2014). Glucocorticoids are used, for example, by people with asthma to reduce inflammation of the airways. In the experiment, four human airway smooth muscle cell lines were treated with 1 micromolar dexamethasone for 18 hours. For each of the four cell lines, we have a treated and an untreated sample. For more description of the experiment see the article, PubMed entry 24926665, and for raw data see the GEO entry GSE52778.
The information which connects the sample information from GEO with the SRA run id is downloaded from SRA using the Send to: File button.
For illustration purposes we will download subsampled fastq files from the course github site.
timeout <- options()$timeout
options(timeout = 300) #prevent timeout
url <- "https://github.com/statOmics/SGA/archive/refs/heads/airwaySeqData.zip"
destFile <- "airway.zip"
download.file(url, destFile)
unzip(destFile, exdir = "./", overwrite = TRUE)
if (file.exists(destFile)) {
#Delete file if it exists
file.remove(destFile)
}## [1] TRUEAssess the fastq files with fastQC.
We will make the input file for the hisat aligner that will be called through the QUASAR package.
fastqFiles <- list.files(path = "SGA-airwaySeqData/fastQ",full.names = TRUE, pattern = "fastq")
forward <- fastqFiles %>% grepl(pattern="1.fastq")
reverse <- fastqFiles %>% grepl(pattern="2.fastq")
samples <- fastqFiles[forward] |>
substr(start=25,stop=34)
fileInfo <- tibble(FileName1=fastqFiles[forward],FileName2=fastqFiles[reverse],SampleName=samples)
fileInfoWe write the fileInfo to disk
All reads in the subsampled fastq files map to chromosome 1. So download genome and annotation data for chromosome 1 only.
timeout <- options()$timeout
options(timeout = 500) #prevent timeout
urlGenome <- "https://ftp.ensembl.org/pub/release-113/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.1.fa.gz"
genomeDestFile <- "Homo_sapiens.GRCh38.dna.chromosome.1.fa.gz"
urlGff <- "https://ftp.ensembl.org/pub/release-113/gff3/homo_sapiens/Homo_sapiens.GRCh38.113.chromosome.1.gff3.gz"
gffDestFile <- "Homo_sapiens.GRCh38.113.chromosome.1.gff3.gz"
download.file(url = urlGenome ,
destfile = genomeDestFile)
download.file(url = urlGff,
destfile = gffDestFile)
gunzip(genomeDestFile, overwrite = TRUE, remove=TRUE)
gunzip(gffDestFile, overwrite = TRUE, remove=TRUE)
list.files()## [1] "_session-info.Rmd"
## [2] "_site.yml"
## [3] "about.html"
## [4] "about.Rmd"
## [5] "airway_files"
## [6] "airway_salmon_edgeR_files"
## [7] "airway_salmon_edgeR.html"
## [8] "airway_salmon_edgeR.Rmd"
## [9] "airway.html"
## [10] "airway.rmd"
## [11] "airwayCountTable.csv"
## [12] "airwayMapping.Rmd"
## [13] "bibliography.bib"
## [14] "Caenorhabditis_elegans.WBcel235.108.gff3"
## [15] "Caenorhabditis_elegans.WBcel235.dna.toplevel.fa"
## [16] "Caenorhabditis_elegans.WBcel235.dna.toplevel.fa.fai"
## [17] "Caenorhabditis_elegans.WBcel235.dna.toplevel.fa.md5"
## [18] "Caenorhabditis_elegans.WBcel235.dna.toplevel.fa.Rhisat2"
## [19] "cancer3x3_files"
## [20] "cancer3x3.Rmd"
## [21] "cancer6x6_files"
## [22] "cancer6x6.Rmd"
## [23] "cancer9x9_files"
## [24] "cancer9x9.Rmd"
## [25] "cptac_maxLFQ_files"
## [26] "cptac_maxLFQ.Rmd"
## [27] "cptac_median_files"
## [28] "cptac_median.Rmd"
## [29] "cptac_robust_files"
## [30] "cptac_robust_gui.Rmd"
## [31] "cptac_robust.Rmd"
## [32] "docs"
## [33] "elegansAlignmentCountTable.Rmd"
## [34] "elegansCountTable.csv"
## [35] "elegansDE.Rmd"
## [36] "fastQfiles.txt"
## [37] "figs"
## [38] "figures"
## [39] "figures.tar.gz"
## [40] "footer.html"
## [41] "heartMainInteraction.Rmd"
## [42] "Homo_sapiens.GRCh38.113.chromosome.1.gff3"
## [43] "Homo_sapiens.GRCh38.dna.chromosome.1.fa"
## [44] "Homo_sapiens.GRCh38.dna.chromosome.1.fa.fai"
## [45] "Homo_sapiens.GRCh38.dna.chromosome.1.fa.md5"
## [46] "Homo_sapiens.GRCh38.dna.chromosome.1.fa.Rhisat2"
## [47] "images_sequencing"
## [48] "index.Rmd"
## [49] "install.R"
## [50] "intro.Rmd"
## [51] "msqrob2.bib"
## [52] "multipleRegression_KPNA2_files"
## [53] "multipleRegression_KPNA2.Rmd"
## [54] "mzIDMsgfSwissprotExample.Rmd"
## [55] "mzIDswissprotTutorial.Rmd"
## [56] "mzIDuniprotTutorial.Rmd"
## [57] "parathyroid_files"
## [58] "parathyroid.Rmd"
## [59] "pda_blocking_wrapup_files"
## [60] "pda_blocking_wrapup.Rmd"
## [61] "pda_quantification_inference_files"
## [62] "pda_quantification_inference.Rmd"
## [63] "pda_quantification_preprocessing_files"
## [64] "pda_quantification_preprocessing.Rmd"
## [65] "pda_robustSummarisation_peptideModels_files"
## [66] "pda_robustSummarisation_peptideModels.Rmd"
## [67] "pda_technicalDetails.Rmd"
## [68] "pda_tutorialDesign.Rmd"
## [69] "pda_tutorialPreprocessing.Rmd"
## [70] "pyrococcusMSGF+.mzid"
## [71] "pyroSwissprot.mzid"
## [72] "pyroUniprot.mzid"
## [73] "QuasR_log_debd2f67fea0.txt"
## [74] "QuasR_log_e34943403f50.txt"
## [75] "README.md"
## [76] "recapGeneralLinearModel_files"
## [77] "recapGeneralLinearModel.Rmd"
## [78] "rmdBub"
## [79] "scRNA-seq-Tutorial.Rmd"
## [80] "sequencing_countData_files"
## [81] "sequencing_countData.log"
## [82] "sequencing_countData.Rmd"
## [83] "sequencing_intro.log"
## [84] "sequencing_intro.Rmd"
## [85] "sequencing_singleCell_files"
## [86] "sequencing_singleCell.Rmd"
## [87] "sequencing_technicalDE_files"
## [88] "sequencing_technicalDE.Rmd"
## [89] "sequencingTutorial1_mapping.Rmd"
## [90] "sequencingTutorial2_DE.Rmd"
## [91] "sequencingTutorial3_DE.Rmd"
## [92] "SGA-airwaySeqData"
## [93] "SGA-airwaySeqData.zip"
## [94] "SGA-elegansFastq"
## [95] "singleCell_MacoskoWorkflow_blank.Rmd"
## [96] "site_libs"
## [97] "style.css"
## [98] "TargetDecoy.Rmd"
## [99] "TargetDecoyTutorial.Rmd"We are aligning RNA-seq reads and have to use a gap aware
alignment modus. We therefore use the argument
splicedAlignment = TRUE
The project is a single-end sequencing project! So we use the
argument paired = "no".
suppressPackageStartupMessages({
library(QuasR)
library(Rhisat2)
})
sampleFile <- "fastQfiles.txt"
genomeFile <- "Homo_sapiens.GRCh38.dna.chromosome.1.fa"
projAirway <- qAlign(
sampleFile = sampleFile,
genome = genomeFile,
splicedAlignment = TRUE,
aligner = "Rhisat2",
paired="fr")## all necessary alignment files foundThe gff3 file contains information on the position of features, e.g. exons and genes, along the chromosome.
suppressPackageStartupMessages({
library(GenomicFeatures)
library(Rsamtools)
})
annotFile <- "Homo_sapiens.GRCh38.113.chromosome.1.gff3"
#chrLen <- scanFaIndex(genomeFile)
#chrominfo <- data.frame(chrom = as.character(seqnames(chrLen)),
# length = width(chrLen),
# is_circular = rep(FALSE, length(chrLen)))
txdb <- makeTxDbFromGFF(file = annotFile,
dataSource = "Ensembl",
organism = "Homo sapiens")## Warning in call_fun_in_txdbmaker("makeTxDbFromGFF", ...): makeTxDbFromGFF() has moved to the txdbmaker package. Please call
## txdbmaker::makeTxDbFromGFF() to get rid of this warning.## Import genomic features from the file as a GRanges object ... OK
## Prepare the 'metadata' data frame ... OK
## Make the TxDb object ... OKWith the qCount function we can count the overlap between the aligned reads and the genomic features of interest.
## extracting gene regions from TxDb...done
## counting alignments...done
## collapsing counts by query name...done## width SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516
## ENSG00000000457 6883 0 0 0 0 0
## ENSG00000000460 5970 0 0 0 0 0
## ENSG00000000938 3382 0 0 0 0 0
## ENSG00000000971 15284 0 0 0 0 0
## ENSG00000001460 8511 0 0 0 0 0
## ENSG00000001461 11672 0 0 0 0 0
## SRR1039517 SRR1039520 SRR1039521
## ENSG00000000457 0 0 0
## ENSG00000000460 0 0 0
## ENSG00000000938 0 0 0
## ENSG00000000971 0 0 0
## ENSG00000001460 0 0 0
## ENSG00000001461 0 0 0Save count table for future use.
With respect to reproducibility, it is highly recommended to include a session info in your script so that readers of your output can see your particular setup of R.
## R version 4.4.0 RC (2024-04-16 r86468)
## Platform: aarch64-apple-darwin20
## Running under: macOS Big Sur 11.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Europe/Brussels
## tzcode source: internal
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] Rsamtools_2.20.0 Biostrings_2.72.1 XVector_0.44.0
## [4] GenomicFeatures_1.56.0 AnnotationDbi_1.66.0 Biobase_2.64.0
## [7] Rhisat2_1.20.0 QuasR_1.44.0 Rbowtie_1.44.0
## [10] GenomicRanges_1.56.1 GenomeInfoDb_1.40.1 IRanges_2.38.1
## [13] S4Vectors_0.42.1 BiocGenerics_0.50.0 R.utils_2.12.3
## [16] R.oo_1.26.0 R.methodsS3_1.8.2 lubridate_1.9.3
## [19] forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
## [22] purrr_1.0.2 readr_2.1.5 tidyr_1.3.1
## [25] tibble_3.2.1 ggplot2_3.5.1 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] DBI_1.2.3 bitops_1.0-8
## [3] deldir_2.0-4 httr2_1.0.5
## [5] biomaRt_2.60.1 rlang_1.1.4
## [7] magrittr_2.0.3 matrixStats_1.4.1
## [9] compiler_4.4.0 RSQLite_2.3.7
## [11] txdbmaker_1.0.1 png_0.1-8
## [13] vctrs_0.6.5 pwalign_1.0.0
## [15] pkgconfig_2.0.3 crayon_1.5.3
## [17] fastmap_1.2.0 dbplyr_2.5.0
## [19] utf8_1.2.4 rmarkdown_2.28
## [21] tzdb_0.4.0 UCSC.utils_1.0.0
## [23] bit_4.5.0 xfun_0.47
## [25] zlibbioc_1.50.0 cachem_1.1.0
## [27] jsonlite_1.8.9 progress_1.2.3
## [29] blob_1.2.4 DelayedArray_0.30.1
## [31] BiocParallel_1.38.0 jpeg_0.1-10
## [33] prettyunits_1.2.0 R6_2.5.1
## [35] VariantAnnotation_1.50.0 bslib_0.8.0
## [37] stringi_1.8.4 RColorBrewer_1.1-3
## [39] rtracklayer_1.64.0 jquerylib_0.1.4
## [41] Rcpp_1.0.13-1 SummarizedExperiment_1.34.0
## [43] knitr_1.48 SGSeq_1.38.0
## [45] igraph_2.0.3 Matrix_1.7-0
## [47] timechange_0.3.0 tidyselect_1.2.1
## [49] rstudioapi_0.16.0 abind_1.4-8
## [51] yaml_2.3.10 codetools_0.2-20
## [53] RUnit_0.4.33 hwriter_1.3.2.1
## [55] curl_5.2.3 lattice_0.22-6
## [57] withr_3.0.1 ShortRead_1.62.0
## [59] KEGGREST_1.44.1 evaluate_1.0.0
## [61] BiocFileCache_2.12.0 xml2_1.3.6
## [63] filelock_1.0.3 pillar_1.9.0
## [65] MatrixGenerics_1.16.0 generics_0.1.3
## [67] vroom_1.6.5 RCurl_1.98-1.16
## [69] hms_1.1.3 munsell_0.5.1
## [71] scales_1.3.0 GenomicFiles_1.40.0
## [73] glue_1.8.0 tools_4.4.0
## [75] interp_1.1-6 BiocIO_1.14.0
## [77] BSgenome_1.72.0 GenomicAlignments_1.40.0
## [79] XML_3.99-0.17 grid_4.4.0
## [81] latticeExtra_0.6-30 colorspace_2.1-1
## [83] GenomeInfoDbData_1.2.12 restfulr_0.0.15
## [85] cli_3.6.3 rappdirs_0.3.3
## [87] fansi_1.0.6 S4Arrays_1.4.1
## [89] gtable_0.3.5 sass_0.4.9
## [91] digest_0.6.37 SparseArray_1.4.8
## [93] rjson_0.2.23 memoise_2.0.1
## [95] htmltools_0.5.8.1 lifecycle_1.0.4
## [97] httr_1.4.7 bit64_4.5.2