```{r setup, include=FALSE}
library(magrittr)
library(tidyverse)
library(vroom)

library(DESeq2)
library(ggdendro)

knitr::opts_chunk$set(echo = TRUE)
```

## 2021-07-20: Getting Started with RNA-seq by Reading in the Data

Please email me your Rmarkdown document at <kkeith@coriell.org> with your **LAST NAME** in the filename by 9AM Wednesday 7/21. **THIS IS DUE TOMORROW, I NEED TO KNOW YOU CAN READ THE DATA IN BEFORE WE KEEP GOING**

---

### Lecture Questions

1. What is the minimum number of samples necessary for RNA-seq?

*WRITE YOUR ANSWER HERE*

3

---

2. When would you need to use ribosomal RNA depletion to make an RNA-seq library instead of using poly-A capture?

*WRITE YOUR ANSWER HERE*

When you want to study non-coding RNAs

---

3. If you have a gene covered by only 5 reads, what should you do with it?

*WRITE YOUR ANSWER HERE*

Remove it/filter it out

---

4. What are 3 qualities of RNA you can analyze using RNA-seq data?

*WRITE YOUR ANSWER HERE*

Change in level of RNA transcripts either at the whole gene or the transcript level, alternative splicing, mutations, non-coding RNAs

---

5. Give one advantage and one disadvantage of single cell RNA-seq compared to bulk RNA-seq.

*WRITE YOUR ANSWER HERE*

ADVANTAGES: Get measurements in many cells, not an average of a cell population, diversity of expression, depending on how it's done can have the spacial location in the tissue
DISADVANTAGES: What cell type is it? Currently can be hard to label the cell type with subject matter expertise, but this is being worked on. Data sparsity, measurement uncertainty, very expensive

---

<br>

### RNA-seq with Full Data Set

#### Read / Wrangle Data

Taking the data in the zipped folder `full_count_files` which contains the full RNA-seq count tables for the human genome (not just the chr21 files we were using for the demo) and use the chunk below to read it in

```{r}
### use this chunk to read the files in

# get list of files
list.files(path = '../full_count_files',
           pattern = '.txt$',
           full.names = T) -> files

### read files in recursively with vroom()
vroom::vroom(files, id = 'file_path', delim = '\t', comment = '#', skip = 2, 
             col_names = c('ensembl_gene_id', 'chr', 'start', 'end', 
                           'strand', 'length', 'count')) %>%
  mutate(sample_id = str_extract(file_path, '[d,s][A-z]{2,3}[0-9]'),
         treatment = str_remove_all(sample_id, '[0-9]')) %>%
  dplyr::select(-file_path) -> data
```

Use this chunk to transform the count tables we just read in into a DESeq object.

```{r}
# transform the data table into a count matrix
data %>%
  dplyr::select(ensembl_gene_id, sample_id, count) %>%
  spread(sample_id, count) %>% 
  as.data.frame() -> count_matrix

# get a metadata table
data %>%
  dplyr::select(sample_id, treatment) %>%
  distinct() -> metadata

# give count matrix and metadata table and transform into a DESeq object
dds <- DESeqDataSetFromMatrix(countData = count_matrix,
                              colData = metadata,
                              tidy = TRUE,
                              design = ~ treatment)
```

Finally, in the last wrangling step, filter the data. Make sure to remove **BOTH** rows that only have zero counts and to filter for genes where at least 3 samples have a count of 10 transcripts or more

```{r}
### filter out rows that contain only zero counts
keep <- rowSums(counts(dds)) >= 1
dds <- dds[keep, ]

# Filter for genes where at least three samples have a count of 10 or higher
keep <- rowSums(counts(dds) >= 10) >= 3
dds <- dds[keep, ]

### OPTIONAL, but recommended: check how many genes are filtered out
table(keep)
```

<br>

















