From 29e800fa80bbd1026f72ec1fe4ac0dc093238bc4 Mon Sep 17 00:00:00 2001
From: Marcel Kempenaar <m.kempenaar@pl.hanze.nl>
Date: Mon, 9 Sep 2024 08:16:00 +0200
Subject: [PATCH] Added project proposal as deliverable

---
 docs/404.html                  | 2 +-
 docs/EDA.html                  | 2 +-
 docs/a2-annotation.html        | 2 +-
 docs/a2-data_loading.html      | 2 +-
 docs/chapter-4.html            | 2 +-
 docs/chapter-5.html            | 2 +-
 docs/datasets.html             | 2 +-
 docs/from-reads-to-counts.html | 2 +-
 docs/index.html                | 8 ++++++--
 docs/search.json               | 2 +-
 index.Rmd                      | 4 +++-
 11 files changed, 18 insertions(+), 12 deletions(-)

diff --git a/docs/404.html b/docs/404.html
index ea6e2d0..f865f72 100644
--- a/docs/404.html
+++ b/docs/404.html
@@ -101,7 +101,7 @@ <h1>Page not found<a class="anchor" aria-label="anchor" href="#page-not-found"><
 <footer class="bg-primary text-light mt-5"><div class="container"><div class="row">
 
   <div class="col-12 col-md-6 mt-3">
-    <p>"<strong>2.1.2 Genomics &amp; Transcriptomics</strong>" was written by Marcel Kempenaar. It was last built on 2024-09-08.</p>
+    <p>"<strong>2.1.2 Genomics &amp; Transcriptomics</strong>" was written by Marcel Kempenaar. It was last built on 2024-09-09.</p>
   </div>
 
   <div class="col-12 col-md-6 mt-3">
diff --git a/docs/EDA.html b/docs/EDA.html
index df5de47..aa1db98 100644
--- a/docs/EDA.html
+++ b/docs/EDA.html
@@ -716,7 +716,7 @@ <h2>
 <footer class="bg-primary text-light mt-5"><div class="container"><div class="row">
 
   <div class="col-12 col-md-6 mt-3">
-    <p>"<strong>2.1.2 Genomics &amp; Transcriptomics</strong>" was written by Marcel Kempenaar. It was last built on 2024-09-08.</p>
+    <p>"<strong>2.1.2 Genomics &amp; Transcriptomics</strong>" was written by Marcel Kempenaar. It was last built on 2024-09-09.</p>
   </div>
 
   <div class="col-12 col-md-6 mt-3">
diff --git a/docs/a2-annotation.html b/docs/a2-annotation.html
index ac0ef00..44ee901 100644
--- a/docs/a2-annotation.html
+++ b/docs/a2-annotation.html
@@ -899,7 +899,7 @@ <h3>
 <footer class="bg-primary text-light mt-5"><div class="container"><div class="row">
 
   <div class="col-12 col-md-6 mt-3">
-    <p>"<strong>2.1.2 Genomics &amp; Transcriptomics</strong>" was written by Marcel Kempenaar. It was last built on 2024-09-08.</p>
+    <p>"<strong>2.1.2 Genomics &amp; Transcriptomics</strong>" was written by Marcel Kempenaar. It was last built on 2024-09-09.</p>
   </div>
 
   <div class="col-12 col-md-6 mt-3">
diff --git a/docs/a2-data_loading.html b/docs/a2-data_loading.html
index 4edb939..88c3598 100644
--- a/docs/a2-data_loading.html
+++ b/docs/a2-data_loading.html
@@ -474,7 +474,7 @@ <h3>
 <footer class="bg-primary text-light mt-5"><div class="container"><div class="row">
 
   <div class="col-12 col-md-6 mt-3">
-    <p>"<strong>2.1.2 Genomics &amp; Transcriptomics</strong>" was written by Marcel Kempenaar. It was last built on 2024-09-08.</p>
+    <p>"<strong>2.1.2 Genomics &amp; Transcriptomics</strong>" was written by Marcel Kempenaar. It was last built on 2024-09-09.</p>
   </div>
 
   <div class="col-12 col-md-6 mt-3">
diff --git a/docs/chapter-4.html b/docs/chapter-4.html
index a1da9b9..a630d38 100644
--- a/docs/chapter-4.html
+++ b/docs/chapter-4.html
@@ -489,7 +489,7 @@ <h3>
 <footer class="bg-primary text-light mt-5"><div class="container"><div class="row">
 
   <div class="col-12 col-md-6 mt-3">
-    <p>"<strong>2.1.2 Genomics &amp; Transcriptomics</strong>" was written by Marcel Kempenaar. It was last built on 2024-09-08.</p>
+    <p>"<strong>2.1.2 Genomics &amp; Transcriptomics</strong>" was written by Marcel Kempenaar. It was last built on 2024-09-09.</p>
   </div>
 
   <div class="col-12 col-md-6 mt-3">
diff --git a/docs/chapter-5.html b/docs/chapter-5.html
index bcf77b2..3351315 100644
--- a/docs/chapter-5.html
+++ b/docs/chapter-5.html
@@ -318,7 +318,7 @@ <h3>
 <footer class="bg-primary text-light mt-5"><div class="container"><div class="row">
 
   <div class="col-12 col-md-6 mt-3">
-    <p>"<strong>2.1.2 Genomics &amp; Transcriptomics</strong>" was written by Marcel Kempenaar. It was last built on 2024-09-08.</p>
+    <p>"<strong>2.1.2 Genomics &amp; Transcriptomics</strong>" was written by Marcel Kempenaar. It was last built on 2024-09-09.</p>
   </div>
 
   <div class="col-12 col-md-6 mt-3">
diff --git a/docs/datasets.html b/docs/datasets.html
index 0791af7..9dffd55 100644
--- a/docs/datasets.html
+++ b/docs/datasets.html
@@ -192,7 +192,7 @@ <h2>
 <footer class="bg-primary text-light mt-5"><div class="container"><div class="row">
 
   <div class="col-12 col-md-6 mt-3">
-    <p>"<strong>2.1.2 Genomics &amp; Transcriptomics</strong>" was written by Marcel Kempenaar. It was last built on 2024-09-08.</p>
+    <p>"<strong>2.1.2 Genomics &amp; Transcriptomics</strong>" was written by Marcel Kempenaar. It was last built on 2024-09-09.</p>
   </div>
 
   <div class="col-12 col-md-6 mt-3">
diff --git a/docs/from-reads-to-counts.html b/docs/from-reads-to-counts.html
index 871b086..415022a 100644
--- a/docs/from-reads-to-counts.html
+++ b/docs/from-reads-to-counts.html
@@ -483,7 +483,7 @@ <h3>
 <footer class="bg-primary text-light mt-5"><div class="container"><div class="row">
 
   <div class="col-12 col-md-6 mt-3">
-    <p>"<strong>2.1.2 Genomics &amp; Transcriptomics</strong>" was written by Marcel Kempenaar. It was last built on 2024-09-08.</p>
+    <p>"<strong>2.1.2 Genomics &amp; Transcriptomics</strong>" was written by Marcel Kempenaar. It was last built on 2024-09-09.</p>
   </div>
 
   <div class="col-12 col-md-6 mt-3">
diff --git a/docs/index.html b/docs/index.html
index 58f64b0..0f17825 100644
--- a/docs/index.html
+++ b/docs/index.html
@@ -98,6 +98,8 @@ <h3>
 </h3>
 <ul>
 <li>
+<strong>Project Proposal</strong>: A short written proposal of the project clearly describing the data set, the research question and the methods + tools to be used.</li>
+<li>
 <strong>Genomics pipeline</strong>: The genomics analysis pipeline is well documented in one or more <em>lab journals</em> and contains the correct tools based on the research question and data set.</li>
 <li>
 <strong>Presentation</strong>: The research results will be presented in an attractive method and proper level given the target audience.</li>
@@ -141,7 +143,7 @@ <h3>
 </li>
 </ul>
 <p>Data is available from the NCBI <a href="https://www.ncbi.nlm.nih.gov/sra">Sequence Read Archive</a> and Gene Expression Omnibus.</p>
-<p>The proper way to search for a data set suitable for both Genomics and Transcriptomics projects is by creating a query that searches for both: <a href="https://www.ncbi.nlm.nih.gov/gds/?term=(%22expression+profiling+by+high+throughput+sequencing%22%5BDataSet+Type%5D)+AND+%22genome+variation+profiling+by+high+throughput+sequencing%22%5BDataSet+Type%5D">(“expression profiling by high throughput sequencing”[DataSet Type]) AND “genome variation profiling by high throughput sequencing”[DataSet Type]</a>.</p>
+<p>The proper way to search for a data set suitable for both Genomics and Transcriptomics projects is by creating a query that searches for both: <a href="https://www.ncbi.nlm.nih.gov/gds/?term=(%22expression+profiling+by+high+throughput+sequencing%22%5BDataSet+Type%5D)+AND+%22genome+variation+profiling+by+high+throughput+sequencing%22%5BDataSet+Type%5D">(“expression profiling by high throughput sequencing”[DataSet Type]) AND “genome variation profiling by high throughput sequencing”[DataSet Type]</a>. If this selection actually contains anything of interest is up to you; if you cannot find something interesting please use a different search strategy. One could be to search for RNA-Seq experiments in the GEO where one of the groups is preferably a mutant.</p>
 <p>Note that for the genomics part, not all available samples have to be processed. Other than for transcriptomics, it is not common that there is a case/control experimental setup but rather a single condition of group with possibly a few replicates. If the data set does include many samples, it is recommended to select a subset of samples to process (do discuss the reasoning behind this selection).</p>
 </div>
 </div>
@@ -157,6 +159,8 @@ <h3>
 </h3>
 <ul>
 <li>
+<strong>Project Proposal</strong>: A short written proposal of the project clearly describing the data set, the research question and the methods + libraries to be used.</li>
+<li>
 <strong>Lab Journal</strong>: The lab journal is readable, complete and aimed at the reproducability of the research. Relevant visualisations are include and statistical methods performed to support the biological interpretation of the results.</li>
 </ul>
 </div>
@@ -574,7 +578,7 @@ <h2>
 <footer class="bg-primary text-light mt-5"><div class="container"><div class="row">
 
   <div class="col-12 col-md-6 mt-3">
-    <p>"<strong>2.1.2 Genomics &amp; Transcriptomics</strong>" was written by Marcel Kempenaar. It was last built on 2024-09-08.</p>
+    <p>"<strong>2.1.2 Genomics &amp; Transcriptomics</strong>" was written by Marcel Kempenaar. It was last built on 2024-09-09.</p>
   </div>
 
   <div class="col-12 col-md-6 mt-3">
diff --git a/docs/search.json b/docs/search.json
index 83d1d8c..833c704 100644
--- a/docs/search.json
+++ b/docs/search.json
@@ -1 +1 @@
-[{"path":"index.html","id":"chapter-1","chapter":" Introduction","heading":" Introduction","text":"Course 2.1.2 involves two separate projects performed single, published, data set. First, start performing Genomics analysis similar case study module 2.1.1. , take result mapping step, process read counts perform subsequent Transcriptomics analysis.","code":""},{"path":"index.html","id":"genomics","chapter":" Introduction","heading":"1.1 2.1.2 Genomics","text":"Genomics analysis going look variants chosen data reference genome. One important aspect analysis often assessing quality, input data set well output certain steps. set variants determined can annotated finally linked biology underlying original experiment, either reproducing findings aiming novel ideas.","code":""},{"path":"index.html","id":"project-deliverables","chapter":" Introduction","heading":"1.1.1 Project Deliverables","text":"Genomics pipeline: genomics analysis pipeline well documented one lab journals contains correct tools based research question data set.Presentation: research results presented attractive method proper level given target audience.","code":""},{"path":"index.html","id":"data-processing","chapter":" Introduction","heading":"1.1.2 Data Processing","text":"Similar 2.1.1 NGS & Mapping course go similar steps process data (.e. quality control, trimming, mapping, variant calling, etc.). time however, using workflow manager Galaxy, execute tools command line. Two reasons number size input files far larger compared previously used data sets Galaxy server capable handling many large data sets. Also, proper terminal experience valuable.Note writing pipeline: steps need documented RMarkdown documents. step must introduced short description tool used, version followed command used execute tool. terminal command, use code cells bash language. code cells must option eval=FALSE set. prevent code executed knitting document.Furthermore, also required run command RStudio server provided https://bioinf.nl/rstudio. RStudio server hosted assemblix server capable accessing data running tools large data sets. workstations provide classrooms capable tasks. Also, server already number tools pre-installed (FastQC, trimmomatic, BWA/STAR mappers, etc.). Lastly, running RMarkdown document RStudio server forces document commands compared simply using ssh.Note storing (intermediate) data: data relevant project must stored /students/2024-2025/Thema05/ directory. , create project directory team create sub-directories go. home folder hard limit 25GB data therefore sufficient storing project files. Use commands chmod chown allow project partners access data .Note processing large amounts data: building correct terminal commands execute instance trimming step large number files can quite challenging. therefore required test-data used steps. Reasons :Spotting errors quickly instead wait several hours see error messageInspecting resource usage (CPU, memory, disk space) prevent overloading system\nsteps can memory intensive (1 TB large data sets!) can crash system\nsteps can memory intensive (1 TB large data sets!) can crash systemCreates small output set can used test set next stepNote use GATK: Genome Analysis Toolkit (GATK) often used set tools variant calling. Many publications cite use GATK unfortunately include needed details reproduce results. Due complexity GATK recommend use, unless publication clear exact workflow settings used. Prepare spend lot time official GATK website browsing manuals can use effectively.","code":""},{"path":"index.html","id":"finding-a-data-set","chapter":" Introduction","heading":"1.1.3 Finding a Data Set","text":"Finding suitable data set can challenging requirements data set needs meet:needs DNA- RNA-Seq data set available.data set needs published (publication open access)RNA-Seq part, minimum amount samples required\nleast two sample groups minimum 3 samples per group (replicates)\nSee appendix chapter A1 example valid experimental setups, chapter also contains tips (mostly relevant Transcriptomics).\nleast two sample groups minimum 3 samples per group (replicates)See appendix chapter A1 example valid experimental setups, chapter also contains tips (mostly relevant Transcriptomics).Data available NCBI Sequence Read Archive Gene Expression Omnibus.proper way search data set suitable Genomics Transcriptomics projects creating query searches : (“expression profiling high throughput sequencing”[DataSet Type]) “genome variation profiling high throughput sequencing”[DataSet Type].Note genomics part, available samples processed. transcriptomics, common case/control experimental setup rather single condition group possibly replicates. data set include many samples, recommended select subset samples process (discuss reasoning behind selection).","code":""},{"path":"index.html","id":"transcriptomics","chapter":" Introduction","heading":"1.2 2.1.2 Transcriptomics","text":"use RNA-Sequencing techniques measuring gene expression relatively new replaces microarrays, though cases microarrays still used. Gene expression data gives valuable insights workings cells certain conditions. Especially comparing instance healthy diseased samples can become clear genes causal influence specific condition. Finding genes interest (genes showing differing expression accross conditions, called Differentially Expressed Genes (DEGs)) goal project.golden standard analyzing RNA-sequencing data sets many tools (manufacturers sequencing equipment also deliver software packages) use proven R libraries processing, visualizing analyzing publicly available data sets. manual describes steps often performed transcriptomics experiment. Use general guideline process data accompanying chosen published research. Note manual originally written processing count data therefore describe processing NGS read files count data part Genomics analysis performed course. RNA-Seq counts describe number reads mapped gene corresponds relative number transcripts (mRNA sequences) gene present cell time sampling.","code":""},{"path":"index.html","id":"project-deliverables-1","chapter":" Introduction","heading":"1.2.1 Project Deliverables","text":"Lab Journal: lab journal readable, complete aimed reproducability research. Relevant visualisations include statistical methods performed support biological interpretation results.","code":""},{"path":"index.html","id":"lab-journal","chapter":" Introduction","heading":"1.3 Lab Journal","text":"may know previous projects, essential keep proper journal detailing every step done analysis. journal kept R markdown file, showing steps taken analysis data set. markdown knitted single PDF-file project completed thus containing text detailing steps decisions ’ve made, R-code (always visible!) resulting output/ images.\nNotes: general advice; wait \nknitting whole document project done knitting \nprone errors trying fix large document \neasy. Give code chunk proper attributes, including name \nminimum. helps spot errors knitting process\nmentiones chunk processed. Note chunk-names\nmust unique. try make proper use chapters \nsections include (optional) table contents.\n","code":""},{"path":"index.html","id":"schedule","chapter":" Introduction","heading":"1.4 Schedule","text":"Genomics Transcriptomics parts schedule listed ","code":""},{"path":"index.html","id":"genomics-1","chapter":" Introduction","heading":"1.4.1 Genomics","text":"Empty rows means previous topic task likely take one session.","code":""},{"path":"index.html","id":"transcriptomics-1","chapter":" Introduction","heading":"1.4.2 Transcriptomics","text":"","code":""},{"path":"index.html","id":"literature","chapter":" Introduction","heading":"1.5 Literature","text":"following (online) documents can used throughout course:Gene Set Enrichment Analysis – Theory. Available : https://www.pathwaycommons.org/guide/primers/data_analysis/gsea/Li D. Statistical Methods RNA Sequencing Data Analysis. : Husi H, editor. Computational Biology [Internet]. Brisbane (AU): Codon Publications; 2019 Nov 21. Chapter 6. Available : https://www.ncbi.nlm.nih.gov/books/NBK550334/Yu, L., Fernandez, S. & Brock, G. Power analysis RNA-Seq differential expression studies. BMC Bioinformatics 18, 234 (2017). https://doi.org/10.1186/s12859-017-1648-2Sun, S., Hood, M., Scott, L., Peng, Q., Mukherjee, S., Tung, J., & Zhou, X. (2017). Differential expression analysis RNAseq using Poisson mixed models. Nucleic acids research, 45(11), e106. https://doi.org/10.1093/nar/gkx204","code":""},{"path":"from-reads-to-counts.html","id":"from-reads-to-counts","chapter":" From Reads to Counts","heading":" From Reads to Counts","text":"Step one RNA-Seq data analysis generating gene transcript counts sample. overlaps partly steps taken genomics workflow:Download SRA files extract FastQ filesPerform quality control raw reads (FastQC)Perform trimming needed (Trimmomatic)Perform quality control post trimming (FastQC)Perform mapping reference genome (often STAR)Generate gene/ transcript counts (either STAR, HTseq featurecounts)Downstream analysis (chapter 4 document)chapter shows steps example data set (GSE149995)","code":""},{"path":"from-reads-to-counts.html","id":"downloading-and-extracting-fastq-files","chapter":" From Reads to Counts","heading":"2.1 Downloading and extracting FastQ files","text":"Given SRA identifier project, select samples download accession file (SRX identifiers). Using identifiers input prefetch system command (requires NCBI SRA-Toolkit installed) downloads SRA files convert using fasterq-dump command fastq files.fasterq-dump default performs split 3 operation, see SRA toolkit manual details.","code":"prefetch $(<GSE149995_SRA_ACC_RNA-Seq.csv) --output-directory \\\n  /local-fs/staff/marcel/GSE149995/SRAfind /local-fs/staff/marcel/GSE149995/SRA/ -name \"*.sra\" | \\\n  parallel 'fasterq-dump -O /local-fs/staff/marcel/GSE149995/fastq/ {}'"},{"path":"from-reads-to-counts.html","id":"optional-annotation","chapter":" From Reads to Counts","heading":"2.2 (Optional) Annotation","text":"NCBI GEO page lists 9 samples 3 conditions (WT, potent arf7). However, experiment includes 18 sequence files. downloading ‘series matrix file’ GEO (using GEOquery library) combining run info downloadable SRA, can list mapping run sample.Note: optional step now useful doubts group/sample setup.","code":"\nlibrary(GEOquery)\ngse <- getGEO(GEO = \"GSE149995\")\n# Combine the condition with a sample number and replicate number\nCondition = paste(rep(gse[[1]]@phenoData@data$`genotype/variation:ch1`, each = 2),\n                  rep(1:3, each = 2),\n                  paste0('r', rep(1:2, 9)), sep = '_')\nrun_info <- read.csv(file = \"data/GSE149995_Sra_RunInfo.csv\")\nsetup <- cbind(run_info %>% dplyr::select(Run, Experiment), Condition)\npander(setup)"},{"path":"from-reads-to-counts.html","id":"quality-control","chapter":" From Reads to Counts","heading":"2.3 Quality Control","text":"perform FastQC -> Trimmomatic -> FastQC steps","code":""},{"path":"from-reads-to-counts.html","id":"fastqc","chapter":" From Reads to Counts","heading":"2.3.1 FastQC","text":"Running FastQCQC-report using multiqc","code":"\nfqdir <- \"/local-fs/staff/marcel/GSE149995/fastq/\"\nfastqc(fq.dir = \"/local-fs/staff/marcel/GSE149995/fastq/\",\n       qc.dir = \"/local-fs/staff/marcel/GSE149995/fqc\",\n       fastqc.path = \"/usr/bin/fastqc\",\n       threads = 32)multiqc \\\n  --filename ~/Development/2.1.2-Transcriptomics/GSE149995_FQC_multiqc_report.html \\\n  /local-fs/staff/marcel/GSE149995/fqc/ "},{"path":"from-reads-to-counts.html","id":"trimming","chapter":" From Reads to Counts","heading":"2.4 Trimming","text":"Proper trimming requires experimentation research regarding used sequencing technology. following command uses Trimmomatic (release 0.39, available Github) trim fastq files. Note example single-end sequencing , use TrimmomaticPE command paired-end data processing TruSeq FASTA file ending PE instead.Re-running FastQC","code":"cat data/GSE149995_Sra_RunInfo.csv | \\\n  parallel 'TrimmomaticSE -threads 4 ' \\\n                  '/local-fs/staff/marcel/GSE149995/fastq/{}.fastq.gz ' \\\n                  '/local-fs/staff/marcel/GSE149995/fastq-trimmed/{}.trimmed.fastq.gz ' \\\n                  'ILLUMINACLIP:/homes/marcelk/Development/2.1.2-Transcriptomics/TruSeq3-SE.fa:2:30:10 ' \\\n                  'MINLEN:40 ' \\\n                  'SLIDINGWINDOW:4:20'\nfastqc(fq.dir = \"/local-fs/staff/marcel/GSE149995/fastq-trimmed/\",\n       qc.dir = \"/local-fs/staff/marcel/GSE149995/fqc-trimmed\",\n       fastqc.path = \"/usr/bin/fastqc\",\n       threads = 32)multiqc --filename \\\n  ~/Development/2.1.2-Transcriptomics/GSE149995_FQC_Trimmed.html \\\n  /local-fs/staff/marcel/GSE149995/fqc-trimmed/"},{"path":"from-reads-to-counts.html","id":"mapping-with-star","chapter":" From Reads to Counts","heading":"2.5 Mapping with STAR","text":"","code":""},{"path":"from-reads-to-counts.html","id":"building-star-index","chapter":" From Reads to Counts","heading":"2.5.1 Building STAR index","text":"bwa, STAR requires extra annotation file (GFF GTF format) describing features genome. data set requires FASTA Arabidopsis thaliana genome accompanying annotation file. genomeDir specifies output location index.","code":"cd /local-fs/staff/marcel/GSE149995/reference\nmkdir star; cd star\nSTAR \\\n  --runThreadN 12 \\\n  --runMode genomeGenerate \\\n  --genomeDir /local-fs/staff/marcel/GSE149995/reference/star \\\n  --genomeFastaFiles /local-fs/staff/marcel/GSE149995/reference/TAIR10_chr_all.fas \\\n  --sjdbGTFfile /local-fs/staff/marcel/GSE149995/reference/TAIR10_GTF_genes.gtf \\\n  --genomeSAindexNbases 12 \\\n  --sjdbOverhang 49"},{"path":"from-reads-to-counts.html","id":"read-mapping-and-transcript-counting","chapter":" From Reads to Counts","heading":"2.5.2 Read mapping and transcript counting","text":"following command performs mapping reference genome, saves sorted BAM file directly performs read count transcripts (NOTE: requires --sjdbGTFfile parameter index building). input file, ReadsPerGenes..tab file generated containing 4 columns:column 1: gene IDcolumn 2: counts unstranded RNA-seqcolumn 3: counts 1st read strand aligned RNA (htseq-count option -s yes)column 4: counts 2nd read strand aligned RNA (htseq-count option -s reverse)See interesting discussion BioStars regarding column use.planning use different tool calculating gene counts, can remove quantMode argument generate sorted BAM files. Please read manual tool choice see indexed BAM file BAM file index.Generating MultiQC report STAR mapping phase","code":"find /local-fs/staff/marcel/GSE149995/fastq-trimmed/ -name \"*.fastq.gz\" -exec basename {} \\; | \\\n  sed 's/\\.fastq.gz$//' | parallel 'STAR --runThreadN 6 ' \\\n    '--genomeDir /local-fs/staff/marcel/GSE149995/reference/star ' \\\n    '--readFilesIn /local-fs/staff/marcel/GSE149995/fastq-trimmed/{}.fastq ' \\\n    '--outSAMtype BAM SortedByCoordinate ' \\\n    '--quantMode GeneCounts ' \\\n    '--outFileNamePrefix /local-fs/staff/marcel/GSE149995/mapping/star/{}_star_'~/miniconda3/bin/multiqc --filename \\\n  ~/Development/2.1.2-Transcriptomics/STAR-mapping.html \\\n  /local-fs/staff/marcel/GSE149995/mapping/star/"},{"path":"from-reads-to-counts.html","id":"combining-star-count-files","chapter":" From Reads to Counts","heading":"2.5.3 Combining STAR Count Files","text":"STAR generates [sample]_star_ReadsPerGene..tab file sample, containing gene IDs count numbers. can merge files single data frame processing (see also Appendix A1 manual). example takes 2nd column count values. Make sure know column use based description . reason ’re merging data might always guaranteed get exact transcripts file.Note: use featurecounts tool, already get single file containing count values samples.Now read count files merge single data frame.now single data frame count values samples transcripts. names however nondescriptive since just SRR identifiers. Since already annotation, ’re going rename columns something meaningful:","code":"\nfile.names <- list.files('/local-fs/staff/marcel/GSE149995/mapping/star/',\n                         pattern = '*_star_ReadsPerGene.out.tab')\n\n## Function for reading in files\nread_sample <- function(file.name) {\n  ## Extract the sample name for naming the column (retaining the 'SRR....' part)\n  sample.name <- strsplit(file.name, \".\", fixed = TRUE)[[1]][1]\n  sample <- read.table(file.name, header = FALSE, sep=\"\\t\", \n                       row.names = NULL, skip = 4)\n  ## Rename the count column\n  names(sample)[2] <- sample.name\n  ## Return a subset containing the transcript ID and sample name columns\n  return(sample[c(1, 2)])\n}\nsetwd('/local-fs/staff/marcel/GSE149995/mapping/star/')\n\n## Read the FIRST sample\ncounts <- read_sample(file.names[1])\n\n## Read the remaining files and merge the contents\nfor (file.name in file.names[2:length(file.names)]) {\n  sample <- read_sample(file.name)\n  counts <- merge(counts, sample, by = 1)\n}\n\n# Set the row names to the transcript IDs\nrownames(counts) <- counts$V1\ncounts <- counts[-1]\n# For each SRR id, find the corresponding column in counts and rename\n# using the name in the 'Condition' column\nfor (i in 1:nrow(setup)) {\n  idx <- grep(setup$Run[i], names(counts))\n  names(counts)[idx] <- setup$Condition[i]\n}\n\n# Show first 6 columns and rows\npander(counts[1:6, 1:6])"},{"path":"EDA.html","id":"EDA","chapter":" Exploratory Data Analysis","heading":" Exploratory Data Analysis","text":"start performing exploratory data analysis steps goal getting grips chosen data set properly identify strategy actual analysis steps. exploration also keep eye quality data. Even though downloadable data ‘processed’, might samples present deviate much group samples (called outlier). Creating basic visualizations data give necessary insight continue.","code":""},{"path":"EDA.html","id":"loading-data","chapter":" Exploratory Data Analysis","heading":"3.1 Loading data into R","text":"one column based text files can read R simply using read.table() function.Follow following steps read data start exploratory data analysis. resulting document treated lab journal log process loading data final analysis steps.Open RStudioCreate new R Markdown document\nGive proper title select PDF format\nGive proper title select PDF formatGive document structure; e.g. create segment (using single hash #) called Exploratory Data Analysis.Whenever add code document make sure readable (keep maximum line length < 100 possible) sufficient documentation either text around code ‘chunks’ (preferred) using comments within code chunk.Read data file(s)\nPreferably use read.table function carefully set arguments. Open file text editor first check contents; header? can set row.names? columns needed? etc.\nNote: data set consists separate files (.e. one per sample) general tips reading data, see Appendix : Batch Loading Expression Data chapter.\nremainder document, try show either contents, structure - case - dimensions relevant R objects\nShow first five lines loaded data set. Including tables markdown document can done using pander function pander(Daróczi Tsegelskyi 2017) R-library.\nGive dimensions (dim() structure (str()) loaded data set.\nCheck output str function see columns expected R data type (e.g. values, factors, character, etc.)\n\nPreferably use read.table function carefully set arguments. Open file text editor first check contents; header? can set row.names? columns needed? etc.Note: data set consists separate files (.e. one per sample) general tips reading data, see Appendix : Batch Loading Expression Data chapter.remainder document, try show either contents, structure - case - dimensions relevant R objects\nShow first five lines loaded data set. Including tables markdown document can done using pander function pander(Daróczi Tsegelskyi 2017) R-library.\nGive dimensions (dim() structure (str()) loaded data set.\nCheck output str function see columns expected R data type (e.g. values, factors, character, etc.)\nShow first five lines loaded data set. Including tables markdown document can done using pander function pander(Daróczi Tsegelskyi 2017) R-library.Give dimensions (dim() structure (str()) loaded data set.Check output str function see columns expected R data type (e.g. values, factors, character, etc.)Examine samples included experiment create many R character objects needed store classification.\nNote: important item often overlooked. See last 4 lines code chunk section Example Data .\ninstance, eight samples divided case/ control columns create variable called case store column indices (numbers!) respective columns data set variable called control remaining four data column indices. later use.\nvariables allows us repeatedly access data rest course, .e. code boxplot(counts[case]) immediately shows plotting case samples counts data clear reading boxplot(counts[c(5:8)]) gives reference indices mean can lead typing mistakes.\nNote: important item often overlooked. See last 4 lines code chunk section Example Data .instance, eight samples divided case/ control columns create variable called case store column indices (numbers!) respective columns data set variable called control remaining four data column indices. later use.variables allows us repeatedly access data rest course, .e. code boxplot(counts[case]) immediately shows plotting case samples counts data clear reading boxplot(counts[c(5:8)]) gives reference indices mean can lead typing mistakes.\ndata sets order samples listed \nloaded data set different order shown GEO\nwebsite. often, names different lacking \ndescription non descriptive IDs. creating\nvariables control need point \ncontrol samples, need make sure right columns\ndata frame. Go Appendix A2;\nannotation order unclear just large number\nsamples might supporting data available can help\nmake sense sample layout.\nwant include external images log better control properties height width individual images can use following code. Note: code example need replace img_to_include.png file path actual image.):See Knitr chunk options page full listing chunk options.","code":"\n# Add the following to the code chunk header to control figure height, width and add a caption:\n#    {r, fig.height=5, fig.width=8, echo=FALSE, fig.cap=\"Image description\"}\n# The `dpi` argument is to remove any scaling\nknitr::include_graphics(\"figures/img_to_include.png\", dpi = NA)"},{"path":"EDA.html","id":"ch3_example_data","chapter":" Exploratory Data Analysis","heading":"3.2 Example Data","text":"section lists (publically available) data set(s) used chapter. chapter contains section new data sets used . Note examples, data different examples one challenges course translating examples data. Keep mind simple copy-pasting code fail reason. examples therefore print input data comparison data.section Normalization onwards, experiment identifier GSE101942 used visualizing normalized raw count values (earlier sections explore unnormalized data). experiment titled “Transcriptome analysis genetically matched human induced pluripotent stem cells disomic trisomic chromosome 21” experimental setup described follows:“12 total polyA selected samples. 6 IPSC samples 3 biological repeats trisomic samples 3 biological repeats disomic samples. 6 IPSC derived neuronal samples 3 biological repeats trisomic samples 3 biological repeats disomic samples.”following code shows load two available files containing raw-count data (one file per 6 samples), stored locally data/gse101942/ folder. simple ‘cleanup’ steps performed aren’t required helpful demonstration purposes (.e., renaming samples identify groups).GSE101942_IPSC_rawCounts.txt: raw count data induced pluripotent stem cells (iPSCs), trisomic disomicGSE101942_Neuron_rawCounts.txt: raw count data IPSC-derived Neurons, trisomic disomic. cells treated remove chromosome 21 iPSCs, described treatment protocol : “Targeted removal CHR 21 IPSC using TKNEO transgene”.See Table 1 Normalization section resulting contents counts dataframe.","code":"\n# Load the two sets of 6 samples\nipsc.data <- read.table('./data/gse101942/GSE101942_IPSC_rawCounts.txt',\n                        header = TRUE)\nneuron.data <- read.table('./data/gse101942/GSE101942_Neuron_rawCounts.txt',\n                          header=TRUE)\n# Merge by rowname (by = 0, see the help of 'merge')\ncounts <- merge(ipsc.data, neuron.data, by = 0, all.x = TRUE, all.y = TRUE,\n                sort = FALSE)\n# Change all NA's introduced by the merge to zeros\ncounts[is.na(counts)] <- 0\n# Show column names coming from the input files\nprint(names(counts))##  [1] \"Row.names\" \"c244A\"     \"c244B\"     \"c243\"      \"c2A\"       \"c2B\"      \n##  [7] \"c2C\"       \"C3_1\"      \"C3_2\"      \"C3_4\"      \"C2_2_1\"    \"C2_2_2\"   \n## [13] \"C2_2_3\"\n# Set the row names to the gene IDs stored in the 'Row.names' column\nrow.names(counts) <- counts$Row.names\n# Remove the gene ID column\ncounts <- counts[-1]\n# Rename samples to include their group name\nnames(counts) <- c(paste0('di_IPSC_r', 1:3),  # Disomic IPSC, replicates 1-3\n                   paste0('tri_IPSC_r', 1:3), # Trisomic IPSC, replicates 1-3\n                   paste0('di_NEUR_r', 1:3),  # Disomic Neuron, replicates 1-3\n                   paste0('tri_NEUR_r', 1:3)) # Trisomic Neuron, replicates 1-3\n# Show results of renaming the samples\nprint(names(counts))##  [1] \"di_IPSC_r1\"  \"di_IPSC_r2\"  \"di_IPSC_r3\"  \"tri_IPSC_r1\" \"tri_IPSC_r2\"\n##  [6] \"tri_IPSC_r3\" \"di_NEUR_r1\"  \"di_NEUR_r2\"  \"di_NEUR_r3\"  \"tri_NEUR_r1\"\n## [11] \"tri_NEUR_r2\" \"tri_NEUR_r3\"\n# Assign column -indices- to variables for later use (i.e. counts[di_NEUR] selects the relevant columns)\ndi_IPSC <- 1:3\ntri_IPSC <- 4:6\ndi_NEUR <- 7:9\ntri_NEUR <- 10:12"},{"path":"EDA.html","id":"EDA_part1","chapter":" Exploratory Data Analysis","heading":"3.3 Visualizing using boxplot and density plot","text":"segment describes basic steps visualization can performed Exploratory Data Analysis. start inspecting unnormalized count data. next steps normalization applied shown Normalization section . mentioned , focus EDA get overview/ perform bit Quality Control data set often requires visualizing data, figures need pretty. Simple figures perfectly fine stage. Try create figures data (keep log) add small description figure pointing anything different expect.Also note steps just selection. Furthermore, make sure every visualization make, add proper axis-labels containing measurement units (important!). might able see, values log-transformed using log2 function often numerical values high range ‘hide’ details plots. See section Fold Change value chapter 4 details. fine use non-log-transformed (simply raw-)data, otherwise use instance boxplot(log2(dataset)) plotting.\nInstead using basic R-plotting library\n(.e. plot, boxplot, hist, etc.)\ncan also opt using (challenging) ggplot2\nlibrary also used boxplot following section.\nconstructing ggplot2 plot feels like learning yet\nanother language, many resources available online can\nfollow.\n","code":""},{"path":"EDA.html","id":"eda_statistics","chapter":" Exploratory Data Analysis","heading":"3.3.1 Statistics","text":"Even basic statistics can give insight data performing 6-number-statistic data columns using summary() function. Note: can also use pander function pretty-print summary knitted output. notice look numbers produced executing function complete data set?","code":""},{"path":"EDA.html","id":"eda_boxplots","chapter":" Exploratory Data Analysis","heading":"3.3.2 Boxplots","text":"visual representation values can shown using boxplot. Boxplots easy create R data frame object just passing data columns. following boxplot shows data experiment separate boxplot sample. allows quick overview spotting irregularities (.e. checking replicates within sample-group show similar expression profiles). course, consider amount data single plot, can hint problems, need look much detail form quality control. Creating boxplot dataframe easy, saw using summary function; data large range maximum average values far apart. create lot outliers plot interesting later , boxplot can either:hide completely using outline = FALSE argument boxplot() (say figure description!), orperform log2() transformation can see .\ncreating plots expression data without performing\nlog-transformations see (instance boxplot) \nrange data large (0 1.000.000) lot\noutliers upper range. makes informative\nfigures, placing simple log2() function call around\ndata almost resolves issue.\n\nHowever, might introduce another problem (taking boxplot \nexample) log-transformation results data \ncontains negative infinity (-Inf) values. \ncaused 0-values data since\nlog2(0) == -Inf. circumvent issue, can add \npseudo count data simply adding value\n1 count values since log2(1) == 0;\nlog2(counts + 1).\n\nsituations boxplot adding \npseudo count value complete dataset can useful. Always add \npseudocount -place, meaning within plotting code instead\noverwriting original data set pseudo count added \nsteps require . cases (see density plot ) can useful \nclearly separate 0-values rest case \n0 < pseudo count < 1 value, 0.1\ncan used generates negative value (-3.3) \nfar away rest data. count values always\n>0 result positive log2 value.\n\nFigure 3.1: Boxplot comparing basic statistics genes across multiple samples\n","code":""},{"path":"EDA.html","id":"eda_density_plot","chapter":" Exploratory Data Analysis","heading":"3.3.3 Density Plots","text":"Another form visualizing data using density plot. method shows distribution (log2-transformed) count data samples allows easy spotting problems. plot commonly used analysing microarrays, still useful comparing complete dataset. code figure show example distribution 12 samples. things note figure huge peak exactly -3.321928 can ignored value calculated log2(0.1) explained box . peak consists 0-values (inactive genes) isn’t important us now. Therefore, added vertical line indicate left-part little interest. Note: lot code extra, simple inspection using line plotDensity(log2(counts + 0.1)) enough, rest extra example code.assess quality data look peaks plot. heigth peak doesn’t matter much, shifted peaks however indicate either lower higher amount reads sequenced sample. can confirmed first plot demonstrated next section library-sizes visualized.\nFigure 3.2: Density plot comparing count distribution 12 samples\n","code":"\n## The affy library has a density plotting function\nlibrary(affy)\n\n## Create a list of 4 colors to use which are the same used throughout this chapter \nlibrary(scales)\nmyColors <- hue_pal()(4)\n\n## Plot the log2-transformed data with a 0.1 pseudocount\nplotDensity(log2(counts + 0.1), col=rep(myColors, each=3),\n            lty=c(1:ncol(counts)), xlab='Log2(count)',\n            main='Expression Distribution')\n\n## Add a legend and vertical line\nlegend('topright', names(counts), lty=c(1:ncol(counts)),\n       col=rep(myColors, each=3))\nabline(v=-1.5, lwd=1, col='red', lty=2)"},{"path":"EDA.html","id":"EDA_part2","chapter":" Exploratory Data Analysis","heading":"3.4 Visualizing using heatmap and MDS","text":"section adds Exploratory Data Analysis techniques measure look distances samples based normalized data. Measuring distances two data objects (samples case) common task cluster analysis compare similarity (low distance indicates similar data). case calculate distances samples visualize heatmap using multidimensional scaling (MDS) technique.previous section used raw count data. might one samples different values (.e. shifted) compared samples. need raw count data use R packages edgeR (Chen et al. 2018) DESeq2 (Love, Anders, Huber 2017), calculating sample distances (used visualizations section) done form normalized data. data can either RPKM/FPKM/TPM/CPM vst-transformed (raw-)read counts. proper method transforming raw read count data using vst method DESeq2 R Bioconductor library shown . ‘variance stabilizing transformation’ normalized data used chapter, chapter 4 normalize using different technique.following code examples shows use library normalize count data GSE101942 experiment vst-normalized data calculate distance metric.new(“standardGeneric”, .Data = function (object, …)standardGeneric(“counts”), generic = structure(“counts”, package = “BiocGenerics”),package = \"BiocGenerics\", group = list(), valueClass = character(0),signature = \"object\", default = NULL, skeleton = (function (object,    ...)stop(gettextf(\"invalid call method dispatch '%s' (default method)\",look Table 1 immediately see huge difference groups genes. gene (A1BG) shows >100 expression tri_NEUR group compared di_IPSC group, within groups high variation indicating might actual expression different. One simple quick method inspection looking total number mapped reads per sample (sum sample/ column), sequencing-depth might vary across samples.Hint:Simple version plot can achieved barplot(colSums(counts) / 1e6)Table 2; Mapped reads per sample (millions) (continued )\nFigure 3.3: Library sizes samples GSE101942\nnumbers shown Table 2 barplot immediately clarify just experimental condition might caused large difference sample counts sequencing depth shows substantial difference . minimum maximum sequencing depth within trisomic neuron group (purple barplot). average expression A1BG gene within group 820, minimum 461 maximum 1010 lot variation. one reason telling us normalize data. following code chunks show DESeq2’s vst method (variance stabilizing transformation). Check Packages tab Rstudio see DESeq2 package installed load library command.","code":"package = \"BiocGenerics\", group = list(), valueClass = character(0),signature = \"object\", default = NULL, skeleton = (function (object,    ...)stop(gettextf(\"invalid call in method dispatch to '%s' (no default method)\","},{"path":"EDA.html","id":"normalization","chapter":" Exploratory Data Analysis","heading":"3.4.1 Normalization","text":"use vst function DESeq2 library must construct DESeqDataSet-object consisting count data combined sample annotation. Since want use (now) performing vst-transformation use basic form simple design sample names annotation (colData argument):now proper DESeqDataSet object can see , containing 56640 rows 12 columns (genes samples) gene symbols row names. Usually object hold data, requirement perform vst transformation good enough now. Next step performing transformation (can take depending size experiment results Large DESeqTransform object) retrieving actual data assay function object contains lot meta-data. table shows updated values now comparable across genes whereas raw count data harder compare.Table 2; VST transformed disomic IPSC count dataNote looking A1BG gene (one displayed figure 3.8), normalized expression values show much less variation across samples (shown, also within sample group). Therefore, assume large difference expression observed earlier might non-existent (might still significant difference though!).","code":"\n# Load the library\nlibrary('DESeq2')\n\n# DESeq2 will construct a SummarizedExperiment object and combine this \n# into a 'DESeqDataSet' object. The 'design' argument usually indicates the \n# experimental design using the condition(s) names as a 'factor', for now we use just '~ 1'\n(ddsMat <- DESeqDataSetFromMatrix(countData = counts,\n                                  colData = data.frame(samples = names(counts)),\n                                  design = ~ 1))## class: DESeqDataSet \n## dim: 56640 12 \n## metadata(1): version\n## assays(1): counts\n## rownames(56640): 5S_rRNA 7SK ... C1orf220 C2orf15\n## rowData names(0):\n## colnames(12): di_IPSC_r1 di_IPSC_r2 ... tri_NEUR_r2 tri_NEUR_r3\n## colData names(1): samples\n# Perform normalization\nrld.dds <- vst(ddsMat)\n# 'Extract' normalized values\nrld <- assay(rld.dds)"},{"path":"EDA.html","id":"eda_distance","chapter":" Exploratory Data Analysis","heading":"3.4.2 Distance Calculation","text":"now normalized data can use distance calculation. standard procedure many data analysis tasks calculates distance metric combination samples use check variation within sample groups. first need transpose matrix (rld) normalized values using t-function, dist function expects different samples rows genes columns. Note output matrix symmetric. Table 3 shows calculated distances within disomic IPSC group distance samples varies 100 140. maximum distance samples 380 demonstrated visualizations next sections.Table 3; Sample distances disomic IPSCs (Euclidean method)","code":"\n# Calculate basic distance metric (using euclidean distance, see '?dist')\nsampledists <- dist( t( rld ))"},{"path":"EDA.html","id":"eda_heatmap","chapter":" Exploratory Data Analysis","heading":"3.4.3 Sample Distances using a Heatmap","text":"(raw-)count data normalized format (TPM, RPKM, etc.), can optionally follow procedure count data create heatmap formats see makes difference. reason RNA-Seq method exists 10 years, still ongoing discussions subject data processing, especially regarding subjects like normalization technque use data analysis.following code block creates heatmap using pheatmap library offers many available heatmap functions. Using annotation data frame (can inspect contents ) identify samples based cell type ploidy. clustering shown heatmap clearly separates data based cell type differences ploidy seems minimal. Looking still, differences within single group minimal meaning - yet - inclined remove outlier-samples. Since 3 samples per category, also lose statistical power eventually remove one samples (also, always check article see remove samples prior data analysis).\nFigure 3.4: Heatmap showing Euclidean distances samples\nresulting heatmap shows interesting comparison across samples order samples purely determined using distance often different order data set. Ideally, want sample groups clustered together see data. particular heatmap clearly shows large difference overal expression induced pluripotent stem cells (iPSCs) differentiated cortical neurons difference - within cell type - disomic trisomic much lower.","code":"\n# We use the 'pheatmap' library (install with install.packages('pheatmap'))\nlibrary(pheatmap)\n\n# Convert the 'dist' object into a matrix for creating a heatmap\nsampleDistMatrix <- as.matrix(sampledists)\n\n# The annotation is an extra layer that will be plotted above the heatmap columns\n# indicating the cell type\nannotation <- data.frame(Cell = factor(rep(1:2, each = 6), \n                                          labels = c(\"IPSC\", \"Neuron\")),\n                         Ploidy = factor(rep(rep(1:2, each = 3), 2),\n                                         labels = c(\"disomic\", \"trisomic\")))\n                         \n# Set the rownames of the annotation dataframe to the sample names (required)\nrownames(annotation) <- names(counts)\n\npheatmap(sampleDistMatrix, show_colnames = FALSE,\n         annotation_col = annotation,\n         clustering_distance_rows = sampledists,\n         clustering_distance_cols = sampledists,\n         main = \"Euclidean Sample Distances\")"},{"path":"EDA.html","id":"eda_mds","chapter":" Exploratory Data Analysis","heading":"3.4.4 Multi-Dimensional Scaling","text":"following code example shows perform Multi-Dimensional Scaling (MDS) displays similarly calculated distances 2D-plot. experiment like two groups samples, hope see two clearly separated clusters formed, however ’ve seen heatmap, sample KO1B showed large deviation also see (confirm) using MDS.\nfigure plotted using ggplot2, \nadvanced method plotting R. plots preferred \nbase-R plotting, always sufficient use just can \nchallenging alter example code shown section. \ndata object plotted shown case contains simple X- \nY-coordinates.\nMDS use slightly different distance metric (Poisson instead default euclidean used dist function). following code shows calculate fitting distance sequencing data called Poisson Distance number optimizations available PoiClaClu library. library specifically designed handle read count data less influenced read-count differences across samples. Note , code usable count data function requires:n--p data matrix observations rows, p features columns. (,j) element x number reads observation mapped feature (e.g. gene exon) j.Table 5: MDS coordinates used plotting distances using Poisson Distance\nFigure 3.5: Multi Dimensional Scaling - Poisson Distance samples\nclearly see separation cell type less clear separation ploidy. spread within disomic groups much lower compared trisomic groups. Note y-axis range 4th x-axis range. sample di_NEUR_r3 shows large difference MDS plot. set genes identified differentially expressed repeat step expectation clear clustering.","code":"\nlibrary('PoiClaClu')\n# Note: uses the raw-count data, PoissonDistance performs normalization\n# set by the 'type' parameter (uses DESeq)\ndds <- assay(ddsMat)\npoisd <- PoissonDistance( t(dds), type = \"deseq\")\n# Extract the matrix with distances\nsamplePoisDistMatrix <- as.matrix(poisd$dd)\n# Calculate the MDS and get the X- and Y-coordinates\nmdsPoisData <- data.frame( cmdscale(samplePoisDistMatrix) )\n\n# And set some better readable names for the columns\nnames(mdsPoisData) <- c('x_coord', 'y_coord')\n# Separate the annotation factor (as the variable name is used as label)\ngroups <- factor(rep(1:4, each=3), \n                 labels = c(\"di_IPSC\", \"tri_IPSC\", \"di_NEUR\", \"tri_NEUR\"))\ncoldata <- names(counts)\n\n# Create the plot using ggplot\nggplot(mdsPoisData, aes(x_coord, y_coord, color = groups, label = coldata)) + \n  geom_text(size = 4) +\n  ggtitle('Multi Dimensional Scaling') +\n  labs(x = \"Poisson Distance\", y = \"Poisson Distance\") +\n  theme_bw()\n## Using base-R plot and legend functions:\n# plot(mdsPoisData, col=rep(myColors, each=3), pch=20, lwd=2)\n# legend(x=-20000, y=11000, legend = levels(groups), col=myColors, pch=20)"},{"path":"EDA.html","id":"cleaning-data","chapter":" Exploratory Data Analysis","heading":"3.5 Cleaning Data","text":"Conclude chapter cleaning dataset visualizations show necessity . Cleaning case means removing complete samples can classified outlier within group. tasks course, clear advice decide remove one samples. clear rule group must retain least three samples. one three samples visibly deviates replicates, even normalization, stay data set mentioned analysis report/ final article noise may introduced sample.\nNote: satisfied data reported findings, return steps next chapters. mentions creating heatmap instance, means heatmap expression data instead sample distances ’ve shown chapter.","code":""},{"path":"chapter-4.html","id":"chapter-4","chapter":" Discovering Differentialy Expressed Genes (DEGs)","heading":" Discovering Differentialy Expressed Genes (DEGs)","text":"first important ‘real’ analysis step finding genes show difference expression sample groups; differentially expressed genes (DEGs). concept might sound rather simple; calculate ratios genes samples determine fold-change (FC) denoting factor change expression groups. , filter genes actually show difference. give list genes showing different behaviour across samples, need focus genes show difference, observed difference also statistically significant!determine whether gene can classified significant DEG going use existing libraries, one article might refer . ’re going brief manual examination count data make sure properly understand ’re looking : significance Fold-Change value.","code":""},{"path":"chapter-4.html","id":"deg-preprocessing","chapter":" Discovering Differentialy Expressed Genes (DEGs)","heading":"4.1 Pre-processing","text":"Given results exploratory data analysis performed chapter 3, might concluded one samples show () deviating expression patterns compared samples group. mentioned , enough (> 3) samples group, might opt remove sample reduce noise statistical tests sensitive . Since performing analysis steps programatically also easy test DEGs without sample(s) question see removal results lower p-values (= higher significance). Remember always properly document choice reasoning!perform manual steps data using libraries, re-normalization different technique. simple method normalizing count data calculate fragments per million mapped fragments (FPM, equal CPM) value transform log2. Opposed FPKM RPKM include gene-length calculation (likely don’t ). dataset also includes FPKM RPKM (pre-normalized), allowed use data .FPM CPM simply dividing count values sample total number reads sample divided 1 million shown :Another step can - might guided article - filter (partially) inactive genes. required using edgeR library (DESeq2 includes low-count filtering methods), good exercise manually use calculating lFC values. data sets available contains lot 0-measurements (see instance density plot previous chapter); genes reads mapped. R-studio Environment tab, click data set (perform View(data) command click one sample columns order data ascending. now likely (unless bacterial data instance) see lot zero values columns.experiment two groups, three replicates , three six samples 0-reads mapped, often advised remove gene completely. can subjective experiment, might expected (thus important) comparing different tissues knockout experiments. Also, genes () low read count (< 5) can give high (artificial) FC value (see left hand side MA-plot). Comparing two samples one value 2 11, reads -regulated gene factor 5.5 might actually just noise!Assignment: Search article advice filter zero values low count genes. nothing stated subject, think tactic (search literature/ online!). perfectly fine discuss peers.\ndescribe aspect; filter data, clearly explain (likely article stated proper reason)\nPerform filtering data set. likely need use one apply, rowSums, rowMeans functions, combined maybe , functions.\nManually verify rows removed correctly filtered.\nProperly document many genes filtered !\ndescribe aspect; filter data, clearly explain (likely article stated proper reason)Perform filtering data set. likely need use one apply, rowSums, rowMeans functions, combined maybe , functions.Manually verify rows removed correctly filtered.Properly document many genes filtered !","code":"\n# Perform a naive FPM normalization\n# Note: log transformation includes a pseudocount of 1\ncounts.fpm <- log2( (counts / (colSums(counts) / 1e6)) + 1 )"},{"path":"chapter-4.html","id":"deg-fc","chapter":" Discovering Differentialy Expressed Genes (DEGs)","heading":"4.2 The Fold Change Value","text":"FC usually given calculated log2 case/control ratio. example, gene average expression 30 mapped reads control group 88 reads experiment group, ratio case/control 2.93. Ratio values > 1 indicate increased expression experiment relation control values 0 1 indicate lower expression. log2 transformed value ratio calculated log2(2.93) results 1.55. counts reversed, ratio 30/88, 0.34. previously calculated value 2.93 means 3-fold regulation 0.34 value means 3-fold regulation can see range numbers different. Comparing log2 values 1.55 -1.55 compares much better.easy calculate FC genes , simple FC value doesn’t mean much, yet! still need use power lie within replicates sample group. Using replicates, want determine observed FC just biological noise sequencing error.Assignment: Create histogram log2-FC values genes. Apply following steps FPM-normalized data least one comparison (need possible combinations):\naverage replicates two groups separately add new columns FPM-normalized data set,\ncalculate log Fold Change (lFC); simply subtract one group (averaged-column) another, usually experiment-control\nNote since already log-transformed data, calculate ratio (experiment/control). Also, Logarithm quotient rule, divide average columns ratio calculation, subtract: \\(log_b{(x/y)} = log_b{(x)} - log_b{(y)}\\)\n\nplot data using hist(logFC_column, breaks=60) function (change breaks argument needed),\nAdd two vertical lines -1 1 (using abline(v=...)) indicate significance (2-fold change).\ndata shows low fold-change values might see genes outside boundary. mean find significantly changed genes (DEGs).\n\naverage replicates two groups separately add new columns FPM-normalized data set,calculate log Fold Change (lFC); simply subtract one group (averaged-column) another, usually experiment-control\nNote since already log-transformed data, calculate ratio (experiment/control). Also, Logarithm quotient rule, divide average columns ratio calculation, subtract: \\(log_b{(x/y)} = log_b{(x)} - log_b{(y)}\\)\nNote since already log-transformed data, calculate ratio (experiment/control). Also, Logarithm quotient rule, divide average columns ratio calculation, subtract: \\(log_b{(x/y)} = log_b{(x)} - log_b{(y)}\\)plot data using hist(logFC_column, breaks=60) function (change breaks argument needed),Add two vertical lines -1 1 (using abline(v=...)) indicate significance (2-fold change).\ndata shows low fold-change values might see genes outside boundary. mean find significantly changed genes (DEGs).\ndata shows low fold-change values might see genes outside boundary. mean find significantly changed genes (DEGs).calculated FC shows large change expression groups means nothing variation within group high. reason use form statistical test checks variation group difference groups samples. simplest form usually comes using t-test:“can used determine two sets data significantly different ”often instead manually performing t-tests ANOVA analysis, article mentions use one R libraries software packages finding DEGs. exactly also chapter. output finding DEG’s always includes - limited - list p-values; one p-value per gene. value indicates whether gene statistically-significant differentially expressed gene (SDEG) find genes interest need get genes p-value threshold (.e. 0.05).","code":""},{"path":"chapter-4.html","id":"deg-bioconductor","chapter":" Discovering Differentialy Expressed Genes (DEGs)","heading":"4.3 Using Bioconductor Packages","text":"section demonstrates use two packages perform DEG-analysis count data. many packages available Bioconductor RNA-Seq analysis, DSS, EBSeq, NOISeq BaySeq, focus edgeR DESeq2 processing count-based data. Chances one two packages mentioned article described use R statistical analysis widely accepted methods processing gene expression data.packages apply normalization methods (described sections ) therefore work raw count data. data consists expected count data, need round values (using floor function). can choose one two packages use, use since increase statistical power hard use understand model data ’ll show next.Assignment: Perform data-analysis using least one two packages. however article mentioned different R package analysis (instance limma package) allowed use one instead. non-R packages can done extra (dataset can always analyzed using edgeR DESeq2). Optionally, use one package can compare results packages.\nStart documenting relevant information article, also used software analysis.\nStart documenting relevant information article, also used software analysis.two listed packages come complete extensive manual. Therefore, course manual (document) guide exactly use ; find . manuals can found PDF BioConductor website, opened within R using vignette function. see available manuals (note; different function-help documentation) execute function without arguments: vignette(). lists vignettes edgeR DESeq2 :vignette(\"DESeq2\") (whole course partly designed document!)vignette(\"edgeR\")\n","code":""},{"path":"chapter-4.html","id":"deg-design","chapter":" Discovering Differentialy Expressed Genes (DEGs)","heading":"4.3.1 The Design (matrix)","text":"packages need properly specify samples grouped. seen examples using R factor object heatmap MDS visualizations tell groups samples group sample belongs. Reading documentation packages shows important part performing DE-analysis. example, following code shown edgeR documentation page 8 two sample groups defined (numbered 1 2), placed factor object used input model.matrix function.can often complex part seems theoretical, practice ’s difficult. really grasp topic strongly adviced watch two videos define design:StatQuest YouTube: Using Linear Models t-tests ANOVA, Clearly ExplainedStatQuest YouTube: Design Matrices Linear Models, Clearly ExplainedThese videos expects basic understanding topics already discussed statistics course, otherwise mentioned topic videos videos available explaining . See video-index handy overview topics discussed channel.mentioned need think question want ask; difference want know? two sample groups shown question rather easy; ‘genes show effect case/ control samples?’. two sample groups however question becomes complex. experiment comparing influence three kinds drugs (thus three groups) combined effect time, want focus influence drugs time? valid questions define way create design matrix.Documentation subject plenty, however often contains overwhelming information. page contains valuable details (can safely start reading Choice design), including following text based example used edgeR documentation:“examples cover , use linear models make comparisons different groups. Hence, design matrices ultimately work least two columns: intercept column, consists column 1’s, second column, specifies samples second group. case, two coefficients fit linear model: intercept, represents population average first group, second coefficient, represents difference population averages second group first group. latter typically coefficient interested performing statistical tests: want know difference two groups.”two sample groups want change question (.e. test influence different group), read section Releveling factor.","code":"\ngroup <- factor(c(1,1,2,2), labels = c(\"Control\", \"Case\"))\n(design <- model.matrix( ~ group))##   (Intercept) groupCase\n## 1           1         0\n## 2           1         0\n## 3           1         1\n## 4           1         1\n## attr(,\"assign\")\n## [1] 0 1\n## attr(,\"contrasts\")\n## attr(,\"contrasts\")$group\n## [1] \"contr.treatment\""},{"path":"chapter-4.html","id":"deg-deseq2","chapter":" Discovering Differentialy Expressed Genes (DEGs)","heading":"4.3.2 DESeq2","text":"used DESeq2 library DEG analysis re-use DESeqDataSet object adviced create new object proper design formula instead ~ 1 used .need normalize data using previously used vst function DESeq2 library normalize count data follows:“DESeq computes scaling factor given sample computing median ratio, gene, read count geometric mean across samples. uses assumption genes DE uses median ratios obtain scaling factor associated sample.”example data four sample groups re-use annotation data frame created EDA: Heatmap section (printed ). allows us comparisons based celltype ploidy samples. One minor note importance first level factor taken reference, expression analysis, control group. example (see table ) IPSC control-celltype disomic control-ploidy. simple case/ control study, first level (check levels() function) control group. Use relevel() function levels need switched.proper DESeqDataSet need run DESeq function object. , using results function can extract DEGs DESeqResults object. Applying summary function object(s) shows number regulated DEGs can seen lists impact trisomy IPSC cell type (~800 DEGs) comparing IPSC Neuron disomic (control) cell types (can seen, huge amount genes affected differentiation stem cell cortical neurons):output results function contains following columns gene:call results function depends heavily experiment. can see output summary function, details given comparison shown (also, default p-value 0.1 used instead 0.05).Depending design used create DESeqDataSet , one comparisons can made (applying DESeq function calculates filter results function). Read help results function carefully; especially regarding contrast argument define comparison retrieve. following code can used example experiment get DEGs comparing celltypes using adjusted p-value 0.05:DESeq2 library contains number plotting functions can applied DESeqResults object (output results function), notable plotMA function. Note first MA-plot shown high range log fold changes (-10, 10) maximum value 22.4 (shown triangle stating outside plotting range). log fold change 22 means >5 million increased expression seems artificially high. DESeq2 includes function perform downstream processing estimated log fold change values called lfcShrink advised always run afterwards. reason executing function described vignette :“useful visualize MA-plot shrunken log2 fold changes, remove noise associated log2 fold changes low count genes without requiring arbitrary filtering thresholds.”\nFigure 4.1: MA-plots standard DESeq2 output (left) shrinking ‘lfcShrink’ (right)\nLinksAnalyzing RNA-seq data DESeq2 comprehensive guide analyzing RNA-Seq data using DESeq2 (part document used material !). adviced read first sections guide take good look index document many interesting sections might help later.Publication, accompanying article showing differences performance compared methods packages.","code":"\nout of 44504 with nonzero total read count\nadjusted p-value < 0.05\nLFC > 0 (up)       : 432, 0.97%\nLFC < 0 (down)     : 417, 0.94%\noutliers [1]       : 219, 0.49%\nlow counts [2]     : 17004, 38%\n(mean count < 6)\n[1] see 'cooksCutoff' argument of ?results\n[2] see 'independentFiltering' argument of ?results\n\n\nout of 44504 with nonzero total read count\nadjusted p-value < 0.05\nLFC > 0 (up)       : 8224, 18%\nLFC < 0 (down)     : 9594, 22%\noutliers [1]       : 219, 0.49%\nlow counts [2]     : 11906, 27%\n(mean count < 2)\n[1] see 'cooksCutoff' argument of ?results\n[2] see 'independentFiltering' argument of ?results\n\n\nres <- results(dds, contrast = c(\"Cell\", \"IPSC\", \"Neuron\"),\n               alpha = 0.05)"},{"path":"chapter-4.html","id":"deg-edger","chapter":" Discovering Differentialy Expressed Genes (DEGs)","heading":"4.3.3 edgeR","text":"One mature libraries RNA-Seq data analysis edgeR library available Bioconductor. complete (sometimes bit complex) manual available need read Chapter 2 focus 2.1 2.7, 2.9 - complex design - 2.10. Section 1.4 (Quickstart) shows code example steps needed DEG analysis using count data using two glm methods (quasi-likelihood F-tests likelihood ratio tests). steps shown identical non-glm method calculating fit object can replaced performing exactTest function shown section 2.9.2.edgeR library normalize count data follows:“trimmed means M values (TMM) Robinson Oshlack, implemented edgeR, computes scaling factor two experiments using weighted average subset genes excluding genes exhibit high average read counts genes large differences expression.”Running edgeR requires raw count data together grouping-factor packaged DGEList object (DGEList() function). Furthermore, proper model.matrix object (see section design) needed input estimateDisp function. exact steps take (variations DESeq2) must searched documentation linked .analysis done can retrieve actual results topTags function:Quitting lines 696-762 [.edger-TEST] (gene_expression.Rmd)\nError h():\n! error evaluating argument ‘x’ selecting method function ‘.data.frame’: object ‘qlf’ found\nBacktrace:\n1. pander::pander(.data.frame(topTags(qlf)), caption = “significant genes given edgeR”)\n4. base::.handleSimpleError(<fn>, “object ‘qlf’ found”, base::quote(topTags(qlf)))\n5. base (local) h(simpleError(msg, call))package also contains plotting methods can use intermediate steps analysis. instance, calculating normalization factors (calcNormFactors), can perform multi-dimensional scaling plotMDS function similar ’ve done manually chapter three:\nFigure 4.2: edgeR MDS plot based calculated log2 fold changes\ndispersion running estimateDisp function plotBCV function:\nFigure 4.3: edgeR plot several dispersion methods\nlog-fold changes genes, output exactTest function (output et DGEExact object) plotSmear function. abline shows log-FC threshold:LinksDifferential Expression Analysis using edgeR tutorialAnother tutorial hosted GitHub","code":"\nplotMDS(y)\nplotBCV(y)\ndeGenes <- decideTestsDGE(qlf, p=0.05)\ndeGenes <- rownames(qlf)[as.logical(deGenes)]\nplotSmear(qlf, de.tags=deGenes)\nabline(h=c(-1, 1), col=2)"},{"path":"chapter-5.html","id":"chapter-5","chapter":" Data Analysis and Visualization","heading":" Data Analysis and Visualization","text":"point one sets DEGs experiment. might just one set might large number comparisons made set DEGs per comparison / results multiple approaches (manual, packages). chapter offers guidance visualize, summarize - finally - connect biological relevance results. Browse chapter first choosing visualization(s) make data! article contains relevant figures like try reproduce (given enough data available) can done well instead visualizations shown . Note however Volcano plot described required.","code":""},{"path":"chapter-5.html","id":"volcano","chapter":" Data Analysis and Visualization","heading":"5.1 Volcano Plot","text":"volcano plot often first visualization data statistical tests completed. plot shows data genes highlight genes considered DEG using thresholds (adjusted) p-value fold-change. Many articles describe values used thresholds methods section, otherwise good default 0.05 adjusted p-value around 1.0 log-FC.use EnhancedVolcano library creates ggplot2 volcano plot shown . packaged function deseq.volcano just create two figures need change input data set plot title. resulting scatterplot places -log10(pvalue) values y-axis log-FC x-axis. often results plot representing volcanic ‘eruption’ fold-change influences spread significance height. Coloring done based thresholds (-log10(0.05) adjusted log2(2) log-FC) follows:(dark) gray: significant; FC adjusted p-value thresholds.green: high FC, adjusted p-value threshold.blue: adjusted p-value indicates significance, FC thresholdred: significant DEGYou can adjust pCutoff FCcutoff values depending requirements, see elaborate documentation examples GitHub. code wraps EnhancedVolcano function custom function (deseq.volcano) can used multiple comparisons make.Notes:used edgeR instead need change names columns (.e. log2FoldChange == logFC padj == FDR) used example.example many DEGs annotate name plots. EnhancedVolcano plot function many options annotating significant genes, see linked manual information.","code":"\nlibrary(EnhancedVolcano)\n\n## Simple function for plotting a Volcano plot, returns a ggplot object\ndeseq.volcano <- function(res, datasetName) {\n  return(EnhancedVolcano(res, x = 'log2FoldChange', y = 'padj',\n                         lab=rownames(res),\n                         title = paste(datasetName, \"trisomic vs disomic\"),\n                         subtitle = bquote(italic('FDR <= 0.05 and absolute FC >= 2')),\n                         # Change text and icon sizes\n                         labSize = 3, pointSize = 1.5, axisLabSize=10, titleLabSize=12,\n                         subtitleLabSize=8, captionLabSize=10,\n                         # Disable legend\n                         legendPosition = \"none\",\n                         # Set cutoffs\n                         pCutoff = 0.05, FCcutoff = 2))\n}\n\n## Note: input data is the corrected DESeq2 output using the 'lfcShrink' function (see chapter 4)\ndeseq.volcano(res = res.ipsc.lfc, datasetName = \"IPSC\")\ndeseq.volcano(res = res.neuron.lfc, datasetName = \"Neuron\")\n## Note that the most simplest case is:\nEnhancedVolcano(res.ipsc, x = 'log2FoldChange', y = 'padj', lab = rownames(res.ipsc))"},{"path":"chapter-5.html","id":"venn-diagram","chapter":" Data Analysis and Visualization","heading":"5.2 Venn Diagram","text":"Another common visualization Venn-diagram. data set ’ve shown two comparisons; trisomic vs disomic two cell types. Using Venn-diagram can show shared unique DEGS IPSC trisomic vs IPSC disomic NEUR trisomic vs NEUR disomic comparisons.Creating Venn-diagram can good method selecting (sub)set DEGs, instance selecting set DEGs shared groups (Intersection, see Wikipedia page Set Theory way around, DEGs unique group (Difference). Keep mind continuing clustering section likely useful cluster combinations might .Code examples:Data preparation (first code-chunk): creating vectors DEGs packages. DEGs gene names (taken row names DEG-analysis).\noutput DESeq2 results lfcShrink functions stored res.ipsc.lfc res.ipsc.neuron objects.\noutput DESeq2 results lfcShrink functions stored res.ipsc.lfc res.ipsc.neuron objects.Using VennDiagram library (second code-chunk): good complex/ large number comparison, highly configurable difficult use. see examples, run following function R console: example('draw.quad.venn')Using gplots library (third code-chunk): Easier use support complex Venn diagrams, good adjusting appearance.See blog post “VENN DIAGRAM R RSTUDIO: MILLION WAYS” methods creating Venn diagrams.two created vectors gene names can used either taking lengths (VennDiagram library) just names (gplots library).\nFigure 5.1: Venn diagram comparing IPSC NEUR DEGs using ‘VennDiagram’ library\nmentioned, following easier alternative (especially > 2 groups) venn function gplots library. needs list containing either row-numbers gene-names DEGs easier (offers less adjustability make prettier).\nFigure 5.2: Venn diagram comparing IPSC NEUR DEGs using ‘gplots’ library\nanalysis result similar amount DEGS (849 IPSC 769 Neuron cell types) 96 DEGs shared cell types thus unique.","code":"\n## Data preparation\npval_threshold <- 0.05\nipsc.degs <- row.names(res.ipsc.lfc[which(res.ipsc.lfc$padj <= pval_threshold), ])\nneur.degs <- row.names(res.neuron.lfc[which(res.neuron.lfc$padj <= pval_threshold), ])\n## Venn-diagram using the `VennDiagram` library (see below for alternative method)\nlibrary(VennDiagram)\n# Arguments for a pairwise (two-sets) venn-diagram are sizes for set1, set2 and overlap (intersect)\n# Many more functions are available for triple, quad and quantuple diagrams (starting with 'draw.***')\nvenn.plot <- draw.pairwise.venn(length(ipsc.degs),\n                                length(neur.degs),\n                                # Calculate the intersection of the two sets\n                                length( intersect(ipsc.degs, neur.degs) ),\n                                category = c(\"IPSC Trisomic vs Disomic\", \"NEURON Trisomic vs Disomic\"), scaled = F,\n                                fill = c(\"light blue\", \"pink\"), alpha = rep(0.5, 2),\n                                cat.pos = c(0, 0))\n\n# Actually plot the plot\ngrid.draw(venn.plot)\n## Alternative method using the `gplots` library\nlibrary(gplots)\n# Create a Venn-diagram given just the list of gene-names for both sets\nvenn(list(\"IPSC Trisomic vs Disomic\" = ipsc.degs,\n          \"NEURON Trisomic vs Disomic\" = neur.degs), )"},{"path":"chapter-5.html","id":"clustering","chapter":" Data Analysis and Visualization","heading":"5.3 Clustering","text":", can create heatmap found DEGs samples colors show fold-change value. Usually heatmap shown two color gradient; red (downregulated) green (upregulated). however optional using default colors pheatmap example perfectly fine long proper legend present.multiple methods creating heatmap (one ’ve already used) directly apply clustering visualization. clustering can applied expression pattern gene (row-clustering), expression pattern sample (column-clustering) (default pheatmap).following example commonly found publications shows log2 Fold Change values comparison done. first fifth columns show log2(FC) calculated DESeq2 two trisomic vs disomic comparisons. remaining columns log fold changes separate trisomic replicates vs average disomic expression (manually calculating Log(FC) done section 4.1 using DESeq2 normalized count data rld.dds object). DESeq2 calculated FC values show greater changes compared manually calculated FC values. value showing manually calculated values see replicates show patterns within groups. however optional unless one comparison create heatmap single column.Note column-clustering turned type plot setting cluster_cols = FALSE (replicates group, clustering add useful information).\nFigure 5.3: Heatmap showing fold changes Trisomic vs Disomic comparisons cell types. shown genes shared DEGS (96 shown Venn-diagrams) also filtered absolute log2(FC) value >2.\nlimitation heatmap show 100 DEGs found; just doesn’t fit well single figure causes clustering especially gene names/ labels unreadable. therefore always adviced filter adjusted p-value, also log-FC value reduce number genes shown. Alternatively data properly annotated, instead showing gene names row, showing GO-term clustering might reveal expression patterns certain gene groups.","code":""},{"path":"chapter-5.html","id":"pathway-analysis","chapter":" Data Analysis and Visualization","heading":"5.4 Pathway Analysis","text":"cases ’ll see hundreds even thousands DEGs result analysis. large amounts DEGs much visualization approaches easily say something biological context many genes.One approach tackling large set DEGs pathway analysis different methods genes grouped pathway get overview affected pathways experiment.section demonstrate two methods analysis, one using online platform gene-annotation enrichment analysis R-method signaling pathway impact analysis.","code":""},{"path":"chapter-5.html","id":"kegg-pathway-visualization-10-degs","chapter":" Data Analysis and Visualization","heading":"5.4.1 KEGG-pathway Visualization (>10 DEGs)","text":"Knowing interesting KEGG-pathway(s) upfront (.e. listed article) allows visualizations applied selected pathway. can use pathview website Bioconductor library:“… maps renders wide variety biological data relevant pathway graphs. users need supply data specify target pathway. Pathview automatically downloads pathway graph data, parses data file, maps user data pathway, render pathway graph mapped data.”section demonstrates use bioconductor package visualizing DEGs onto pathways, also change expression gene.Table continues belowSYMBOL, GENENAME, ENSEMBL, ENSEMBLPROT, UNIGENE, UNIPROT, ACCNUM, ENSEMBLTRANS, REFSEQ, ENZYME, TAIR, PROSITE ORF","code":"\n## Bioconductor library (install with bioclite() if missing)\nlibrary(pathview)\n\n## Example pathway IDs (for human, change organism key for other organisms)\ndata(\"paths.hsa\")\npander(head(paths.hsa, n=5))\n## Check 'gene.idtype' argument possibilities\ndata(gene.idtype.list); \npander(gene.idtype.list)\n## Prepare data for visualization\ndeseq.degs.logfc <- subset(deseq.results, padj < pval_threshold, select = log2FoldChange)\n\npander(head(deseq.degs.logfc))\npathview(gene.data=deseq.degs.logfc,\n         pathway.id=\"00020\", # TCA-cycle\n         species=\"hsa\"     # Organism key\n)"},{"path":"datasets.html","id":"datasets","chapter":"A Discovering RNA-Seq Data Sets","heading":"A Discovering RNA-Seq Data Sets","text":"chapter describes method finding public RNA-Seq data sets interest. data set context course refers data belonging certain gene-expression experiment, usually consisting number sequencing-samples combined meta-data describing experiment.","code":""},{"path":"datasets.html","id":"finding-public-data","chapter":"A Discovering RNA-Seq Data Sets","heading":"A.1 Finding Public Data Sets of Interest","text":"number large databases online offer access thousands - published - experiments focus searching downloading gene expression data high-throughput (NGS) sequencing techniques (RNA-Seq).thousands freely available data sets possible start performing research without need actually performing lab-experiment. common conditions, organisms tissues can download samples compare multiple studies find novel relations genes conditions verify experiments performed laboratory. instance, power public data sets demonstrated jointly three alumni article called Calling genotypes public RNA-sequencing data enables identification genetic variants affect gene-expression levels.(al. 2015).’s method section begins sentence “downloaded raw reads available human RNA-seq data sets …” amount data work become clear later . interested expression Quantitative Trait Loci (eQTLs) Allele-Specific Expression (ASE) analysis, please read interesting paper.start diving large complex databases potentially containing thousands experiments millions samples terabytes data, need get rough idea like found something interest order know looking first place. definitive guide protocol can followed processing analyzing gene-expression data set, following goals can considered:Re-(part ) analysis described accompanying publication.\ndata sets - except ones added recently - accompanied publication, can gather lot information regarding experimental design results just paper. often, researchers already performed interesting research data set therefore good exercise try reproduce results. Challenges getting proper understanding techniques used researchers (article) translating analysis steps. Often though, details clear article requires interpretation results performing analysis steps arrive conclusion.\ndata sets - except ones added recently - accompanied publication, can gather lot information regarding experimental design results just paper. often, researchers already performed interesting research data set therefore good exercise try reproduce results. Challenges getting proper understanding techniques used researchers (article) translating analysis steps. Often though, details clear article requires interpretation results performing analysis steps arrive conclusion.Alternatively can come ideas (biological) question(s) try answer given data set instead redoing published work.\ncourse challenging suitable data sets.\nmight notice publications focus small set genes instead whole transcriptome. instance, perform experiment yeast want measure activity genes involved alcohol fermentation, researchers might look genes Glycolysis pathway. leaves 6000+ yeast genes explore possibly come novel relations gene expression experimental setup.\nanalysis approach explained paper focusing pre-selected list genes, depending experiment might come comparisons done original research. instance, experiment maternal age (birth) correlated autism offspring, paternal age also known addressed publication. subject research (spoiler alert; clear conclusion can drawn including extra factor…).\ncourse challenging suitable data sets.might notice publications focus small set genes instead whole transcriptome. instance, perform experiment yeast want measure activity genes involved alcohol fermentation, researchers might look genes Glycolysis pathway. leaves 6000+ yeast genes explore possibly come novel relations gene expression experimental setup.analysis approach explained paper focusing pre-selected list genes, depending experiment might come comparisons done original research. instance, experiment maternal age (birth) correlated autism offspring, paternal age also known addressed publication. subject research (spoiler alert; clear conclusion can drawn including extra factor…).Another type project evaluate different methods either normalization statistical analysis either confirm published results find novel genes involved experiment.\ncan viewed technical research subject ok pursue, however final report mainly concern biological impact findings.\nmeans accompanying publication shows conclusive results, probably come acknowledging results (therefore, extend point 1 list). conclusive result hope find results case might become hard formulate good conclusion end.\nHowever, publication isn’t specific states conclusion ‘Benzene exposure shows increased risk leukemia’ followed list genes might involved, try see can find genes might involved changing analysis approach.\n\ncan viewed technical research subject ok pursue, however final report mainly concern biological impact findings.\nmeans accompanying publication shows conclusive results, probably come acknowledging results (therefore, extend point 1 list). conclusive result hope find results case might become hard formulate good conclusion end.\nHowever, publication isn’t specific states conclusion ‘Benzene exposure shows increased risk leukemia’ followed list genes might involved, try see can find genes might involved changing analysis approach.\nmeans accompanying publication shows conclusive results, probably come acknowledging results (therefore, extend point 1 list). conclusive result hope find results case might become hard formulate good conclusion end.However, publication isn’t specific states conclusion ‘Benzene exposure shows increased risk leukemia’ followed list genes might involved, try see can find genes might involved changing analysis approach.Although last project goal many risks involved need sound project proposal, might result interesting project. shown linked article, possible combine data multiple experiments.\ninstance find two related experiments (.e. researching effect certain drug) measuring expression different tissue types (.e. liver brain tissue). analyzing experiments, present set genes show effect tissues - interestingly - genes showing effect one tissue.\ninstance find two related experiments (.e. researching effect certain drug) measuring expression different tissue types (.e. liver brain tissue). analyzing experiments, present set genes show effect tissues - interestingly - genes showing effect one tissue.list guideline kept mind browsing suitable data sets extended week three. Also, don’t worry terms unclear, getting good understanding used terminology part first weeks","code":""},{"path":"datasets.html","id":"ncbi-geo","chapter":"A Discovering RNA-Seq Data Sets","heading":"A.2 The NCBI Gene Expression Omnibus","text":"easiest method exploring finding interesting data sets simply using web-browser access data-repository. One repository NCBI Gene Expression Omnibus (GEO). repository primarily used store micro array data sets describe experiments, linking raw data, processed data accompanying publication. however many sources data browsable GEO, project limit searches experiments using RNA-Seq technology 60.000 experiments listed. Following “Expression profiling high throughput sequencing” link summary table get complete listing can find interesting experiment. Besides finding experiment interested , important requirements need account browsing listed . Note point need download data yet, described chapter 3.experiment must published publication must fully available, free otherwise Hanze University library.sample group, minimum three replicates must present. one groups less three replicates, group used (rest data might still usable minimum two groups 3 replicates).\nReplicates often easy spot. lucky samples named “WildType - rep 2” instance often need read “Overall design” part page click samples (“GSM…” ID) make sure.\nReplicates often easy spot. lucky samples named “WildType - rep 2” instance often need read “Overall design” part page click samples (“GSM…” ID) make sure.available data must contain least count data.\nData offered multiple formats, always including RAW data (actual reads don’t use), also processed data FPKM, RPKM TPM (see interesting blog-post discussing use data). tools R analysis libraries use downstream analysis rely unnormalized unprocessed data count data. counts simply represent number reads mapped transcript.\noften hardest part assess; sometimes clear provided data actual raw count values. Usually, can find (sometimes rather cryptic) information clicking sample ID ‘Data processing’ section. paragraph text often mention ‘Supplementary_files_format_and_content’ describes actual content downloadable files. Following number examples valid invalid data provided experiment. single-sample pages can find description. Unfortunately, authors experiment free explain data please, likely encounter exact descriptions:\nValid: Global effects SUPT4H1 RNAi gene expression HEK293 cells; description: Supplementary_files_format_and_content: tab-delimited, gzip-compressed text file matrix raw quantitation results sample. Mentions raw-anything likely valid.\nValid: Transcriptome profiling normal transformed mouse urothelial cells; description: Supplementary_files_format_and_content: raw_counts_unfiltered.txt.gz: Raw counts.. Even though 4 different file types available download, one contains actual raw count data.\nValid: Quantitative Analysis PPARD RNA-Seq Transcriptomes Mouse Gastric Corpus Epithenial Cells Next Generation Sequencing (NGS); description: Supplementary_files_format_and_content: raw counts sequencing reads genes. Single file download containing raw counts.\nValid: Impact shRNA-mediated KLF4 -modulation transcriptome profile human keratinocyte precursor cells; description: Supplementary_files_format_and_content: raw normalized counts.\nValid: Role TAF5L TAF6L controlling mouse ESC state; description: Supplementary_files_format_and_content: Comma-separated text files containing raw counts, normalized counts, FPKM expression values TPM expression values.. common include normalized raw values.\nInvalid: Co-chaperone Mzb1 key effector Blimp1 plasma cell differentiation b1 integrin function; description: Supplementary_files_format_and_content: normalized gene expression values (FPKM) calculated Cuffnorm (v2.2.1) provided TXT format. clearly states normalized data (FPKM values) available.\nInvalid: EGFR-targeted therapy-induced resistance mechanism glioblastoma; description: Supplementary_files_format_and_content: .txt file including gene name, RefSeq id, locus, RPKM value gene sample. RPKM values available normalized raw count values, usable.\nInvalid: RNA-seq subcellular fractions nESC EpiLC; description: Supplementary_files_format_and_content: _genes.fpkm_table.txt: Tab-delimited text files include FPKM values replicate.*. FPKM data available.\n\nData offered multiple formats, always including RAW data (actual reads don’t use), also processed data FPKM, RPKM TPM (see interesting blog-post discussing use data). tools R analysis libraries use downstream analysis rely unnormalized unprocessed data count data. counts simply represent number reads mapped transcript.often hardest part assess; sometimes clear provided data actual raw count values. Usually, can find (sometimes rather cryptic) information clicking sample ID ‘Data processing’ section. paragraph text often mention ‘Supplementary_files_format_and_content’ describes actual content downloadable files. Following number examples valid invalid data provided experiment. single-sample pages can find description. Unfortunately, authors experiment free explain data please, likely encounter exact descriptions:\nValid: Global effects SUPT4H1 RNAi gene expression HEK293 cells; description: Supplementary_files_format_and_content: tab-delimited, gzip-compressed text file matrix raw quantitation results sample. Mentions raw-anything likely valid.\nValid: Transcriptome profiling normal transformed mouse urothelial cells; description: Supplementary_files_format_and_content: raw_counts_unfiltered.txt.gz: Raw counts.. Even though 4 different file types available download, one contains actual raw count data.\nValid: Quantitative Analysis PPARD RNA-Seq Transcriptomes Mouse Gastric Corpus Epithenial Cells Next Generation Sequencing (NGS); description: Supplementary_files_format_and_content: raw counts sequencing reads genes. Single file download containing raw counts.\nValid: Impact shRNA-mediated KLF4 -modulation transcriptome profile human keratinocyte precursor cells; description: Supplementary_files_format_and_content: raw normalized counts.\nValid: Role TAF5L TAF6L controlling mouse ESC state; description: Supplementary_files_format_and_content: Comma-separated text files containing raw counts, normalized counts, FPKM expression values TPM expression values.. common include normalized raw values.\nInvalid: Co-chaperone Mzb1 key effector Blimp1 plasma cell differentiation b1 integrin function; description: Supplementary_files_format_and_content: normalized gene expression values (FPKM) calculated Cuffnorm (v2.2.1) provided TXT format. clearly states normalized data (FPKM values) available.\nInvalid: EGFR-targeted therapy-induced resistance mechanism glioblastoma; description: Supplementary_files_format_and_content: .txt file including gene name, RefSeq id, locus, RPKM value gene sample. RPKM values available normalized raw count values, usable.\nInvalid: RNA-seq subcellular fractions nESC EpiLC; description: Supplementary_files_format_and_content: _genes.fpkm_table.txt: Tab-delimited text files include FPKM values replicate.*. FPKM data available.\nValid: Global effects SUPT4H1 RNAi gene expression HEK293 cells; description: Supplementary_files_format_and_content: tab-delimited, gzip-compressed text file matrix raw quantitation results sample. Mentions raw-anything likely valid.Valid: Transcriptome profiling normal transformed mouse urothelial cells; description: Supplementary_files_format_and_content: raw_counts_unfiltered.txt.gz: Raw counts.. Even though 4 different file types available download, one contains actual raw count data.Valid: Quantitative Analysis PPARD RNA-Seq Transcriptomes Mouse Gastric Corpus Epithenial Cells Next Generation Sequencing (NGS); description: Supplementary_files_format_and_content: raw counts sequencing reads genes. Single file download containing raw counts.Valid: Impact shRNA-mediated KLF4 -modulation transcriptome profile human keratinocyte precursor cells; description: Supplementary_files_format_and_content: raw normalized counts.Valid: Role TAF5L TAF6L controlling mouse ESC state; description: Supplementary_files_format_and_content: Comma-separated text files containing raw counts, normalized counts, FPKM expression values TPM expression values.. common include normalized raw values.Invalid: Co-chaperone Mzb1 key effector Blimp1 plasma cell differentiation b1 integrin function; description: Supplementary_files_format_and_content: normalized gene expression values (FPKM) calculated Cuffnorm (v2.2.1) provided TXT format. clearly states normalized data (FPKM values) available.Invalid: EGFR-targeted therapy-induced resistance mechanism glioblastoma; description: Supplementary_files_format_and_content: .txt file including gene name, RefSeq id, locus, RPKM value gene sample. RPKM values available normalized raw count values, usable.Invalid: RNA-seq subcellular fractions nESC EpiLC; description: Supplementary_files_format_and_content: _genes.fpkm_table.txt: Tab-delimited text files include FPKM values replicate.*. FPKM data available.","code":""},{"path":"a2-data_loading.html","id":"a2-data_loading","chapter":"B Batch Loading Expression Data in R","heading":"B Batch Loading Expression Data in R","text":"code example shows batch-load multiple files containing expression (count) data single sample. data example can found GEO ID GSE109798. Depending tool used convert mapping data per-sample count file, contents files can () different. Downloading data example experiment GEO gives us single .tar file called GSE109798_RAW.tar. Extracting archive file nets us folder following files:","code":"\nfile.names <- list.files('./data/GSE109798_RAW/')"},{"path":"a2-data_loading.html","id":"decompressing","chapter":"B Batch Loading Expression Data in R","heading":"B.0.1 Decompressing","text":"file extension files .txt.gz means files compressed using gzip need unpacked can loaded. easiest method using system gunzip command files can done within R applying gunzip command using system function file.Now can update file.names variable since file name changed.","code":"\n## Change directory to where the files are stored\nsetwd('./data/GSE109798_RAW/')\nsapply(file.names, FUN = function(file.name) {\n  system(paste(\"gunzip\", file.name)) \n  })\nfile.names <- list.files('./data/GSE109798_RAW/')"},{"path":"a2-data_loading.html","id":"determining-data-format","chapter":"B Batch Loading Expression Data in R","heading":"B.0.2 Determining Data Format","text":"Next, can inspect contents files, assuming layout/ column names etc. decide need use analysis.files contain (much) just count value gene can see columns (transcript) length, TPM, FPKM, etc. Also, count-column called expected_count raises questions well.\nexpected count value usable contains \ninformation - compared raw count - actually require. \nexpected part results multimapped reads \nsingle read mapped multiple positions genome. \ntranscript originates one location, multimapped read \nusually discarded. expected count though, instead \ndiscarding read completely estimated \noriginates added fraction count\nvalue. value 1.32 see line 5 \nexample means true count 1 (uniquely mapped\nread) .32 (estimated part) results \nalgorithm can mean multiple things.\n\nmentioned , require integer count data use \npackages DESeq2 edgeR \ntwo methods convert expected count raw count data: +\nround value nearest integer (widely accepted method\nwell within expected sampling variation), + discard \nfraction part using example floor()\nfunction.\n","code":"\n## Call the system 'head' tool to 'peek' inside the file\nsystem(paste0(\"head \", \"./data/GSE109798_RAW/\", file.names[1]))"},{"path":"a2-data_loading.html","id":"loading-data-1","chapter":"B Batch Loading Expression Data in R","heading":"B.0.3 Loading Data","text":"columns want keep transcript_id expected_count columns ignore rest (might interested data later though). need lead file separately can define function reads data, keeping columns interest returning dataframe data. Note first line file used header, check setting header argument TRUE, sometimes expression data starts line 1. file name also used name column dataframe know column sample. done splitting file name (using strsplit) using dot (‘.’) keeping first part (.e. ‘GSM2970156_HCC1806E224A’) discarding second part (‘isoforms.results.txt’). strsplit function however always returns list, case containing vector 5 splitted elements:Applying read_sample function file names gives us set data frames can merge together using merge function. merge data based transcript id defined = 1 argument pointing first column. start reading just one file ‘base’ dataframe merge files.processing seemed data set divided two groups also listed GEO website project:GPL11154 Illumina HiSeq 2000 (Homo sapiens)GPL13112 Illumina HiSeq 2000 (Mus musculus)first 6 files human source last 6 mouse. Therefore, following code shows read first 6 samples merge single dataframe. Repeating process 6 files result another dataframe samples.Table continues belowThe dataset variable now contains data first 6 samples experiment. advisable compare number rows data set number rows single sample. guaranteed samples exactly number genes/transcripts present (.e., 0-values might discarded) results final data set many rows smallest sample. See help merge case argument can used introduce extra rows missing data.","code":"\n## String splitting in R\n## (the fixed = TRUE is required as the dot is a special character, see '?strsplit')\nstrsplit('GSM2970155_HCC1806E224B.isoforms.results.txt.gz', '.', fixed = TRUE)## [[1]]\n## [1] \"GSM2970155_HCC1806E224B\" \"isoforms\"               \n## [3] \"results\"                 \"txt\"                    \n## [5] \"gz\"\n## Keeping the sample identifier\nstrsplit('GSM2970155_HCC1806E224B.isoforms.results.txt.gz', '.', fixed = TRUE)[[1]][1]## [1] \"GSM2970155_HCC1806E224B\"\n## Function for reading in files\nread_sample <- function(file.name) {\n  ## Extract the sample name for naming the column\n  sample.name <- strsplit(file.name, \".\", fixed = TRUE)[[1]][1]\n  ## Read the data, setting the 'transcript_id' as row.names (column 1)\n  sample <- read.table(file.name, header = TRUE, sep=\"\\t\", row.names = NULL)\n  ## Rename the count column\n  names(sample)[5] <- sample.name\n  ## Return a subset containing the 'transcript_id' and sample name columns\n  return(sample[c(1, 5)])\n}\nsetwd('./data/GSE109798_RAW/')\n\n## Read the FIRST sample\ndataset <- read_sample(file.names[1])\n\n## Read first sample group (6)\nfor (file.name in file.names[2:6]) {\n  sample <- read_sample(file.name)\n  dataset <- merge(dataset, sample, by = 1)\n}\n\npander(head(dataset))"},{"path":"a2-annotation.html","id":"a2-annotation","chapter":"C Annotating an RNA-Seq Experiment","heading":"C Annotating an RNA-Seq Experiment","text":"chapter describes annotate data, meaning assigning names functions Differentially Expressed Genes. step can either easy bit challenging depending data source project. first example relevant data set originated NCBI GEO, ‘’ always annotated default. found data set elsewhere (GEO identifier project), skip Manually annotate data section .goal chapter - least - find Gene Symbol common ID (NCBI/ Ensembl) gene. Using information much easier find relevant information online sources say something functionality impact DEGs.Even though might lucky either already data simple GEO query results everything want/ need, tasked also perform manual annotation phase even use found annotation data. written code access sources good excercise using complex Bioconductor packages. Also, can always use find information require later .","code":""},{"path":"a2-annotation.html","id":"download-geo","chapter":"C Annotating an RNA-Seq Experiment","heading":"C.1 Downloading annotation data from GEO","text":"Given GSE**** id experiment can download available (annotation)data GEO using GEOquery library. example experiment used previously, following code downloads two files; series_matrix.txt.gz file also available download NCBI GEO website. Note destdir argument storing known location instead /tmp folder later use.downloaded, can load data using:gse object may contain useful information samples genes. Unfortunately example data set information available regarding samples can access phenoData slot. Slots R objects method storing multiple data objects single R object, case object called ExpressionSet (see ?ExpressionSet description contents structure). storing information experiment handy; can store data frame actual expression data, list containing laboratory information, chunk text abstract, etc. single object.common data sets found ExpressionSet object:\n* @assayData: actual expression data. common microarray experiments () rare RNA-seq experiments\n* @phenoData: sample information, present GEO data set\n* @featureData: data frame holding feature (=gene) information, slightly rare RNA-seq data sets\n* @experimentData: information lab performed experimentYou already know access column data frame using $ notation, accessing complete data frame gse object done using @ symbol:see sample described using sample identifier (GSM****), information sample group (title column) number characteristics. actual information included depends data set.Contents gse@phenoData@data information (selection) (continued )data set properly annotated gene information, accessible @featureData slot. instance, data set GSE20489 annotation available. loaded, see following annotation columns present gene featureData slot (following small subset, total can well 20 columns information per gene). Also, metadata listing shows many features (=genes, 54675) samples (54) expression data present (@assayData slot):Gene information stored ExpressionSet object experiment GSE20489 (continued )","code":"\nlibrary(GEOquery)\ngetGEO(GEO = 'GSE101942', destdir = \"./data/\")\nlibrary(GEOquery)\ngse <- getGEO(filename = './data/GSE101942_series_matrix.txt.gz')\n# Access a data frame containing phenotype information\ngse@phenoData@data\n# This gets a list with a single 'ExpressionSet' object\nGSE20489 <- getGEO(GEO = \"GSE20489\", destdir = './data')\n# Print some metadata of the featureData slot\nGSE20489[[1]]## ExpressionSet (storageMode: lockedEnvironment)\n## assayData: 54675 features, 54 samples \n##   element names: exprs \n## protocolData: none\n## phenoData\n##   sampleNames: GSM514737 GSM514738 ... GSM514790 (54 total)\n##   varLabels: title geo_accession ... tissue:ch1 (39 total)\n##   varMetadata: labelDescription\n## featureData\n##   featureNames: 1007_s_at 1053_at ... AFFX-TrpnX-M_at (54675 total)\n##   fvarLabels: ID GB_ACC ... Gene Ontology Molecular Function (16 total)\n##   fvarMetadata: Column Description labelDescription\n## experimentData: use 'experimentData(object)'\n##   pubMedIds: 23883607 \n## Annotation: GPL570"},{"path":"a2-annotation.html","id":"manually-annotate","chapter":"C Annotating an RNA-Seq Experiment","heading":"C.2 Manual Data Annotation","text":"Luckily us, R offers number libraries automatically retrieve information large sets genes, unless chosen organism much data available. section demonstrates use two libraries, starting AnnotationDbi followed biomaRt.AnnotationDbi library downloads local copy organism-specific database gene information biomaRt uses online databases retrieve data given query. biomaRt offers far data (1000 data fields per organism) complex use.Since biomaRt also relying online databases might good strategy annotate genes interest (DEGs) instead querying >20.000 genes might retain 20 statistical analysis. planning use biomaRt, skip Discovering Differentialy Expressed Genes (DEGs) chapter first return Using R Bioconductors biomaRt section annotate data.","code":""},{"path":"a2-annotation.html","id":"annotationDBI","chapter":"C Annotating an RNA-Seq Experiment","heading":"C.2.1 Using AnnotationDBI","text":"AnnotationDbi offers data sets many organisms form installable libraries depending experiment need find proper library. example experiment contains samples house mouse (Mus musculus) therefore select data set.following information types available database (use columns function inspect).Available fields databaseThe table shows data available information types given randomly chosen EntrezID ‘1080’. Note table split show 24 data types values.Example values field ‘org.Hs.eg.db’Retrieving data locally stored annotation database can done using mapIds function takes number arguments:x: local database querykeys: IDs data set, often ENSEMBL IDs gene SYMBOLS. example data rownames data set contain gene SYMBOLS (see Table 1 page 2)column: data column retrieve databasekeytype: type data provide. example data gene SYMBOLmultiVals: number columns contain one entry single gene, case want store first one.output mapIds function single character vector can add data frame containing DEGs (case DEGs IPSC trisomic vs disomic comparison containing 828 genes stored ipsc.tri.vs.di data frame).Unfortunately, organisms offer access information information always complete. instance, following table shows number available records Ensembl, Entrez en KEGG IDs downloaded:Statistics annotation columns using AnnotationDbiFor data sets might luck, otherwise continue biomaRt method explained .","code":"\n# Load the AnnotationDbi interface library\nlibrary(AnnotationDbi)\n# Load the Bioconductor installation library (contains 'biocLite()')\nlibrary(BiocInstaller)\n\n# Install and load the organism specific gene database\n# 'org' for Organism\n# 'Hs' for Homo sapiens\n# 'eg' for Entrez Gene IDs\n# (try to load (using the library function) first before installing, it might already be present)\nbiocLite('org.Hs.eg.db')\nlibrary(org.Hs.eg.db)\n.Hs.data.table <- cbind(.Hs.data[1:12,], .Hs.data[13:24,])\ncolnames(.Hs.data.table) <- c('Data Type', 'Example Value',\n                              'Data Type', 'Example Value')\npander(.Hs.data.table, caption=\"Example values for each field in 'org.Hs.eg.db'\")\n# Retrieve the ENSEMBL gene ID, e.g. 'ENSG00000001626'\nipsc.tri.vs.di$Ensembl <- mapIds(x = org.Hs.eg.db,\n                           keys=row.names(ipsc.tri.vs.di),\n                           column=\"ENSEMBL\",\n                           keytype=\"SYMBOL\",\n                           multiVals=\"first\")\n\n# Retrieve the ENTREZ gene ID, e.g. '12323'\nipsc.tri.vs.di$EntrezID <- mapIds(x = org.Hs.eg.db,\n                            keys=row.names(ipsc.tri.vs.di),\n                            column=\"ENTREZID\",\n                            keytype=\"SYMBOL\",\n                            multiVals=\"first\")\n\n# Retrieve the KEGG enzyme code, e.g. 2.7.11.17\nipsc.tri.vs.di$Enzyme <- mapIds(x = org.Hs.eg.db,\n                           keys=row.names(ipsc.tri.vs.di),\n                           column=\"ENZYME\",\n                           keytype=\"SYMBOL\",\n                           multiVals=\"first\")"},{"path":"a2-annotation.html","id":"using-biomart","chapter":"C Annotating an RNA-Seq Experiment","heading":"C.2.2 Using biomaRt","text":"short explanation biomart (source: Wikipedia):“purpose BioMarts Ensembl Genomes allow user mine download tables containing genes single species, genes specific region chromosome genes one region chromosome associated InterPro domain. BioMarts also include filters refine data extracted attributes (Variant ID, Chromosome name, Ensembl ID, location, etc.) appear final table file can selected user.”text mentiones species, attributes filters, need combine elements query biomart databases annotation. following code ‘chunks’ show possible values elements gather store relevant data.biomaRt library interfaces biomart.org online database. Sometimes biomart website (offers browsable database) due maintenance, many mirrors can still used, instance Ensembl.required explore number objects given example code likely suitable data/ organism. project well documented Bioconductor biomaRt website links biomaRt users guide.Using listDatasets function can get full list available datasets. Store list R object ‘browse’ object RStudio see organism included. Copy name dataset, species element use.Subset 214 datasets availableNext deciding filter. can use listFilters function ensembl object, storing full list filters. wise view RStudio find filter use.filter specifies use search . instance, AnnotationDbi queries gave us Ensembl gene ID’s use ensembl_gene_id filter. viewing list filters RStudio can use search text-box top-right view, example ensembl look filters use type identifier.Subset 437 available filtersFinally get point determine data like get database. attributes can get listAttributes function ensembl object. - especially attributes since often >1000 selectable options - store attributes view RStudio look data want.Subset 3172 available attributesIf want gene chromosome, start- end-position well description (note, just example, , interesting information available !) combine character vector (attrs.get, see ). one caveat though, order get back results input order! means simple combine data original data set need merge together. However, able merge data need know record belongs gene therefore add selected filter (case ensembl_gene_id) list attributes get. , get results dataframe can use merge function R combine gene information actual data:now three needed elements (species, filter attributes) ready query database getBM function (read help using ?getBM).example retrieves data using set five Ensembl gene IDs since genes Ensembl ID ’ve seen , filter first use values parameter .results object data.frame 5 columns can merge data set giving us following annotation columns (combined AnnotationDBI biomaRt libraries).just example use biomaRt library comes selecting correct filter looking interesting attributes retrieve. information can found documentation avaialble `vignette(‘biomaRt’)","code":"\n# Load the library\nlibrary(biomaRt)\n\n# Use an alternative database server as the regular one sometimes has issues..\nensembl=useMart(\"ENSEMBL_MART_ENSEMBL\", host=\"https://www.ensembl.org\")\n# Select the 'ensembl' database\nensembl <- useMart(\"ensembl\")\nmart.datasets <- listDatasets(ensembl)\n\n# Select the correct dataset, for the example data we select the 'hsapiens_gene_ensembl' \nensembl <- useDataset('hsapiens_gene_ensembl', mart = ensembl)\nfilters <- listFilters(ensembl)\nattributes <- listAttributes(ensembl)\nmerge(x = ipsc.tri.vs.di, y = results, by.x = 'Ensembl', by.y = 'ensembl_gene_id')\n# Set the 'attributes' values\nattrs.get <- c(\"ensembl_gene_id\", \"chromosome_name\", \n               \"start_position\",\"end_position\", \"description\")\n\n# Perform a biomaRt query using 'getBM'\nresults <- getBM(attributes = attrs.get,\n                 filters = \"ensembl_gene_id\",\n                 values = ipsc.tri.vs.di$Ensembl[1:5], \n                 mart = ensembl)\n\nresults$gene_length <- abs(results$end_position - results$start_position)"}]
+[{"path":"index.html","id":"chapter-1","chapter":" Introduction","heading":" Introduction","text":"Course 2.1.2 involves two separate projects performed single, published, data set. First, start performing Genomics analysis similar case study module 2.1.1. , take result mapping step, process read counts perform subsequent Transcriptomics analysis.","code":""},{"path":"index.html","id":"genomics","chapter":" Introduction","heading":"1.1 2.1.2 Genomics","text":"Genomics analysis going look variants chosen data reference genome. One important aspect analysis often assessing quality, input data set well output certain steps. set variants determined can annotated finally linked biology underlying original experiment, either reproducing findings aiming novel ideas.","code":""},{"path":"index.html","id":"project-deliverables","chapter":" Introduction","heading":"1.1.1 Project Deliverables","text":"Project Proposal: short written proposal project clearly describing data set, research question methods + tools used.Genomics pipeline: genomics analysis pipeline well documented one lab journals contains correct tools based research question data set.Presentation: research results presented attractive method proper level given target audience.","code":""},{"path":"index.html","id":"data-processing","chapter":" Introduction","heading":"1.1.2 Data Processing","text":"Similar 2.1.1 NGS & Mapping course go similar steps process data (.e. quality control, trimming, mapping, variant calling, etc.). time however, using workflow manager Galaxy, execute tools command line. Two reasons number size input files far larger compared previously used data sets Galaxy server capable handling many large data sets. Also, proper terminal experience valuable.Note writing pipeline: steps need documented RMarkdown documents. step must introduced short description tool used, version followed command used execute tool. terminal command, use code cells bash language. code cells must option eval=FALSE set. prevent code executed knitting document.Furthermore, also required run command RStudio server provided https://bioinf.nl/rstudio. RStudio server hosted assemblix server capable accessing data running tools large data sets. workstations provide classrooms capable tasks. Also, server already number tools pre-installed (FastQC, trimmomatic, BWA/STAR mappers, etc.). Lastly, running RMarkdown document RStudio server forces document commands compared simply using ssh.Note storing (intermediate) data: data relevant project must stored /students/2024-2025/Thema05/ directory. , create project directory team create sub-directories go. home folder hard limit 25GB data therefore sufficient storing project files. Use commands chmod chown allow project partners access data .Note processing large amounts data: building correct terminal commands execute instance trimming step large number files can quite challenging. therefore required test-data used steps. Reasons :Spotting errors quickly instead wait several hours see error messageInspecting resource usage (CPU, memory, disk space) prevent overloading system\nsteps can memory intensive (1 TB large data sets!) can crash system\nsteps can memory intensive (1 TB large data sets!) can crash systemCreates small output set can used test set next stepNote use GATK: Genome Analysis Toolkit (GATK) often used set tools variant calling. Many publications cite use GATK unfortunately include needed details reproduce results. Due complexity GATK recommend use, unless publication clear exact workflow settings used. Prepare spend lot time official GATK website browsing manuals can use effectively.","code":""},{"path":"index.html","id":"finding-a-data-set","chapter":" Introduction","heading":"1.1.3 Finding a Data Set","text":"Finding suitable data set can challenging requirements data set needs meet:needs DNA- RNA-Seq data set available.data set needs published (publication open access)RNA-Seq part, minimum amount samples required\nleast two sample groups minimum 3 samples per group (replicates)\nSee appendix chapter A1 example valid experimental setups, chapter also contains tips (mostly relevant Transcriptomics).\nleast two sample groups minimum 3 samples per group (replicates)See appendix chapter A1 example valid experimental setups, chapter also contains tips (mostly relevant Transcriptomics).Data available NCBI Sequence Read Archive Gene Expression Omnibus.proper way search data set suitable Genomics Transcriptomics projects creating query searches : (“expression profiling high throughput sequencing”[DataSet Type]) “genome variation profiling high throughput sequencing”[DataSet Type]. selection actually contains anything interest ; find something interesting please use different search strategy. One search RNA-Seq experiments GEO one groups preferably mutant.Note genomics part, available samples processed. transcriptomics, common case/control experimental setup rather single condition group possibly replicates. data set include many samples, recommended select subset samples process (discuss reasoning behind selection).","code":""},{"path":"index.html","id":"transcriptomics","chapter":" Introduction","heading":"1.2 2.1.2 Transcriptomics","text":"use RNA-Sequencing techniques measuring gene expression relatively new replaces microarrays, though cases microarrays still used. Gene expression data gives valuable insights workings cells certain conditions. Especially comparing instance healthy diseased samples can become clear genes causal influence specific condition. Finding genes interest (genes showing differing expression accross conditions, called Differentially Expressed Genes (DEGs)) goal project.golden standard analyzing RNA-sequencing data sets many tools (manufacturers sequencing equipment also deliver software packages) use proven R libraries processing, visualizing analyzing publicly available data sets. manual describes steps often performed transcriptomics experiment. Use general guideline process data accompanying chosen published research. Note manual originally written processing count data therefore describe processing NGS read files count data part Genomics analysis performed course. RNA-Seq counts describe number reads mapped gene corresponds relative number transcripts (mRNA sequences) gene present cell time sampling.","code":""},{"path":"index.html","id":"project-deliverables-1","chapter":" Introduction","heading":"1.2.1 Project Deliverables","text":"Project Proposal: short written proposal project clearly describing data set, research question methods + libraries used.Lab Journal: lab journal readable, complete aimed reproducability research. Relevant visualisations include statistical methods performed support biological interpretation results.","code":""},{"path":"index.html","id":"lab-journal","chapter":" Introduction","heading":"1.3 Lab Journal","text":"may know previous projects, essential keep proper journal detailing every step done analysis. journal kept R markdown file, showing steps taken analysis data set. markdown knitted single PDF-file project completed thus containing text detailing steps decisions ’ve made, R-code (always visible!) resulting output/ images.\nNotes: general advice; wait \nknitting whole document project done knitting \nprone errors trying fix large document \neasy. Give code chunk proper attributes, including name \nminimum. helps spot errors knitting process\nmentiones chunk processed. Note chunk-names\nmust unique. try make proper use chapters \nsections include (optional) table contents.\n","code":""},{"path":"index.html","id":"schedule","chapter":" Introduction","heading":"1.4 Schedule","text":"Genomics Transcriptomics parts schedule listed ","code":""},{"path":"index.html","id":"genomics-1","chapter":" Introduction","heading":"1.4.1 Genomics","text":"Empty rows means previous topic task likely take one session.","code":""},{"path":"index.html","id":"transcriptomics-1","chapter":" Introduction","heading":"1.4.2 Transcriptomics","text":"","code":""},{"path":"index.html","id":"literature","chapter":" Introduction","heading":"1.5 Literature","text":"following (online) documents can used throughout course:Gene Set Enrichment Analysis – Theory. Available : https://www.pathwaycommons.org/guide/primers/data_analysis/gsea/Li D. Statistical Methods RNA Sequencing Data Analysis. : Husi H, editor. Computational Biology [Internet]. Brisbane (AU): Codon Publications; 2019 Nov 21. Chapter 6. Available : https://www.ncbi.nlm.nih.gov/books/NBK550334/Yu, L., Fernandez, S. & Brock, G. Power analysis RNA-Seq differential expression studies. BMC Bioinformatics 18, 234 (2017). https://doi.org/10.1186/s12859-017-1648-2Sun, S., Hood, M., Scott, L., Peng, Q., Mukherjee, S., Tung, J., & Zhou, X. (2017). Differential expression analysis RNAseq using Poisson mixed models. Nucleic acids research, 45(11), e106. https://doi.org/10.1093/nar/gkx204","code":""},{"path":"from-reads-to-counts.html","id":"from-reads-to-counts","chapter":" From Reads to Counts","heading":" From Reads to Counts","text":"Step one RNA-Seq data analysis generating gene transcript counts sample. overlaps partly steps taken genomics workflow:Download SRA files extract FastQ filesPerform quality control raw reads (FastQC)Perform trimming needed (Trimmomatic)Perform quality control post trimming (FastQC)Perform mapping reference genome (often STAR)Generate gene/ transcript counts (either STAR, HTseq featurecounts)Downstream analysis (chapter 4 document)chapter shows steps example data set (GSE149995)","code":""},{"path":"from-reads-to-counts.html","id":"downloading-and-extracting-fastq-files","chapter":" From Reads to Counts","heading":"2.1 Downloading and extracting FastQ files","text":"Given SRA identifier project, select samples download accession file (SRX identifiers). Using identifiers input prefetch system command (requires NCBI SRA-Toolkit installed) downloads SRA files convert using fasterq-dump command fastq files.fasterq-dump default performs split 3 operation, see SRA toolkit manual details.","code":"prefetch $(<GSE149995_SRA_ACC_RNA-Seq.csv) --output-directory \\\n  /local-fs/staff/marcel/GSE149995/SRAfind /local-fs/staff/marcel/GSE149995/SRA/ -name \"*.sra\" | \\\n  parallel 'fasterq-dump -O /local-fs/staff/marcel/GSE149995/fastq/ {}'"},{"path":"from-reads-to-counts.html","id":"optional-annotation","chapter":" From Reads to Counts","heading":"2.2 (Optional) Annotation","text":"NCBI GEO page lists 9 samples 3 conditions (WT, potent arf7). However, experiment includes 18 sequence files. downloading ‘series matrix file’ GEO (using GEOquery library) combining run info downloadable SRA, can list mapping run sample.Note: optional step now useful doubts group/sample setup.","code":"\nlibrary(GEOquery)\ngse <- getGEO(GEO = \"GSE149995\")\n# Combine the condition with a sample number and replicate number\nCondition = paste(rep(gse[[1]]@phenoData@data$`genotype/variation:ch1`, each = 2),\n                  rep(1:3, each = 2),\n                  paste0('r', rep(1:2, 9)), sep = '_')\nrun_info <- read.csv(file = \"data/GSE149995_Sra_RunInfo.csv\")\nsetup <- cbind(run_info %>% dplyr::select(Run, Experiment), Condition)\npander(setup)"},{"path":"from-reads-to-counts.html","id":"quality-control","chapter":" From Reads to Counts","heading":"2.3 Quality Control","text":"perform FastQC -> Trimmomatic -> FastQC steps","code":""},{"path":"from-reads-to-counts.html","id":"fastqc","chapter":" From Reads to Counts","heading":"2.3.1 FastQC","text":"Running FastQCQC-report using multiqc","code":"\nfqdir <- \"/local-fs/staff/marcel/GSE149995/fastq/\"\nfastqc(fq.dir = \"/local-fs/staff/marcel/GSE149995/fastq/\",\n       qc.dir = \"/local-fs/staff/marcel/GSE149995/fqc\",\n       fastqc.path = \"/usr/bin/fastqc\",\n       threads = 32)multiqc \\\n  --filename ~/Development/2.1.2-Transcriptomics/GSE149995_FQC_multiqc_report.html \\\n  /local-fs/staff/marcel/GSE149995/fqc/ "},{"path":"from-reads-to-counts.html","id":"trimming","chapter":" From Reads to Counts","heading":"2.4 Trimming","text":"Proper trimming requires experimentation research regarding used sequencing technology. following command uses Trimmomatic (release 0.39, available Github) trim fastq files. Note example single-end sequencing , use TrimmomaticPE command paired-end data processing TruSeq FASTA file ending PE instead.Re-running FastQC","code":"cat data/GSE149995_Sra_RunInfo.csv | \\\n  parallel 'TrimmomaticSE -threads 4 ' \\\n                  '/local-fs/staff/marcel/GSE149995/fastq/{}.fastq.gz ' \\\n                  '/local-fs/staff/marcel/GSE149995/fastq-trimmed/{}.trimmed.fastq.gz ' \\\n                  'ILLUMINACLIP:/homes/marcelk/Development/2.1.2-Transcriptomics/TruSeq3-SE.fa:2:30:10 ' \\\n                  'MINLEN:40 ' \\\n                  'SLIDINGWINDOW:4:20'\nfastqc(fq.dir = \"/local-fs/staff/marcel/GSE149995/fastq-trimmed/\",\n       qc.dir = \"/local-fs/staff/marcel/GSE149995/fqc-trimmed\",\n       fastqc.path = \"/usr/bin/fastqc\",\n       threads = 32)multiqc --filename \\\n  ~/Development/2.1.2-Transcriptomics/GSE149995_FQC_Trimmed.html \\\n  /local-fs/staff/marcel/GSE149995/fqc-trimmed/"},{"path":"from-reads-to-counts.html","id":"mapping-with-star","chapter":" From Reads to Counts","heading":"2.5 Mapping with STAR","text":"","code":""},{"path":"from-reads-to-counts.html","id":"building-star-index","chapter":" From Reads to Counts","heading":"2.5.1 Building STAR index","text":"bwa, STAR requires extra annotation file (GFF GTF format) describing features genome. data set requires FASTA Arabidopsis thaliana genome accompanying annotation file. genomeDir specifies output location index.","code":"cd /local-fs/staff/marcel/GSE149995/reference\nmkdir star; cd star\nSTAR \\\n  --runThreadN 12 \\\n  --runMode genomeGenerate \\\n  --genomeDir /local-fs/staff/marcel/GSE149995/reference/star \\\n  --genomeFastaFiles /local-fs/staff/marcel/GSE149995/reference/TAIR10_chr_all.fas \\\n  --sjdbGTFfile /local-fs/staff/marcel/GSE149995/reference/TAIR10_GTF_genes.gtf \\\n  --genomeSAindexNbases 12 \\\n  --sjdbOverhang 49"},{"path":"from-reads-to-counts.html","id":"read-mapping-and-transcript-counting","chapter":" From Reads to Counts","heading":"2.5.2 Read mapping and transcript counting","text":"following command performs mapping reference genome, saves sorted BAM file directly performs read count transcripts (NOTE: requires --sjdbGTFfile parameter index building). input file, ReadsPerGenes..tab file generated containing 4 columns:column 1: gene IDcolumn 2: counts unstranded RNA-seqcolumn 3: counts 1st read strand aligned RNA (htseq-count option -s yes)column 4: counts 2nd read strand aligned RNA (htseq-count option -s reverse)See interesting discussion BioStars regarding column use.planning use different tool calculating gene counts, can remove quantMode argument generate sorted BAM files. Please read manual tool choice see indexed BAM file BAM file index.Generating MultiQC report STAR mapping phase","code":"find /local-fs/staff/marcel/GSE149995/fastq-trimmed/ -name \"*.fastq.gz\" -exec basename {} \\; | \\\n  sed 's/\\.fastq.gz$//' | parallel 'STAR --runThreadN 6 ' \\\n    '--genomeDir /local-fs/staff/marcel/GSE149995/reference/star ' \\\n    '--readFilesIn /local-fs/staff/marcel/GSE149995/fastq-trimmed/{}.fastq ' \\\n    '--outSAMtype BAM SortedByCoordinate ' \\\n    '--quantMode GeneCounts ' \\\n    '--outFileNamePrefix /local-fs/staff/marcel/GSE149995/mapping/star/{}_star_'~/miniconda3/bin/multiqc --filename \\\n  ~/Development/2.1.2-Transcriptomics/STAR-mapping.html \\\n  /local-fs/staff/marcel/GSE149995/mapping/star/"},{"path":"from-reads-to-counts.html","id":"combining-star-count-files","chapter":" From Reads to Counts","heading":"2.5.3 Combining STAR Count Files","text":"STAR generates [sample]_star_ReadsPerGene..tab file sample, containing gene IDs count numbers. can merge files single data frame processing (see also Appendix A1 manual). example takes 2nd column count values. Make sure know column use based description . reason ’re merging data might always guaranteed get exact transcripts file.Note: use featurecounts tool, already get single file containing count values samples.Now read count files merge single data frame.now single data frame count values samples transcripts. names however nondescriptive since just SRR identifiers. Since already annotation, ’re going rename columns something meaningful:","code":"\nfile.names <- list.files('/local-fs/staff/marcel/GSE149995/mapping/star/',\n                         pattern = '*_star_ReadsPerGene.out.tab')\n\n## Function for reading in files\nread_sample <- function(file.name) {\n  ## Extract the sample name for naming the column (retaining the 'SRR....' part)\n  sample.name <- strsplit(file.name, \".\", fixed = TRUE)[[1]][1]\n  sample <- read.table(file.name, header = FALSE, sep=\"\\t\", \n                       row.names = NULL, skip = 4)\n  ## Rename the count column\n  names(sample)[2] <- sample.name\n  ## Return a subset containing the transcript ID and sample name columns\n  return(sample[c(1, 2)])\n}\nsetwd('/local-fs/staff/marcel/GSE149995/mapping/star/')\n\n## Read the FIRST sample\ncounts <- read_sample(file.names[1])\n\n## Read the remaining files and merge the contents\nfor (file.name in file.names[2:length(file.names)]) {\n  sample <- read_sample(file.name)\n  counts <- merge(counts, sample, by = 1)\n}\n\n# Set the row names to the transcript IDs\nrownames(counts) <- counts$V1\ncounts <- counts[-1]\n# For each SRR id, find the corresponding column in counts and rename\n# using the name in the 'Condition' column\nfor (i in 1:nrow(setup)) {\n  idx <- grep(setup$Run[i], names(counts))\n  names(counts)[idx] <- setup$Condition[i]\n}\n\n# Show first 6 columns and rows\npander(counts[1:6, 1:6])"},{"path":"EDA.html","id":"EDA","chapter":" Exploratory Data Analysis","heading":" Exploratory Data Analysis","text":"start performing exploratory data analysis steps goal getting grips chosen data set properly identify strategy actual analysis steps. exploration also keep eye quality data. Even though downloadable data ‘processed’, might samples present deviate much group samples (called outlier). Creating basic visualizations data give necessary insight continue.","code":""},{"path":"EDA.html","id":"loading-data","chapter":" Exploratory Data Analysis","heading":"3.1 Loading data into R","text":"one column based text files can read R simply using read.table() function.Follow following steps read data start exploratory data analysis. resulting document treated lab journal log process loading data final analysis steps.Open RStudioCreate new R Markdown document\nGive proper title select PDF format\nGive proper title select PDF formatGive document structure; e.g. create segment (using single hash #) called Exploratory Data Analysis.Whenever add code document make sure readable (keep maximum line length < 100 possible) sufficient documentation either text around code ‘chunks’ (preferred) using comments within code chunk.Read data file(s)\nPreferably use read.table function carefully set arguments. Open file text editor first check contents; header? can set row.names? columns needed? etc.\nNote: data set consists separate files (.e. one per sample) general tips reading data, see Appendix : Batch Loading Expression Data chapter.\nremainder document, try show either contents, structure - case - dimensions relevant R objects\nShow first five lines loaded data set. Including tables markdown document can done using pander function pander(Daróczi Tsegelskyi 2017) R-library.\nGive dimensions (dim() structure (str()) loaded data set.\nCheck output str function see columns expected R data type (e.g. values, factors, character, etc.)\n\nPreferably use read.table function carefully set arguments. Open file text editor first check contents; header? can set row.names? columns needed? etc.Note: data set consists separate files (.e. one per sample) general tips reading data, see Appendix : Batch Loading Expression Data chapter.remainder document, try show either contents, structure - case - dimensions relevant R objects\nShow first five lines loaded data set. Including tables markdown document can done using pander function pander(Daróczi Tsegelskyi 2017) R-library.\nGive dimensions (dim() structure (str()) loaded data set.\nCheck output str function see columns expected R data type (e.g. values, factors, character, etc.)\nShow first five lines loaded data set. Including tables markdown document can done using pander function pander(Daróczi Tsegelskyi 2017) R-library.Give dimensions (dim() structure (str()) loaded data set.Check output str function see columns expected R data type (e.g. values, factors, character, etc.)Examine samples included experiment create many R character objects needed store classification.\nNote: important item often overlooked. See last 4 lines code chunk section Example Data .\ninstance, eight samples divided case/ control columns create variable called case store column indices (numbers!) respective columns data set variable called control remaining four data column indices. later use.\nvariables allows us repeatedly access data rest course, .e. code boxplot(counts[case]) immediately shows plotting case samples counts data clear reading boxplot(counts[c(5:8)]) gives reference indices mean can lead typing mistakes.\nNote: important item often overlooked. See last 4 lines code chunk section Example Data .instance, eight samples divided case/ control columns create variable called case store column indices (numbers!) respective columns data set variable called control remaining four data column indices. later use.variables allows us repeatedly access data rest course, .e. code boxplot(counts[case]) immediately shows plotting case samples counts data clear reading boxplot(counts[c(5:8)]) gives reference indices mean can lead typing mistakes.\ndata sets order samples listed \nloaded data set different order shown GEO\nwebsite. often, names different lacking \ndescription non descriptive IDs. creating\nvariables control need point \ncontrol samples, need make sure right columns\ndata frame. Go Appendix A2;\nannotation order unclear just large number\nsamples might supporting data available can help\nmake sense sample layout.\nwant include external images log better control properties height width individual images can use following code. Note: code example need replace img_to_include.png file path actual image.):See Knitr chunk options page full listing chunk options.","code":"\n# Add the following to the code chunk header to control figure height, width and add a caption:\n#    {r, fig.height=5, fig.width=8, echo=FALSE, fig.cap=\"Image description\"}\n# The `dpi` argument is to remove any scaling\nknitr::include_graphics(\"figures/img_to_include.png\", dpi = NA)"},{"path":"EDA.html","id":"ch3_example_data","chapter":" Exploratory Data Analysis","heading":"3.2 Example Data","text":"section lists (publically available) data set(s) used chapter. chapter contains section new data sets used . Note examples, data different examples one challenges course translating examples data. Keep mind simple copy-pasting code fail reason. examples therefore print input data comparison data.section Normalization onwards, experiment identifier GSE101942 used visualizing normalized raw count values (earlier sections explore unnormalized data). experiment titled “Transcriptome analysis genetically matched human induced pluripotent stem cells disomic trisomic chromosome 21” experimental setup described follows:“12 total polyA selected samples. 6 IPSC samples 3 biological repeats trisomic samples 3 biological repeats disomic samples. 6 IPSC derived neuronal samples 3 biological repeats trisomic samples 3 biological repeats disomic samples.”following code shows load two available files containing raw-count data (one file per 6 samples), stored locally data/gse101942/ folder. simple ‘cleanup’ steps performed aren’t required helpful demonstration purposes (.e., renaming samples identify groups).GSE101942_IPSC_rawCounts.txt: raw count data induced pluripotent stem cells (iPSCs), trisomic disomicGSE101942_Neuron_rawCounts.txt: raw count data IPSC-derived Neurons, trisomic disomic. cells treated remove chromosome 21 iPSCs, described treatment protocol : “Targeted removal CHR 21 IPSC using TKNEO transgene”.See Table 1 Normalization section resulting contents counts dataframe.","code":"\n# Load the two sets of 6 samples\nipsc.data <- read.table('./data/gse101942/GSE101942_IPSC_rawCounts.txt',\n                        header = TRUE)\nneuron.data <- read.table('./data/gse101942/GSE101942_Neuron_rawCounts.txt',\n                          header=TRUE)\n# Merge by rowname (by = 0, see the help of 'merge')\ncounts <- merge(ipsc.data, neuron.data, by = 0, all.x = TRUE, all.y = TRUE,\n                sort = FALSE)\n# Change all NA's introduced by the merge to zeros\ncounts[is.na(counts)] <- 0\n# Show column names coming from the input files\nprint(names(counts))##  [1] \"Row.names\" \"c244A\"     \"c244B\"     \"c243\"      \"c2A\"       \"c2B\"      \n##  [7] \"c2C\"       \"C3_1\"      \"C3_2\"      \"C3_4\"      \"C2_2_1\"    \"C2_2_2\"   \n## [13] \"C2_2_3\"\n# Set the row names to the gene IDs stored in the 'Row.names' column\nrow.names(counts) <- counts$Row.names\n# Remove the gene ID column\ncounts <- counts[-1]\n# Rename samples to include their group name\nnames(counts) <- c(paste0('di_IPSC_r', 1:3),  # Disomic IPSC, replicates 1-3\n                   paste0('tri_IPSC_r', 1:3), # Trisomic IPSC, replicates 1-3\n                   paste0('di_NEUR_r', 1:3),  # Disomic Neuron, replicates 1-3\n                   paste0('tri_NEUR_r', 1:3)) # Trisomic Neuron, replicates 1-3\n# Show results of renaming the samples\nprint(names(counts))##  [1] \"di_IPSC_r1\"  \"di_IPSC_r2\"  \"di_IPSC_r3\"  \"tri_IPSC_r1\" \"tri_IPSC_r2\"\n##  [6] \"tri_IPSC_r3\" \"di_NEUR_r1\"  \"di_NEUR_r2\"  \"di_NEUR_r3\"  \"tri_NEUR_r1\"\n## [11] \"tri_NEUR_r2\" \"tri_NEUR_r3\"\n# Assign column -indices- to variables for later use (i.e. counts[di_NEUR] selects the relevant columns)\ndi_IPSC <- 1:3\ntri_IPSC <- 4:6\ndi_NEUR <- 7:9\ntri_NEUR <- 10:12"},{"path":"EDA.html","id":"EDA_part1","chapter":" Exploratory Data Analysis","heading":"3.3 Visualizing using boxplot and density plot","text":"segment describes basic steps visualization can performed Exploratory Data Analysis. start inspecting unnormalized count data. next steps normalization applied shown Normalization section . mentioned , focus EDA get overview/ perform bit Quality Control data set often requires visualizing data, figures need pretty. Simple figures perfectly fine stage. Try create figures data (keep log) add small description figure pointing anything different expect.Also note steps just selection. Furthermore, make sure every visualization make, add proper axis-labels containing measurement units (important!). might able see, values log-transformed using log2 function often numerical values high range ‘hide’ details plots. See section Fold Change value chapter 4 details. fine use non-log-transformed (simply raw-)data, otherwise use instance boxplot(log2(dataset)) plotting.\nInstead using basic R-plotting library\n(.e. plot, boxplot, hist, etc.)\ncan also opt using (challenging) ggplot2\nlibrary also used boxplot following section.\nconstructing ggplot2 plot feels like learning yet\nanother language, many resources available online can\nfollow.\n","code":""},{"path":"EDA.html","id":"eda_statistics","chapter":" Exploratory Data Analysis","heading":"3.3.1 Statistics","text":"Even basic statistics can give insight data performing 6-number-statistic data columns using summary() function. Note: can also use pander function pretty-print summary knitted output. notice look numbers produced executing function complete data set?","code":""},{"path":"EDA.html","id":"eda_boxplots","chapter":" Exploratory Data Analysis","heading":"3.3.2 Boxplots","text":"visual representation values can shown using boxplot. Boxplots easy create R data frame object just passing data columns. following boxplot shows data experiment separate boxplot sample. allows quick overview spotting irregularities (.e. checking replicates within sample-group show similar expression profiles). course, consider amount data single plot, can hint problems, need look much detail form quality control. Creating boxplot dataframe easy, saw using summary function; data large range maximum average values far apart. create lot outliers plot interesting later , boxplot can either:hide completely using outline = FALSE argument boxplot() (say figure description!), orperform log2() transformation can see .\ncreating plots expression data without performing\nlog-transformations see (instance boxplot) \nrange data large (0 1.000.000) lot\noutliers upper range. makes informative\nfigures, placing simple log2() function call around\ndata almost resolves issue.\n\nHowever, might introduce another problem (taking boxplot \nexample) log-transformation results data \ncontains negative infinity (-Inf) values. \ncaused 0-values data since\nlog2(0) == -Inf. circumvent issue, can add \npseudo count data simply adding value\n1 count values since log2(1) == 0;\nlog2(counts + 1).\n\nsituations boxplot adding \npseudo count value complete dataset can useful. Always add \npseudocount -place, meaning within plotting code instead\noverwriting original data set pseudo count added \nsteps require . cases (see density plot ) can useful \nclearly separate 0-values rest case \n0 < pseudo count < 1 value, 0.1\ncan used generates negative value (-3.3) \nfar away rest data. count values always\n>0 result positive log2 value.\n\nFigure 3.1: Boxplot comparing basic statistics genes across multiple samples\n","code":""},{"path":"EDA.html","id":"eda_density_plot","chapter":" Exploratory Data Analysis","heading":"3.3.3 Density Plots","text":"Another form visualizing data using density plot. method shows distribution (log2-transformed) count data samples allows easy spotting problems. plot commonly used analysing microarrays, still useful comparing complete dataset. code figure show example distribution 12 samples. things note figure huge peak exactly -3.321928 can ignored value calculated log2(0.1) explained box . peak consists 0-values (inactive genes) isn’t important us now. Therefore, added vertical line indicate left-part little interest. Note: lot code extra, simple inspection using line plotDensity(log2(counts + 0.1)) enough, rest extra example code.assess quality data look peaks plot. heigth peak doesn’t matter much, shifted peaks however indicate either lower higher amount reads sequenced sample. can confirmed first plot demonstrated next section library-sizes visualized.\nFigure 3.2: Density plot comparing count distribution 12 samples\n","code":"\n## The affy library has a density plotting function\nlibrary(affy)\n\n## Create a list of 4 colors to use which are the same used throughout this chapter \nlibrary(scales)\nmyColors <- hue_pal()(4)\n\n## Plot the log2-transformed data with a 0.1 pseudocount\nplotDensity(log2(counts + 0.1), col=rep(myColors, each=3),\n            lty=c(1:ncol(counts)), xlab='Log2(count)',\n            main='Expression Distribution')\n\n## Add a legend and vertical line\nlegend('topright', names(counts), lty=c(1:ncol(counts)),\n       col=rep(myColors, each=3))\nabline(v=-1.5, lwd=1, col='red', lty=2)"},{"path":"EDA.html","id":"EDA_part2","chapter":" Exploratory Data Analysis","heading":"3.4 Visualizing using heatmap and MDS","text":"section adds Exploratory Data Analysis techniques measure look distances samples based normalized data. Measuring distances two data objects (samples case) common task cluster analysis compare similarity (low distance indicates similar data). case calculate distances samples visualize heatmap using multidimensional scaling (MDS) technique.previous section used raw count data. might one samples different values (.e. shifted) compared samples. need raw count data use R packages edgeR (Chen et al. 2018) DESeq2 (Love, Anders, Huber 2017), calculating sample distances (used visualizations section) done form normalized data. data can either RPKM/FPKM/TPM/CPM vst-transformed (raw-)read counts. proper method transforming raw read count data using vst method DESeq2 R Bioconductor library shown . ‘variance stabilizing transformation’ normalized data used chapter, chapter 4 normalize using different technique.following code examples shows use library normalize count data GSE101942 experiment vst-normalized data calculate distance metric.new(“standardGeneric”, .Data = function (object, …)standardGeneric(“counts”), generic = structure(“counts”, package = “BiocGenerics”),package = \"BiocGenerics\", group = list(), valueClass = character(0),signature = \"object\", default = NULL, skeleton = (function (object,    ...)stop(gettextf(\"invalid call method dispatch '%s' (default method)\",look Table 1 immediately see huge difference groups genes. gene (A1BG) shows >100 expression tri_NEUR group compared di_IPSC group, within groups high variation indicating might actual expression different. One simple quick method inspection looking total number mapped reads per sample (sum sample/ column), sequencing-depth might vary across samples.Hint:Simple version plot can achieved barplot(colSums(counts) / 1e6)Table 2; Mapped reads per sample (millions) (continued )\nFigure 3.3: Library sizes samples GSE101942\nnumbers shown Table 2 barplot immediately clarify just experimental condition might caused large difference sample counts sequencing depth shows substantial difference . minimum maximum sequencing depth within trisomic neuron group (purple barplot). average expression A1BG gene within group 820, minimum 461 maximum 1010 lot variation. one reason telling us normalize data. following code chunks show DESeq2’s vst method (variance stabilizing transformation). Check Packages tab Rstudio see DESeq2 package installed load library command.","code":"package = \"BiocGenerics\", group = list(), valueClass = character(0),signature = \"object\", default = NULL, skeleton = (function (object,    ...)stop(gettextf(\"invalid call in method dispatch to '%s' (no default method)\","},{"path":"EDA.html","id":"normalization","chapter":" Exploratory Data Analysis","heading":"3.4.1 Normalization","text":"use vst function DESeq2 library must construct DESeqDataSet-object consisting count data combined sample annotation. Since want use (now) performing vst-transformation use basic form simple design sample names annotation (colData argument):now proper DESeqDataSet object can see , containing 56640 rows 12 columns (genes samples) gene symbols row names. Usually object hold data, requirement perform vst transformation good enough now. Next step performing transformation (can take depending size experiment results Large DESeqTransform object) retrieving actual data assay function object contains lot meta-data. table shows updated values now comparable across genes whereas raw count data harder compare.Table 2; VST transformed disomic IPSC count dataNote looking A1BG gene (one displayed figure 3.8), normalized expression values show much less variation across samples (shown, also within sample group). Therefore, assume large difference expression observed earlier might non-existent (might still significant difference though!).","code":"\n# Load the library\nlibrary('DESeq2')\n\n# DESeq2 will construct a SummarizedExperiment object and combine this \n# into a 'DESeqDataSet' object. The 'design' argument usually indicates the \n# experimental design using the condition(s) names as a 'factor', for now we use just '~ 1'\n(ddsMat <- DESeqDataSetFromMatrix(countData = counts,\n                                  colData = data.frame(samples = names(counts)),\n                                  design = ~ 1))## class: DESeqDataSet \n## dim: 56640 12 \n## metadata(1): version\n## assays(1): counts\n## rownames(56640): 5S_rRNA 7SK ... C1orf220 C2orf15\n## rowData names(0):\n## colnames(12): di_IPSC_r1 di_IPSC_r2 ... tri_NEUR_r2 tri_NEUR_r3\n## colData names(1): samples\n# Perform normalization\nrld.dds <- vst(ddsMat)\n# 'Extract' normalized values\nrld <- assay(rld.dds)"},{"path":"EDA.html","id":"eda_distance","chapter":" Exploratory Data Analysis","heading":"3.4.2 Distance Calculation","text":"now normalized data can use distance calculation. standard procedure many data analysis tasks calculates distance metric combination samples use check variation within sample groups. first need transpose matrix (rld) normalized values using t-function, dist function expects different samples rows genes columns. Note output matrix symmetric. Table 3 shows calculated distances within disomic IPSC group distance samples varies 100 140. maximum distance samples 380 demonstrated visualizations next sections.Table 3; Sample distances disomic IPSCs (Euclidean method)","code":"\n# Calculate basic distance metric (using euclidean distance, see '?dist')\nsampledists <- dist( t( rld ))"},{"path":"EDA.html","id":"eda_heatmap","chapter":" Exploratory Data Analysis","heading":"3.4.3 Sample Distances using a Heatmap","text":"(raw-)count data normalized format (TPM, RPKM, etc.), can optionally follow procedure count data create heatmap formats see makes difference. reason RNA-Seq method exists 10 years, still ongoing discussions subject data processing, especially regarding subjects like normalization technque use data analysis.following code block creates heatmap using pheatmap library offers many available heatmap functions. Using annotation data frame (can inspect contents ) identify samples based cell type ploidy. clustering shown heatmap clearly separates data based cell type differences ploidy seems minimal. Looking still, differences within single group minimal meaning - yet - inclined remove outlier-samples. Since 3 samples per category, also lose statistical power eventually remove one samples (also, always check article see remove samples prior data analysis).\nFigure 3.4: Heatmap showing Euclidean distances samples\nresulting heatmap shows interesting comparison across samples order samples purely determined using distance often different order data set. Ideally, want sample groups clustered together see data. particular heatmap clearly shows large difference overal expression induced pluripotent stem cells (iPSCs) differentiated cortical neurons difference - within cell type - disomic trisomic much lower.","code":"\n# We use the 'pheatmap' library (install with install.packages('pheatmap'))\nlibrary(pheatmap)\n\n# Convert the 'dist' object into a matrix for creating a heatmap\nsampleDistMatrix <- as.matrix(sampledists)\n\n# The annotation is an extra layer that will be plotted above the heatmap columns\n# indicating the cell type\nannotation <- data.frame(Cell = factor(rep(1:2, each = 6), \n                                          labels = c(\"IPSC\", \"Neuron\")),\n                         Ploidy = factor(rep(rep(1:2, each = 3), 2),\n                                         labels = c(\"disomic\", \"trisomic\")))\n                         \n# Set the rownames of the annotation dataframe to the sample names (required)\nrownames(annotation) <- names(counts)\n\npheatmap(sampleDistMatrix, show_colnames = FALSE,\n         annotation_col = annotation,\n         clustering_distance_rows = sampledists,\n         clustering_distance_cols = sampledists,\n         main = \"Euclidean Sample Distances\")"},{"path":"EDA.html","id":"eda_mds","chapter":" Exploratory Data Analysis","heading":"3.4.4 Multi-Dimensional Scaling","text":"following code example shows perform Multi-Dimensional Scaling (MDS) displays similarly calculated distances 2D-plot. experiment like two groups samples, hope see two clearly separated clusters formed, however ’ve seen heatmap, sample KO1B showed large deviation also see (confirm) using MDS.\nfigure plotted using ggplot2, \nadvanced method plotting R. plots preferred \nbase-R plotting, always sufficient use just can \nchallenging alter example code shown section. \ndata object plotted shown case contains simple X- \nY-coordinates.\nMDS use slightly different distance metric (Poisson instead default euclidean used dist function). following code shows calculate fitting distance sequencing data called Poisson Distance number optimizations available PoiClaClu library. library specifically designed handle read count data less influenced read-count differences across samples. Note , code usable count data function requires:n--p data matrix observations rows, p features columns. (,j) element x number reads observation mapped feature (e.g. gene exon) j.Table 5: MDS coordinates used plotting distances using Poisson Distance\nFigure 3.5: Multi Dimensional Scaling - Poisson Distance samples\nclearly see separation cell type less clear separation ploidy. spread within disomic groups much lower compared trisomic groups. Note y-axis range 4th x-axis range. sample di_NEUR_r3 shows large difference MDS plot. set genes identified differentially expressed repeat step expectation clear clustering.","code":"\nlibrary('PoiClaClu')\n# Note: uses the raw-count data, PoissonDistance performs normalization\n# set by the 'type' parameter (uses DESeq)\ndds <- assay(ddsMat)\npoisd <- PoissonDistance( t(dds), type = \"deseq\")\n# Extract the matrix with distances\nsamplePoisDistMatrix <- as.matrix(poisd$dd)\n# Calculate the MDS and get the X- and Y-coordinates\nmdsPoisData <- data.frame( cmdscale(samplePoisDistMatrix) )\n\n# And set some better readable names for the columns\nnames(mdsPoisData) <- c('x_coord', 'y_coord')\n# Separate the annotation factor (as the variable name is used as label)\ngroups <- factor(rep(1:4, each=3), \n                 labels = c(\"di_IPSC\", \"tri_IPSC\", \"di_NEUR\", \"tri_NEUR\"))\ncoldata <- names(counts)\n\n# Create the plot using ggplot\nggplot(mdsPoisData, aes(x_coord, y_coord, color = groups, label = coldata)) + \n  geom_text(size = 4) +\n  ggtitle('Multi Dimensional Scaling') +\n  labs(x = \"Poisson Distance\", y = \"Poisson Distance\") +\n  theme_bw()\n## Using base-R plot and legend functions:\n# plot(mdsPoisData, col=rep(myColors, each=3), pch=20, lwd=2)\n# legend(x=-20000, y=11000, legend = levels(groups), col=myColors, pch=20)"},{"path":"EDA.html","id":"cleaning-data","chapter":" Exploratory Data Analysis","heading":"3.5 Cleaning Data","text":"Conclude chapter cleaning dataset visualizations show necessity . Cleaning case means removing complete samples can classified outlier within group. tasks course, clear advice decide remove one samples. clear rule group must retain least three samples. one three samples visibly deviates replicates, even normalization, stay data set mentioned analysis report/ final article noise may introduced sample.\nNote: satisfied data reported findings, return steps next chapters. mentions creating heatmap instance, means heatmap expression data instead sample distances ’ve shown chapter.","code":""},{"path":"chapter-4.html","id":"chapter-4","chapter":" Discovering Differentialy Expressed Genes (DEGs)","heading":" Discovering Differentialy Expressed Genes (DEGs)","text":"first important ‘real’ analysis step finding genes show difference expression sample groups; differentially expressed genes (DEGs). concept might sound rather simple; calculate ratios genes samples determine fold-change (FC) denoting factor change expression groups. , filter genes actually show difference. give list genes showing different behaviour across samples, need focus genes show difference, observed difference also statistically significant!determine whether gene can classified significant DEG going use existing libraries, one article might refer . ’re going brief manual examination count data make sure properly understand ’re looking : significance Fold-Change value.","code":""},{"path":"chapter-4.html","id":"deg-preprocessing","chapter":" Discovering Differentialy Expressed Genes (DEGs)","heading":"4.1 Pre-processing","text":"Given results exploratory data analysis performed chapter 3, might concluded one samples show () deviating expression patterns compared samples group. mentioned , enough (> 3) samples group, might opt remove sample reduce noise statistical tests sensitive . Since performing analysis steps programatically also easy test DEGs without sample(s) question see removal results lower p-values (= higher significance). Remember always properly document choice reasoning!perform manual steps data using libraries, re-normalization different technique. simple method normalizing count data calculate fragments per million mapped fragments (FPM, equal CPM) value transform log2. Opposed FPKM RPKM include gene-length calculation (likely don’t ). dataset also includes FPKM RPKM (pre-normalized), allowed use data .FPM CPM simply dividing count values sample total number reads sample divided 1 million shown :Another step can - might guided article - filter (partially) inactive genes. required using edgeR library (DESeq2 includes low-count filtering methods), good exercise manually use calculating lFC values. data sets available contains lot 0-measurements (see instance density plot previous chapter); genes reads mapped. R-studio Environment tab, click data set (perform View(data) command click one sample columns order data ascending. now likely (unless bacterial data instance) see lot zero values columns.experiment two groups, three replicates , three six samples 0-reads mapped, often advised remove gene completely. can subjective experiment, might expected (thus important) comparing different tissues knockout experiments. Also, genes () low read count (< 5) can give high (artificial) FC value (see left hand side MA-plot). Comparing two samples one value 2 11, reads -regulated gene factor 5.5 might actually just noise!Assignment: Search article advice filter zero values low count genes. nothing stated subject, think tactic (search literature/ online!). perfectly fine discuss peers.\ndescribe aspect; filter data, clearly explain (likely article stated proper reason)\nPerform filtering data set. likely need use one apply, rowSums, rowMeans functions, combined maybe , functions.\nManually verify rows removed correctly filtered.\nProperly document many genes filtered !\ndescribe aspect; filter data, clearly explain (likely article stated proper reason)Perform filtering data set. likely need use one apply, rowSums, rowMeans functions, combined maybe , functions.Manually verify rows removed correctly filtered.Properly document many genes filtered !","code":"\n# Perform a naive FPM normalization\n# Note: log transformation includes a pseudocount of 1\ncounts.fpm <- log2( (counts / (colSums(counts) / 1e6)) + 1 )"},{"path":"chapter-4.html","id":"deg-fc","chapter":" Discovering Differentialy Expressed Genes (DEGs)","heading":"4.2 The Fold Change Value","text":"FC usually given calculated log2 case/control ratio. example, gene average expression 30 mapped reads control group 88 reads experiment group, ratio case/control 2.93. Ratio values > 1 indicate increased expression experiment relation control values 0 1 indicate lower expression. log2 transformed value ratio calculated log2(2.93) results 1.55. counts reversed, ratio 30/88, 0.34. previously calculated value 2.93 means 3-fold regulation 0.34 value means 3-fold regulation can see range numbers different. Comparing log2 values 1.55 -1.55 compares much better.easy calculate FC genes , simple FC value doesn’t mean much, yet! still need use power lie within replicates sample group. Using replicates, want determine observed FC just biological noise sequencing error.Assignment: Create histogram log2-FC values genes. Apply following steps FPM-normalized data least one comparison (need possible combinations):\naverage replicates two groups separately add new columns FPM-normalized data set,\ncalculate log Fold Change (lFC); simply subtract one group (averaged-column) another, usually experiment-control\nNote since already log-transformed data, calculate ratio (experiment/control). Also, Logarithm quotient rule, divide average columns ratio calculation, subtract: \\(log_b{(x/y)} = log_b{(x)} - log_b{(y)}\\)\n\nplot data using hist(logFC_column, breaks=60) function (change breaks argument needed),\nAdd two vertical lines -1 1 (using abline(v=...)) indicate significance (2-fold change).\ndata shows low fold-change values might see genes outside boundary. mean find significantly changed genes (DEGs).\n\naverage replicates two groups separately add new columns FPM-normalized data set,calculate log Fold Change (lFC); simply subtract one group (averaged-column) another, usually experiment-control\nNote since already log-transformed data, calculate ratio (experiment/control). Also, Logarithm quotient rule, divide average columns ratio calculation, subtract: \\(log_b{(x/y)} = log_b{(x)} - log_b{(y)}\\)\nNote since already log-transformed data, calculate ratio (experiment/control). Also, Logarithm quotient rule, divide average columns ratio calculation, subtract: \\(log_b{(x/y)} = log_b{(x)} - log_b{(y)}\\)plot data using hist(logFC_column, breaks=60) function (change breaks argument needed),Add two vertical lines -1 1 (using abline(v=...)) indicate significance (2-fold change).\ndata shows low fold-change values might see genes outside boundary. mean find significantly changed genes (DEGs).\ndata shows low fold-change values might see genes outside boundary. mean find significantly changed genes (DEGs).calculated FC shows large change expression groups means nothing variation within group high. reason use form statistical test checks variation group difference groups samples. simplest form usually comes using t-test:“can used determine two sets data significantly different ”often instead manually performing t-tests ANOVA analysis, article mentions use one R libraries software packages finding DEGs. exactly also chapter. output finding DEG’s always includes - limited - list p-values; one p-value per gene. value indicates whether gene statistically-significant differentially expressed gene (SDEG) find genes interest need get genes p-value threshold (.e. 0.05).","code":""},{"path":"chapter-4.html","id":"deg-bioconductor","chapter":" Discovering Differentialy Expressed Genes (DEGs)","heading":"4.3 Using Bioconductor Packages","text":"section demonstrates use two packages perform DEG-analysis count data. many packages available Bioconductor RNA-Seq analysis, DSS, EBSeq, NOISeq BaySeq, focus edgeR DESeq2 processing count-based data. Chances one two packages mentioned article described use R statistical analysis widely accepted methods processing gene expression data.packages apply normalization methods (described sections ) therefore work raw count data. data consists expected count data, need round values (using floor function). can choose one two packages use, use since increase statistical power hard use understand model data ’ll show next.Assignment: Perform data-analysis using least one two packages. however article mentioned different R package analysis (instance limma package) allowed use one instead. non-R packages can done extra (dataset can always analyzed using edgeR DESeq2). Optionally, use one package can compare results packages.\nStart documenting relevant information article, also used software analysis.\nStart documenting relevant information article, also used software analysis.two listed packages come complete extensive manual. Therefore, course manual (document) guide exactly use ; find . manuals can found PDF BioConductor website, opened within R using vignette function. see available manuals (note; different function-help documentation) execute function without arguments: vignette(). lists vignettes edgeR DESeq2 :vignette(\"DESeq2\") (whole course partly designed document!)vignette(\"edgeR\")\n","code":""},{"path":"chapter-4.html","id":"deg-design","chapter":" Discovering Differentialy Expressed Genes (DEGs)","heading":"4.3.1 The Design (matrix)","text":"packages need properly specify samples grouped. seen examples using R factor object heatmap MDS visualizations tell groups samples group sample belongs. Reading documentation packages shows important part performing DE-analysis. example, following code shown edgeR documentation page 8 two sample groups defined (numbered 1 2), placed factor object used input model.matrix function.can often complex part seems theoretical, practice ’s difficult. really grasp topic strongly adviced watch two videos define design:StatQuest YouTube: Using Linear Models t-tests ANOVA, Clearly ExplainedStatQuest YouTube: Design Matrices Linear Models, Clearly ExplainedThese videos expects basic understanding topics already discussed statistics course, otherwise mentioned topic videos videos available explaining . See video-index handy overview topics discussed channel.mentioned need think question want ask; difference want know? two sample groups shown question rather easy; ‘genes show effect case/ control samples?’. two sample groups however question becomes complex. experiment comparing influence three kinds drugs (thus three groups) combined effect time, want focus influence drugs time? valid questions define way create design matrix.Documentation subject plenty, however often contains overwhelming information. page contains valuable details (can safely start reading Choice design), including following text based example used edgeR documentation:“examples cover , use linear models make comparisons different groups. Hence, design matrices ultimately work least two columns: intercept column, consists column 1’s, second column, specifies samples second group. case, two coefficients fit linear model: intercept, represents population average first group, second coefficient, represents difference population averages second group first group. latter typically coefficient interested performing statistical tests: want know difference two groups.”two sample groups want change question (.e. test influence different group), read section Releveling factor.","code":"\ngroup <- factor(c(1,1,2,2), labels = c(\"Control\", \"Case\"))\n(design <- model.matrix( ~ group))##   (Intercept) groupCase\n## 1           1         0\n## 2           1         0\n## 3           1         1\n## 4           1         1\n## attr(,\"assign\")\n## [1] 0 1\n## attr(,\"contrasts\")\n## attr(,\"contrasts\")$group\n## [1] \"contr.treatment\""},{"path":"chapter-4.html","id":"deg-deseq2","chapter":" Discovering Differentialy Expressed Genes (DEGs)","heading":"4.3.2 DESeq2","text":"used DESeq2 library DEG analysis re-use DESeqDataSet object adviced create new object proper design formula instead ~ 1 used .need normalize data using previously used vst function DESeq2 library normalize count data follows:“DESeq computes scaling factor given sample computing median ratio, gene, read count geometric mean across samples. uses assumption genes DE uses median ratios obtain scaling factor associated sample.”example data four sample groups re-use annotation data frame created EDA: Heatmap section (printed ). allows us comparisons based celltype ploidy samples. One minor note importance first level factor taken reference, expression analysis, control group. example (see table ) IPSC control-celltype disomic control-ploidy. simple case/ control study, first level (check levels() function) control group. Use relevel() function levels need switched.proper DESeqDataSet need run DESeq function object. , using results function can extract DEGs DESeqResults object. Applying summary function object(s) shows number regulated DEGs can seen lists impact trisomy IPSC cell type (~800 DEGs) comparing IPSC Neuron disomic (control) cell types (can seen, huge amount genes affected differentiation stem cell cortical neurons):output results function contains following columns gene:call results function depends heavily experiment. can see output summary function, details given comparison shown (also, default p-value 0.1 used instead 0.05).Depending design used create DESeqDataSet , one comparisons can made (applying DESeq function calculates filter results function). Read help results function carefully; especially regarding contrast argument define comparison retrieve. following code can used example experiment get DEGs comparing celltypes using adjusted p-value 0.05:DESeq2 library contains number plotting functions can applied DESeqResults object (output results function), notable plotMA function. Note first MA-plot shown high range log fold changes (-10, 10) maximum value 22.4 (shown triangle stating outside plotting range). log fold change 22 means >5 million increased expression seems artificially high. DESeq2 includes function perform downstream processing estimated log fold change values called lfcShrink advised always run afterwards. reason executing function described vignette :“useful visualize MA-plot shrunken log2 fold changes, remove noise associated log2 fold changes low count genes without requiring arbitrary filtering thresholds.”\nFigure 4.1: MA-plots standard DESeq2 output (left) shrinking ‘lfcShrink’ (right)\nLinksAnalyzing RNA-seq data DESeq2 comprehensive guide analyzing RNA-Seq data using DESeq2 (part document used material !). adviced read first sections guide take good look index document many interesting sections might help later.Publication, accompanying article showing differences performance compared methods packages.","code":"\nout of 44504 with nonzero total read count\nadjusted p-value < 0.05\nLFC > 0 (up)       : 432, 0.97%\nLFC < 0 (down)     : 417, 0.94%\noutliers [1]       : 219, 0.49%\nlow counts [2]     : 17004, 38%\n(mean count < 6)\n[1] see 'cooksCutoff' argument of ?results\n[2] see 'independentFiltering' argument of ?results\n\n\nout of 44504 with nonzero total read count\nadjusted p-value < 0.05\nLFC > 0 (up)       : 8224, 18%\nLFC < 0 (down)     : 9594, 22%\noutliers [1]       : 219, 0.49%\nlow counts [2]     : 11906, 27%\n(mean count < 2)\n[1] see 'cooksCutoff' argument of ?results\n[2] see 'independentFiltering' argument of ?results\n\n\nres <- results(dds, contrast = c(\"Cell\", \"IPSC\", \"Neuron\"),\n               alpha = 0.05)"},{"path":"chapter-4.html","id":"deg-edger","chapter":" Discovering Differentialy Expressed Genes (DEGs)","heading":"4.3.3 edgeR","text":"One mature libraries RNA-Seq data analysis edgeR library available Bioconductor. complete (sometimes bit complex) manual available need read Chapter 2 focus 2.1 2.7, 2.9 - complex design - 2.10. Section 1.4 (Quickstart) shows code example steps needed DEG analysis using count data using two glm methods (quasi-likelihood F-tests likelihood ratio tests). steps shown identical non-glm method calculating fit object can replaced performing exactTest function shown section 2.9.2.edgeR library normalize count data follows:“trimmed means M values (TMM) Robinson Oshlack, implemented edgeR, computes scaling factor two experiments using weighted average subset genes excluding genes exhibit high average read counts genes large differences expression.”Running edgeR requires raw count data together grouping-factor packaged DGEList object (DGEList() function). Furthermore, proper model.matrix object (see section design) needed input estimateDisp function. exact steps take (variations DESeq2) must searched documentation linked .analysis done can retrieve actual results topTags function:Quitting lines 696-762 [.edger-TEST] (gene_expression.Rmd)\nError h():\n! error evaluating argument ‘x’ selecting method function ‘.data.frame’: object ‘qlf’ found\nBacktrace:\n1. pander::pander(.data.frame(topTags(qlf)), caption = “significant genes given edgeR”)\n4. base::.handleSimpleError(<fn>, “object ‘qlf’ found”, base::quote(topTags(qlf)))\n5. base (local) h(simpleError(msg, call))package also contains plotting methods can use intermediate steps analysis. instance, calculating normalization factors (calcNormFactors), can perform multi-dimensional scaling plotMDS function similar ’ve done manually chapter three:\nFigure 4.2: edgeR MDS plot based calculated log2 fold changes\ndispersion running estimateDisp function plotBCV function:\nFigure 4.3: edgeR plot several dispersion methods\nlog-fold changes genes, output exactTest function (output et DGEExact object) plotSmear function. abline shows log-FC threshold:LinksDifferential Expression Analysis using edgeR tutorialAnother tutorial hosted GitHub","code":"\nplotMDS(y)\nplotBCV(y)\ndeGenes <- decideTestsDGE(qlf, p=0.05)\ndeGenes <- rownames(qlf)[as.logical(deGenes)]\nplotSmear(qlf, de.tags=deGenes)\nabline(h=c(-1, 1), col=2)"},{"path":"chapter-5.html","id":"chapter-5","chapter":" Data Analysis and Visualization","heading":" Data Analysis and Visualization","text":"point one sets DEGs experiment. might just one set might large number comparisons made set DEGs per comparison / results multiple approaches (manual, packages). chapter offers guidance visualize, summarize - finally - connect biological relevance results. Browse chapter first choosing visualization(s) make data! article contains relevant figures like try reproduce (given enough data available) can done well instead visualizations shown . Note however Volcano plot described required.","code":""},{"path":"chapter-5.html","id":"volcano","chapter":" Data Analysis and Visualization","heading":"5.1 Volcano Plot","text":"volcano plot often first visualization data statistical tests completed. plot shows data genes highlight genes considered DEG using thresholds (adjusted) p-value fold-change. Many articles describe values used thresholds methods section, otherwise good default 0.05 adjusted p-value around 1.0 log-FC.use EnhancedVolcano library creates ggplot2 volcano plot shown . packaged function deseq.volcano just create two figures need change input data set plot title. resulting scatterplot places -log10(pvalue) values y-axis log-FC x-axis. often results plot representing volcanic ‘eruption’ fold-change influences spread significance height. Coloring done based thresholds (-log10(0.05) adjusted log2(2) log-FC) follows:(dark) gray: significant; FC adjusted p-value thresholds.green: high FC, adjusted p-value threshold.blue: adjusted p-value indicates significance, FC thresholdred: significant DEGYou can adjust pCutoff FCcutoff values depending requirements, see elaborate documentation examples GitHub. code wraps EnhancedVolcano function custom function (deseq.volcano) can used multiple comparisons make.Notes:used edgeR instead need change names columns (.e. log2FoldChange == logFC padj == FDR) used example.example many DEGs annotate name plots. EnhancedVolcano plot function many options annotating significant genes, see linked manual information.","code":"\nlibrary(EnhancedVolcano)\n\n## Simple function for plotting a Volcano plot, returns a ggplot object\ndeseq.volcano <- function(res, datasetName) {\n  return(EnhancedVolcano(res, x = 'log2FoldChange', y = 'padj',\n                         lab=rownames(res),\n                         title = paste(datasetName, \"trisomic vs disomic\"),\n                         subtitle = bquote(italic('FDR <= 0.05 and absolute FC >= 2')),\n                         # Change text and icon sizes\n                         labSize = 3, pointSize = 1.5, axisLabSize=10, titleLabSize=12,\n                         subtitleLabSize=8, captionLabSize=10,\n                         # Disable legend\n                         legendPosition = \"none\",\n                         # Set cutoffs\n                         pCutoff = 0.05, FCcutoff = 2))\n}\n\n## Note: input data is the corrected DESeq2 output using the 'lfcShrink' function (see chapter 4)\ndeseq.volcano(res = res.ipsc.lfc, datasetName = \"IPSC\")\ndeseq.volcano(res = res.neuron.lfc, datasetName = \"Neuron\")\n## Note that the most simplest case is:\nEnhancedVolcano(res.ipsc, x = 'log2FoldChange', y = 'padj', lab = rownames(res.ipsc))"},{"path":"chapter-5.html","id":"venn-diagram","chapter":" Data Analysis and Visualization","heading":"5.2 Venn Diagram","text":"Another common visualization Venn-diagram. data set ’ve shown two comparisons; trisomic vs disomic two cell types. Using Venn-diagram can show shared unique DEGS IPSC trisomic vs IPSC disomic NEUR trisomic vs NEUR disomic comparisons.Creating Venn-diagram can good method selecting (sub)set DEGs, instance selecting set DEGs shared groups (Intersection, see Wikipedia page Set Theory way around, DEGs unique group (Difference). Keep mind continuing clustering section likely useful cluster combinations might .Code examples:Data preparation (first code-chunk): creating vectors DEGs packages. DEGs gene names (taken row names DEG-analysis).\noutput DESeq2 results lfcShrink functions stored res.ipsc.lfc res.ipsc.neuron objects.\noutput DESeq2 results lfcShrink functions stored res.ipsc.lfc res.ipsc.neuron objects.Using VennDiagram library (second code-chunk): good complex/ large number comparison, highly configurable difficult use. see examples, run following function R console: example('draw.quad.venn')Using gplots library (third code-chunk): Easier use support complex Venn diagrams, good adjusting appearance.See blog post “VENN DIAGRAM R RSTUDIO: MILLION WAYS” methods creating Venn diagrams.two created vectors gene names can used either taking lengths (VennDiagram library) just names (gplots library).\nFigure 5.1: Venn diagram comparing IPSC NEUR DEGs using ‘VennDiagram’ library\nmentioned, following easier alternative (especially > 2 groups) venn function gplots library. needs list containing either row-numbers gene-names DEGs easier (offers less adjustability make prettier).\nFigure 5.2: Venn diagram comparing IPSC NEUR DEGs using ‘gplots’ library\nanalysis result similar amount DEGS (849 IPSC 769 Neuron cell types) 96 DEGs shared cell types thus unique.","code":"\n## Data preparation\npval_threshold <- 0.05\nipsc.degs <- row.names(res.ipsc.lfc[which(res.ipsc.lfc$padj <= pval_threshold), ])\nneur.degs <- row.names(res.neuron.lfc[which(res.neuron.lfc$padj <= pval_threshold), ])\n## Venn-diagram using the `VennDiagram` library (see below for alternative method)\nlibrary(VennDiagram)\n# Arguments for a pairwise (two-sets) venn-diagram are sizes for set1, set2 and overlap (intersect)\n# Many more functions are available for triple, quad and quantuple diagrams (starting with 'draw.***')\nvenn.plot <- draw.pairwise.venn(length(ipsc.degs),\n                                length(neur.degs),\n                                # Calculate the intersection of the two sets\n                                length( intersect(ipsc.degs, neur.degs) ),\n                                category = c(\"IPSC Trisomic vs Disomic\", \"NEURON Trisomic vs Disomic\"), scaled = F,\n                                fill = c(\"light blue\", \"pink\"), alpha = rep(0.5, 2),\n                                cat.pos = c(0, 0))\n\n# Actually plot the plot\ngrid.draw(venn.plot)\n## Alternative method using the `gplots` library\nlibrary(gplots)\n# Create a Venn-diagram given just the list of gene-names for both sets\nvenn(list(\"IPSC Trisomic vs Disomic\" = ipsc.degs,\n          \"NEURON Trisomic vs Disomic\" = neur.degs), )"},{"path":"chapter-5.html","id":"clustering","chapter":" Data Analysis and Visualization","heading":"5.3 Clustering","text":", can create heatmap found DEGs samples colors show fold-change value. Usually heatmap shown two color gradient; red (downregulated) green (upregulated). however optional using default colors pheatmap example perfectly fine long proper legend present.multiple methods creating heatmap (one ’ve already used) directly apply clustering visualization. clustering can applied expression pattern gene (row-clustering), expression pattern sample (column-clustering) (default pheatmap).following example commonly found publications shows log2 Fold Change values comparison done. first fifth columns show log2(FC) calculated DESeq2 two trisomic vs disomic comparisons. remaining columns log fold changes separate trisomic replicates vs average disomic expression (manually calculating Log(FC) done section 4.1 using DESeq2 normalized count data rld.dds object). DESeq2 calculated FC values show greater changes compared manually calculated FC values. value showing manually calculated values see replicates show patterns within groups. however optional unless one comparison create heatmap single column.Note column-clustering turned type plot setting cluster_cols = FALSE (replicates group, clustering add useful information).\nFigure 5.3: Heatmap showing fold changes Trisomic vs Disomic comparisons cell types. shown genes shared DEGS (96 shown Venn-diagrams) also filtered absolute log2(FC) value >2.\nlimitation heatmap show 100 DEGs found; just doesn’t fit well single figure causes clustering especially gene names/ labels unreadable. therefore always adviced filter adjusted p-value, also log-FC value reduce number genes shown. Alternatively data properly annotated, instead showing gene names row, showing GO-term clustering might reveal expression patterns certain gene groups.","code":""},{"path":"chapter-5.html","id":"pathway-analysis","chapter":" Data Analysis and Visualization","heading":"5.4 Pathway Analysis","text":"cases ’ll see hundreds even thousands DEGs result analysis. large amounts DEGs much visualization approaches easily say something biological context many genes.One approach tackling large set DEGs pathway analysis different methods genes grouped pathway get overview affected pathways experiment.section demonstrate two methods analysis, one using online platform gene-annotation enrichment analysis R-method signaling pathway impact analysis.","code":""},{"path":"chapter-5.html","id":"kegg-pathway-visualization-10-degs","chapter":" Data Analysis and Visualization","heading":"5.4.1 KEGG-pathway Visualization (>10 DEGs)","text":"Knowing interesting KEGG-pathway(s) upfront (.e. listed article) allows visualizations applied selected pathway. can use pathview website Bioconductor library:“… maps renders wide variety biological data relevant pathway graphs. users need supply data specify target pathway. Pathview automatically downloads pathway graph data, parses data file, maps user data pathway, render pathway graph mapped data.”section demonstrates use bioconductor package visualizing DEGs onto pathways, also change expression gene.Table continues belowSYMBOL, GENENAME, ENSEMBL, ENSEMBLPROT, UNIGENE, UNIPROT, ACCNUM, ENSEMBLTRANS, REFSEQ, ENZYME, TAIR, PROSITE ORF","code":"\n## Bioconductor library (install with bioclite() if missing)\nlibrary(pathview)\n\n## Example pathway IDs (for human, change organism key for other organisms)\ndata(\"paths.hsa\")\npander(head(paths.hsa, n=5))\n## Check 'gene.idtype' argument possibilities\ndata(gene.idtype.list); \npander(gene.idtype.list)\n## Prepare data for visualization\ndeseq.degs.logfc <- subset(deseq.results, padj < pval_threshold, select = log2FoldChange)\n\npander(head(deseq.degs.logfc))\npathview(gene.data=deseq.degs.logfc,\n         pathway.id=\"00020\", # TCA-cycle\n         species=\"hsa\"     # Organism key\n)"},{"path":"datasets.html","id":"datasets","chapter":"A Discovering RNA-Seq Data Sets","heading":"A Discovering RNA-Seq Data Sets","text":"chapter describes method finding public RNA-Seq data sets interest. data set context course refers data belonging certain gene-expression experiment, usually consisting number sequencing-samples combined meta-data describing experiment.","code":""},{"path":"datasets.html","id":"finding-public-data","chapter":"A Discovering RNA-Seq Data Sets","heading":"A.1 Finding Public Data Sets of Interest","text":"number large databases online offer access thousands - published - experiments focus searching downloading gene expression data high-throughput (NGS) sequencing techniques (RNA-Seq).thousands freely available data sets possible start performing research without need actually performing lab-experiment. common conditions, organisms tissues can download samples compare multiple studies find novel relations genes conditions verify experiments performed laboratory. instance, power public data sets demonstrated jointly three alumni article called Calling genotypes public RNA-sequencing data enables identification genetic variants affect gene-expression levels.(al. 2015).’s method section begins sentence “downloaded raw reads available human RNA-seq data sets …” amount data work become clear later . interested expression Quantitative Trait Loci (eQTLs) Allele-Specific Expression (ASE) analysis, please read interesting paper.start diving large complex databases potentially containing thousands experiments millions samples terabytes data, need get rough idea like found something interest order know looking first place. definitive guide protocol can followed processing analyzing gene-expression data set, following goals can considered:Re-(part ) analysis described accompanying publication.\ndata sets - except ones added recently - accompanied publication, can gather lot information regarding experimental design results just paper. often, researchers already performed interesting research data set therefore good exercise try reproduce results. Challenges getting proper understanding techniques used researchers (article) translating analysis steps. Often though, details clear article requires interpretation results performing analysis steps arrive conclusion.\ndata sets - except ones added recently - accompanied publication, can gather lot information regarding experimental design results just paper. often, researchers already performed interesting research data set therefore good exercise try reproduce results. Challenges getting proper understanding techniques used researchers (article) translating analysis steps. Often though, details clear article requires interpretation results performing analysis steps arrive conclusion.Alternatively can come ideas (biological) question(s) try answer given data set instead redoing published work.\ncourse challenging suitable data sets.\nmight notice publications focus small set genes instead whole transcriptome. instance, perform experiment yeast want measure activity genes involved alcohol fermentation, researchers might look genes Glycolysis pathway. leaves 6000+ yeast genes explore possibly come novel relations gene expression experimental setup.\nanalysis approach explained paper focusing pre-selected list genes, depending experiment might come comparisons done original research. instance, experiment maternal age (birth) correlated autism offspring, paternal age also known addressed publication. subject research (spoiler alert; clear conclusion can drawn including extra factor…).\ncourse challenging suitable data sets.might notice publications focus small set genes instead whole transcriptome. instance, perform experiment yeast want measure activity genes involved alcohol fermentation, researchers might look genes Glycolysis pathway. leaves 6000+ yeast genes explore possibly come novel relations gene expression experimental setup.analysis approach explained paper focusing pre-selected list genes, depending experiment might come comparisons done original research. instance, experiment maternal age (birth) correlated autism offspring, paternal age also known addressed publication. subject research (spoiler alert; clear conclusion can drawn including extra factor…).Another type project evaluate different methods either normalization statistical analysis either confirm published results find novel genes involved experiment.\ncan viewed technical research subject ok pursue, however final report mainly concern biological impact findings.\nmeans accompanying publication shows conclusive results, probably come acknowledging results (therefore, extend point 1 list). conclusive result hope find results case might become hard formulate good conclusion end.\nHowever, publication isn’t specific states conclusion ‘Benzene exposure shows increased risk leukemia’ followed list genes might involved, try see can find genes might involved changing analysis approach.\n\ncan viewed technical research subject ok pursue, however final report mainly concern biological impact findings.\nmeans accompanying publication shows conclusive results, probably come acknowledging results (therefore, extend point 1 list). conclusive result hope find results case might become hard formulate good conclusion end.\nHowever, publication isn’t specific states conclusion ‘Benzene exposure shows increased risk leukemia’ followed list genes might involved, try see can find genes might involved changing analysis approach.\nmeans accompanying publication shows conclusive results, probably come acknowledging results (therefore, extend point 1 list). conclusive result hope find results case might become hard formulate good conclusion end.However, publication isn’t specific states conclusion ‘Benzene exposure shows increased risk leukemia’ followed list genes might involved, try see can find genes might involved changing analysis approach.Although last project goal many risks involved need sound project proposal, might result interesting project. shown linked article, possible combine data multiple experiments.\ninstance find two related experiments (.e. researching effect certain drug) measuring expression different tissue types (.e. liver brain tissue). analyzing experiments, present set genes show effect tissues - interestingly - genes showing effect one tissue.\ninstance find two related experiments (.e. researching effect certain drug) measuring expression different tissue types (.e. liver brain tissue). analyzing experiments, present set genes show effect tissues - interestingly - genes showing effect one tissue.list guideline kept mind browsing suitable data sets extended week three. Also, don’t worry terms unclear, getting good understanding used terminology part first weeks","code":""},{"path":"datasets.html","id":"ncbi-geo","chapter":"A Discovering RNA-Seq Data Sets","heading":"A.2 The NCBI Gene Expression Omnibus","text":"easiest method exploring finding interesting data sets simply using web-browser access data-repository. One repository NCBI Gene Expression Omnibus (GEO). repository primarily used store micro array data sets describe experiments, linking raw data, processed data accompanying publication. however many sources data browsable GEO, project limit searches experiments using RNA-Seq technology 60.000 experiments listed. Following “Expression profiling high throughput sequencing” link summary table get complete listing can find interesting experiment. Besides finding experiment interested , important requirements need account browsing listed . Note point need download data yet, described chapter 3.experiment must published publication must fully available, free otherwise Hanze University library.sample group, minimum three replicates must present. one groups less three replicates, group used (rest data might still usable minimum two groups 3 replicates).\nReplicates often easy spot. lucky samples named “WildType - rep 2” instance often need read “Overall design” part page click samples (“GSM…” ID) make sure.\nReplicates often easy spot. lucky samples named “WildType - rep 2” instance often need read “Overall design” part page click samples (“GSM…” ID) make sure.available data must contain least count data.\nData offered multiple formats, always including RAW data (actual reads don’t use), also processed data FPKM, RPKM TPM (see interesting blog-post discussing use data). tools R analysis libraries use downstream analysis rely unnormalized unprocessed data count data. counts simply represent number reads mapped transcript.\noften hardest part assess; sometimes clear provided data actual raw count values. Usually, can find (sometimes rather cryptic) information clicking sample ID ‘Data processing’ section. paragraph text often mention ‘Supplementary_files_format_and_content’ describes actual content downloadable files. Following number examples valid invalid data provided experiment. single-sample pages can find description. Unfortunately, authors experiment free explain data please, likely encounter exact descriptions:\nValid: Global effects SUPT4H1 RNAi gene expression HEK293 cells; description: Supplementary_files_format_and_content: tab-delimited, gzip-compressed text file matrix raw quantitation results sample. Mentions raw-anything likely valid.\nValid: Transcriptome profiling normal transformed mouse urothelial cells; description: Supplementary_files_format_and_content: raw_counts_unfiltered.txt.gz: Raw counts.. Even though 4 different file types available download, one contains actual raw count data.\nValid: Quantitative Analysis PPARD RNA-Seq Transcriptomes Mouse Gastric Corpus Epithenial Cells Next Generation Sequencing (NGS); description: Supplementary_files_format_and_content: raw counts sequencing reads genes. Single file download containing raw counts.\nValid: Impact shRNA-mediated KLF4 -modulation transcriptome profile human keratinocyte precursor cells; description: Supplementary_files_format_and_content: raw normalized counts.\nValid: Role TAF5L TAF6L controlling mouse ESC state; description: Supplementary_files_format_and_content: Comma-separated text files containing raw counts, normalized counts, FPKM expression values TPM expression values.. common include normalized raw values.\nInvalid: Co-chaperone Mzb1 key effector Blimp1 plasma cell differentiation b1 integrin function; description: Supplementary_files_format_and_content: normalized gene expression values (FPKM) calculated Cuffnorm (v2.2.1) provided TXT format. clearly states normalized data (FPKM values) available.\nInvalid: EGFR-targeted therapy-induced resistance mechanism glioblastoma; description: Supplementary_files_format_and_content: .txt file including gene name, RefSeq id, locus, RPKM value gene sample. RPKM values available normalized raw count values, usable.\nInvalid: RNA-seq subcellular fractions nESC EpiLC; description: Supplementary_files_format_and_content: _genes.fpkm_table.txt: Tab-delimited text files include FPKM values replicate.*. FPKM data available.\n\nData offered multiple formats, always including RAW data (actual reads don’t use), also processed data FPKM, RPKM TPM (see interesting blog-post discussing use data). tools R analysis libraries use downstream analysis rely unnormalized unprocessed data count data. counts simply represent number reads mapped transcript.often hardest part assess; sometimes clear provided data actual raw count values. Usually, can find (sometimes rather cryptic) information clicking sample ID ‘Data processing’ section. paragraph text often mention ‘Supplementary_files_format_and_content’ describes actual content downloadable files. Following number examples valid invalid data provided experiment. single-sample pages can find description. Unfortunately, authors experiment free explain data please, likely encounter exact descriptions:\nValid: Global effects SUPT4H1 RNAi gene expression HEK293 cells; description: Supplementary_files_format_and_content: tab-delimited, gzip-compressed text file matrix raw quantitation results sample. Mentions raw-anything likely valid.\nValid: Transcriptome profiling normal transformed mouse urothelial cells; description: Supplementary_files_format_and_content: raw_counts_unfiltered.txt.gz: Raw counts.. Even though 4 different file types available download, one contains actual raw count data.\nValid: Quantitative Analysis PPARD RNA-Seq Transcriptomes Mouse Gastric Corpus Epithenial Cells Next Generation Sequencing (NGS); description: Supplementary_files_format_and_content: raw counts sequencing reads genes. Single file download containing raw counts.\nValid: Impact shRNA-mediated KLF4 -modulation transcriptome profile human keratinocyte precursor cells; description: Supplementary_files_format_and_content: raw normalized counts.\nValid: Role TAF5L TAF6L controlling mouse ESC state; description: Supplementary_files_format_and_content: Comma-separated text files containing raw counts, normalized counts, FPKM expression values TPM expression values.. common include normalized raw values.\nInvalid: Co-chaperone Mzb1 key effector Blimp1 plasma cell differentiation b1 integrin function; description: Supplementary_files_format_and_content: normalized gene expression values (FPKM) calculated Cuffnorm (v2.2.1) provided TXT format. clearly states normalized data (FPKM values) available.\nInvalid: EGFR-targeted therapy-induced resistance mechanism glioblastoma; description: Supplementary_files_format_and_content: .txt file including gene name, RefSeq id, locus, RPKM value gene sample. RPKM values available normalized raw count values, usable.\nInvalid: RNA-seq subcellular fractions nESC EpiLC; description: Supplementary_files_format_and_content: _genes.fpkm_table.txt: Tab-delimited text files include FPKM values replicate.*. FPKM data available.\nValid: Global effects SUPT4H1 RNAi gene expression HEK293 cells; description: Supplementary_files_format_and_content: tab-delimited, gzip-compressed text file matrix raw quantitation results sample. Mentions raw-anything likely valid.Valid: Transcriptome profiling normal transformed mouse urothelial cells; description: Supplementary_files_format_and_content: raw_counts_unfiltered.txt.gz: Raw counts.. Even though 4 different file types available download, one contains actual raw count data.Valid: Quantitative Analysis PPARD RNA-Seq Transcriptomes Mouse Gastric Corpus Epithenial Cells Next Generation Sequencing (NGS); description: Supplementary_files_format_and_content: raw counts sequencing reads genes. Single file download containing raw counts.Valid: Impact shRNA-mediated KLF4 -modulation transcriptome profile human keratinocyte precursor cells; description: Supplementary_files_format_and_content: raw normalized counts.Valid: Role TAF5L TAF6L controlling mouse ESC state; description: Supplementary_files_format_and_content: Comma-separated text files containing raw counts, normalized counts, FPKM expression values TPM expression values.. common include normalized raw values.Invalid: Co-chaperone Mzb1 key effector Blimp1 plasma cell differentiation b1 integrin function; description: Supplementary_files_format_and_content: normalized gene expression values (FPKM) calculated Cuffnorm (v2.2.1) provided TXT format. clearly states normalized data (FPKM values) available.Invalid: EGFR-targeted therapy-induced resistance mechanism glioblastoma; description: Supplementary_files_format_and_content: .txt file including gene name, RefSeq id, locus, RPKM value gene sample. RPKM values available normalized raw count values, usable.Invalid: RNA-seq subcellular fractions nESC EpiLC; description: Supplementary_files_format_and_content: _genes.fpkm_table.txt: Tab-delimited text files include FPKM values replicate.*. FPKM data available.","code":""},{"path":"a2-data_loading.html","id":"a2-data_loading","chapter":"B Batch Loading Expression Data in R","heading":"B Batch Loading Expression Data in R","text":"code example shows batch-load multiple files containing expression (count) data single sample. data example can found GEO ID GSE109798. Depending tool used convert mapping data per-sample count file, contents files can () different. Downloading data example experiment GEO gives us single .tar file called GSE109798_RAW.tar. Extracting archive file nets us folder following files:","code":"\nfile.names <- list.files('./data/GSE109798_RAW/')"},{"path":"a2-data_loading.html","id":"decompressing","chapter":"B Batch Loading Expression Data in R","heading":"B.0.1 Decompressing","text":"file extension files .txt.gz means files compressed using gzip need unpacked can loaded. easiest method using system gunzip command files can done within R applying gunzip command using system function file.Now can update file.names variable since file name changed.","code":"\n## Change directory to where the files are stored\nsetwd('./data/GSE109798_RAW/')\nsapply(file.names, FUN = function(file.name) {\n  system(paste(\"gunzip\", file.name)) \n  })\nfile.names <- list.files('./data/GSE109798_RAW/')"},{"path":"a2-data_loading.html","id":"determining-data-format","chapter":"B Batch Loading Expression Data in R","heading":"B.0.2 Determining Data Format","text":"Next, can inspect contents files, assuming layout/ column names etc. decide need use analysis.files contain (much) just count value gene can see columns (transcript) length, TPM, FPKM, etc. Also, count-column called expected_count raises questions well.\nexpected count value usable contains \ninformation - compared raw count - actually require. \nexpected part results multimapped reads \nsingle read mapped multiple positions genome. \ntranscript originates one location, multimapped read \nusually discarded. expected count though, instead \ndiscarding read completely estimated \noriginates added fraction count\nvalue. value 1.32 see line 5 \nexample means true count 1 (uniquely mapped\nread) .32 (estimated part) results \nalgorithm can mean multiple things.\n\nmentioned , require integer count data use \npackages DESeq2 edgeR \ntwo methods convert expected count raw count data: +\nround value nearest integer (widely accepted method\nwell within expected sampling variation), + discard \nfraction part using example floor()\nfunction.\n","code":"\n## Call the system 'head' tool to 'peek' inside the file\nsystem(paste0(\"head \", \"./data/GSE109798_RAW/\", file.names[1]))"},{"path":"a2-data_loading.html","id":"loading-data-1","chapter":"B Batch Loading Expression Data in R","heading":"B.0.3 Loading Data","text":"columns want keep transcript_id expected_count columns ignore rest (might interested data later though). need lead file separately can define function reads data, keeping columns interest returning dataframe data. Note first line file used header, check setting header argument TRUE, sometimes expression data starts line 1. file name also used name column dataframe know column sample. done splitting file name (using strsplit) using dot (‘.’) keeping first part (.e. ‘GSM2970156_HCC1806E224A’) discarding second part (‘isoforms.results.txt’). strsplit function however always returns list, case containing vector 5 splitted elements:Applying read_sample function file names gives us set data frames can merge together using merge function. merge data based transcript id defined = 1 argument pointing first column. start reading just one file ‘base’ dataframe merge files.processing seemed data set divided two groups also listed GEO website project:GPL11154 Illumina HiSeq 2000 (Homo sapiens)GPL13112 Illumina HiSeq 2000 (Mus musculus)first 6 files human source last 6 mouse. Therefore, following code shows read first 6 samples merge single dataframe. Repeating process 6 files result another dataframe samples.Table continues belowThe dataset variable now contains data first 6 samples experiment. advisable compare number rows data set number rows single sample. guaranteed samples exactly number genes/transcripts present (.e., 0-values might discarded) results final data set many rows smallest sample. See help merge case argument can used introduce extra rows missing data.","code":"\n## String splitting in R\n## (the fixed = TRUE is required as the dot is a special character, see '?strsplit')\nstrsplit('GSM2970155_HCC1806E224B.isoforms.results.txt.gz', '.', fixed = TRUE)## [[1]]\n## [1] \"GSM2970155_HCC1806E224B\" \"isoforms\"               \n## [3] \"results\"                 \"txt\"                    \n## [5] \"gz\"\n## Keeping the sample identifier\nstrsplit('GSM2970155_HCC1806E224B.isoforms.results.txt.gz', '.', fixed = TRUE)[[1]][1]## [1] \"GSM2970155_HCC1806E224B\"\n## Function for reading in files\nread_sample <- function(file.name) {\n  ## Extract the sample name for naming the column\n  sample.name <- strsplit(file.name, \".\", fixed = TRUE)[[1]][1]\n  ## Read the data, setting the 'transcript_id' as row.names (column 1)\n  sample <- read.table(file.name, header = TRUE, sep=\"\\t\", row.names = NULL)\n  ## Rename the count column\n  names(sample)[5] <- sample.name\n  ## Return a subset containing the 'transcript_id' and sample name columns\n  return(sample[c(1, 5)])\n}\nsetwd('./data/GSE109798_RAW/')\n\n## Read the FIRST sample\ndataset <- read_sample(file.names[1])\n\n## Read first sample group (6)\nfor (file.name in file.names[2:6]) {\n  sample <- read_sample(file.name)\n  dataset <- merge(dataset, sample, by = 1)\n}\n\npander(head(dataset))"},{"path":"a2-annotation.html","id":"a2-annotation","chapter":"C Annotating an RNA-Seq Experiment","heading":"C Annotating an RNA-Seq Experiment","text":"chapter describes annotate data, meaning assigning names functions Differentially Expressed Genes. step can either easy bit challenging depending data source project. first example relevant data set originated NCBI GEO, ‘’ always annotated default. found data set elsewhere (GEO identifier project), skip Manually annotate data section .goal chapter - least - find Gene Symbol common ID (NCBI/ Ensembl) gene. Using information much easier find relevant information online sources say something functionality impact DEGs.Even though might lucky either already data simple GEO query results everything want/ need, tasked also perform manual annotation phase even use found annotation data. written code access sources good excercise using complex Bioconductor packages. Also, can always use find information require later .","code":""},{"path":"a2-annotation.html","id":"download-geo","chapter":"C Annotating an RNA-Seq Experiment","heading":"C.1 Downloading annotation data from GEO","text":"Given GSE**** id experiment can download available (annotation)data GEO using GEOquery library. example experiment used previously, following code downloads two files; series_matrix.txt.gz file also available download NCBI GEO website. Note destdir argument storing known location instead /tmp folder later use.downloaded, can load data using:gse object may contain useful information samples genes. Unfortunately example data set information available regarding samples can access phenoData slot. Slots R objects method storing multiple data objects single R object, case object called ExpressionSet (see ?ExpressionSet description contents structure). storing information experiment handy; can store data frame actual expression data, list containing laboratory information, chunk text abstract, etc. single object.common data sets found ExpressionSet object:\n* @assayData: actual expression data. common microarray experiments () rare RNA-seq experiments\n* @phenoData: sample information, present GEO data set\n* @featureData: data frame holding feature (=gene) information, slightly rare RNA-seq data sets\n* @experimentData: information lab performed experimentYou already know access column data frame using $ notation, accessing complete data frame gse object done using @ symbol:see sample described using sample identifier (GSM****), information sample group (title column) number characteristics. actual information included depends data set.Contents gse@phenoData@data information (selection) (continued )data set properly annotated gene information, accessible @featureData slot. instance, data set GSE20489 annotation available. loaded, see following annotation columns present gene featureData slot (following small subset, total can well 20 columns information per gene). Also, metadata listing shows many features (=genes, 54675) samples (54) expression data present (@assayData slot):Gene information stored ExpressionSet object experiment GSE20489 (continued )","code":"\nlibrary(GEOquery)\ngetGEO(GEO = 'GSE101942', destdir = \"./data/\")\nlibrary(GEOquery)\ngse <- getGEO(filename = './data/GSE101942_series_matrix.txt.gz')\n# Access a data frame containing phenotype information\ngse@phenoData@data\n# This gets a list with a single 'ExpressionSet' object\nGSE20489 <- getGEO(GEO = \"GSE20489\", destdir = './data')\n# Print some metadata of the featureData slot\nGSE20489[[1]]## ExpressionSet (storageMode: lockedEnvironment)\n## assayData: 54675 features, 54 samples \n##   element names: exprs \n## protocolData: none\n## phenoData\n##   sampleNames: GSM514737 GSM514738 ... GSM514790 (54 total)\n##   varLabels: title geo_accession ... tissue:ch1 (39 total)\n##   varMetadata: labelDescription\n## featureData\n##   featureNames: 1007_s_at 1053_at ... AFFX-TrpnX-M_at (54675 total)\n##   fvarLabels: ID GB_ACC ... Gene Ontology Molecular Function (16 total)\n##   fvarMetadata: Column Description labelDescription\n## experimentData: use 'experimentData(object)'\n##   pubMedIds: 23883607 \n## Annotation: GPL570"},{"path":"a2-annotation.html","id":"manually-annotate","chapter":"C Annotating an RNA-Seq Experiment","heading":"C.2 Manual Data Annotation","text":"Luckily us, R offers number libraries automatically retrieve information large sets genes, unless chosen organism much data available. section demonstrates use two libraries, starting AnnotationDbi followed biomaRt.AnnotationDbi library downloads local copy organism-specific database gene information biomaRt uses online databases retrieve data given query. biomaRt offers far data (1000 data fields per organism) complex use.Since biomaRt also relying online databases might good strategy annotate genes interest (DEGs) instead querying >20.000 genes might retain 20 statistical analysis. planning use biomaRt, skip Discovering Differentialy Expressed Genes (DEGs) chapter first return Using R Bioconductors biomaRt section annotate data.","code":""},{"path":"a2-annotation.html","id":"annotationDBI","chapter":"C Annotating an RNA-Seq Experiment","heading":"C.2.1 Using AnnotationDBI","text":"AnnotationDbi offers data sets many organisms form installable libraries depending experiment need find proper library. example experiment contains samples house mouse (Mus musculus) therefore select data set.following information types available database (use columns function inspect).Available fields databaseThe table shows data available information types given randomly chosen EntrezID ‘1080’. Note table split show 24 data types values.Example values field ‘org.Hs.eg.db’Retrieving data locally stored annotation database can done using mapIds function takes number arguments:x: local database querykeys: IDs data set, often ENSEMBL IDs gene SYMBOLS. example data rownames data set contain gene SYMBOLS (see Table 1 page 2)column: data column retrieve databasekeytype: type data provide. example data gene SYMBOLmultiVals: number columns contain one entry single gene, case want store first one.output mapIds function single character vector can add data frame containing DEGs (case DEGs IPSC trisomic vs disomic comparison containing 828 genes stored ipsc.tri.vs.di data frame).Unfortunately, organisms offer access information information always complete. instance, following table shows number available records Ensembl, Entrez en KEGG IDs downloaded:Statistics annotation columns using AnnotationDbiFor data sets might luck, otherwise continue biomaRt method explained .","code":"\n# Load the AnnotationDbi interface library\nlibrary(AnnotationDbi)\n# Load the Bioconductor installation library (contains 'biocLite()')\nlibrary(BiocInstaller)\n\n# Install and load the organism specific gene database\n# 'org' for Organism\n# 'Hs' for Homo sapiens\n# 'eg' for Entrez Gene IDs\n# (try to load (using the library function) first before installing, it might already be present)\nbiocLite('org.Hs.eg.db')\nlibrary(org.Hs.eg.db)\n.Hs.data.table <- cbind(.Hs.data[1:12,], .Hs.data[13:24,])\ncolnames(.Hs.data.table) <- c('Data Type', 'Example Value',\n                              'Data Type', 'Example Value')\npander(.Hs.data.table, caption=\"Example values for each field in 'org.Hs.eg.db'\")\n# Retrieve the ENSEMBL gene ID, e.g. 'ENSG00000001626'\nipsc.tri.vs.di$Ensembl <- mapIds(x = org.Hs.eg.db,\n                           keys=row.names(ipsc.tri.vs.di),\n                           column=\"ENSEMBL\",\n                           keytype=\"SYMBOL\",\n                           multiVals=\"first\")\n\n# Retrieve the ENTREZ gene ID, e.g. '12323'\nipsc.tri.vs.di$EntrezID <- mapIds(x = org.Hs.eg.db,\n                            keys=row.names(ipsc.tri.vs.di),\n                            column=\"ENTREZID\",\n                            keytype=\"SYMBOL\",\n                            multiVals=\"first\")\n\n# Retrieve the KEGG enzyme code, e.g. 2.7.11.17\nipsc.tri.vs.di$Enzyme <- mapIds(x = org.Hs.eg.db,\n                           keys=row.names(ipsc.tri.vs.di),\n                           column=\"ENZYME\",\n                           keytype=\"SYMBOL\",\n                           multiVals=\"first\")"},{"path":"a2-annotation.html","id":"using-biomart","chapter":"C Annotating an RNA-Seq Experiment","heading":"C.2.2 Using biomaRt","text":"short explanation biomart (source: Wikipedia):“purpose BioMarts Ensembl Genomes allow user mine download tables containing genes single species, genes specific region chromosome genes one region chromosome associated InterPro domain. BioMarts also include filters refine data extracted attributes (Variant ID, Chromosome name, Ensembl ID, location, etc.) appear final table file can selected user.”text mentiones species, attributes filters, need combine elements query biomart databases annotation. following code ‘chunks’ show possible values elements gather store relevant data.biomaRt library interfaces biomart.org online database. Sometimes biomart website (offers browsable database) due maintenance, many mirrors can still used, instance Ensembl.required explore number objects given example code likely suitable data/ organism. project well documented Bioconductor biomaRt website links biomaRt users guide.Using listDatasets function can get full list available datasets. Store list R object ‘browse’ object RStudio see organism included. Copy name dataset, species element use.Subset 214 datasets availableNext deciding filter. can use listFilters function ensembl object, storing full list filters. wise view RStudio find filter use.filter specifies use search . instance, AnnotationDbi queries gave us Ensembl gene ID’s use ensembl_gene_id filter. viewing list filters RStudio can use search text-box top-right view, example ensembl look filters use type identifier.Subset 437 available filtersFinally get point determine data like get database. attributes can get listAttributes function ensembl object. - especially attributes since often >1000 selectable options - store attributes view RStudio look data want.Subset 3172 available attributesIf want gene chromosome, start- end-position well description (note, just example, , interesting information available !) combine character vector (attrs.get, see ). one caveat though, order get back results input order! means simple combine data original data set need merge together. However, able merge data need know record belongs gene therefore add selected filter (case ensembl_gene_id) list attributes get. , get results dataframe can use merge function R combine gene information actual data:now three needed elements (species, filter attributes) ready query database getBM function (read help using ?getBM).example retrieves data using set five Ensembl gene IDs since genes Ensembl ID ’ve seen , filter first use values parameter .results object data.frame 5 columns can merge data set giving us following annotation columns (combined AnnotationDBI biomaRt libraries).just example use biomaRt library comes selecting correct filter looking interesting attributes retrieve. information can found documentation avaialble `vignette(‘biomaRt’)","code":"\n# Load the library\nlibrary(biomaRt)\n\n# Use an alternative database server as the regular one sometimes has issues..\nensembl=useMart(\"ENSEMBL_MART_ENSEMBL\", host=\"https://www.ensembl.org\")\n# Select the 'ensembl' database\nensembl <- useMart(\"ensembl\")\nmart.datasets <- listDatasets(ensembl)\n\n# Select the correct dataset, for the example data we select the 'hsapiens_gene_ensembl' \nensembl <- useDataset('hsapiens_gene_ensembl', mart = ensembl)\nfilters <- listFilters(ensembl)\nattributes <- listAttributes(ensembl)\nmerge(x = ipsc.tri.vs.di, y = results, by.x = 'Ensembl', by.y = 'ensembl_gene_id')\n# Set the 'attributes' values\nattrs.get <- c(\"ensembl_gene_id\", \"chromosome_name\", \n               \"start_position\",\"end_position\", \"description\")\n\n# Perform a biomaRt query using 'getBM'\nresults <- getBM(attributes = attrs.get,\n                 filters = \"ensembl_gene_id\",\n                 values = ipsc.tri.vs.di$Ensembl[1:5], \n                 mart = ensembl)\n\nresults$gene_length <- abs(results$end_position - results$start_position)"}]
diff --git a/index.Rmd b/index.Rmd
index 8cab780..9b3b931 100644
--- a/index.Rmd
+++ b/index.Rmd
@@ -23,6 +23,7 @@ For this Genomics analysis we are going to look at variants between the chosen d
 
 ### Project Deliverables
 
+* **Project Proposal**: A short written proposal of the project clearly describing the data set, the research question and the methods + tools to be used.
 * **Genomics pipeline**: The genomics analysis pipeline is well documented in one or more *lab journals* and contains the correct tools based on the research question and data set.
 * **Presentation**: The research results will be presented in an attractive method and proper level given the target audience.
 
@@ -61,7 +62,7 @@ Finding a suitable data set can be challenging as there are a few requirements s
 
 Data is available from the NCBI [Sequence Read Archive](https://www.ncbi.nlm.nih.gov/sra) and Gene Expression Omnibus.
 
-The proper way to search for a data set suitable for both Genomics and Transcriptomics projects is by creating a query that searches for both: [("expression profiling by high throughput sequencing"[DataSet Type]) AND "genome variation profiling by high throughput sequencing"[DataSet Type]](https://www.ncbi.nlm.nih.gov/gds/?term=(%22expression+profiling+by+high+throughput+sequencing%22%5BDataSet+Type%5D)+AND+%22genome+variation+profiling+by+high+throughput+sequencing%22%5BDataSet+Type%5D).
+The proper way to search for a data set suitable for both Genomics and Transcriptomics projects is by creating a query that searches for both: [("expression profiling by high throughput sequencing"[DataSet Type]) AND "genome variation profiling by high throughput sequencing"[DataSet Type]](https://www.ncbi.nlm.nih.gov/gds/?term=(%22expression+profiling+by+high+throughput+sequencing%22%5BDataSet+Type%5D)+AND+%22genome+variation+profiling+by+high+throughput+sequencing%22%5BDataSet+Type%5D). If this selection actually contains anything of interest is up to you; if you cannot find something interesting please use a different search strategy. One could be to search for RNA-Seq experiments in the GEO where one of the groups is preferably a mutant. 
 
 Note that for the genomics part, not all available samples have to be processed. Other than for transcriptomics, it is not common that there is a case/control experimental setup but rather a single condition of group with possibly a few replicates. If the data set does include many samples, it is recommended to select a subset of samples to process (do discuss the reasoning behind this selection).
 
@@ -73,6 +74,7 @@ While there is no golden standard for analyzing RNA-sequencing data sets as ther
 
 ### Project Deliverables
 
+* **Project Proposal**: A short written proposal of the project clearly describing the data set, the research question and the methods + libraries to be used.
 * **Lab Journal**: The lab journal is readable, complete and aimed at the reproducability of the research. Relevant visualisations are include and statistical methods performed to support the biological interpretation of the results.
 
 ## *Lab Journal*