rna-seq-survey-MS-v2021_02_22.Rmd

---
title: "`r gsub('.Rmd', ' ', gsub('_', ' ', knitr::current_input()))`"
author: "`r Sys.getenv('USER')`"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  html_document:
    toc: true
    toc_float: true
---


# Libraries
```{r rna-seq-survey-MS-v20-1 }

library(tidyverse)
library(janitor)
library(knitr)
library(cowplot) # for plot_grid and draw_label
library(RColorBrewer)
library(pander) # to make session info
library(kableExtra)
library(snakecase)
library(corrr)
```

# Config
```{r rna-seq-survey-MS-v20-2 }

inputs_dir <- "inputs"
file_pattern <- "survey_data"
data_file <- tibble(files = list.files(inputs_dir, file_pattern)) %>%
  arrange(desc(files)) %>%
  top_n(1, wt = files) %>%
  pull(files)

print(data_file)

```

# Settings
```{r rna-seq-survey-MS-v20-3 }

min_genes_gt_0 <- 100

```

# Functions
```{r rna-seq-survey-MS-v20-4 }

is_outlier <- function(x) {
  (x > quantile(x, 0.75) + 1.5*IQR(x)) |
    (x < quantile(x, 0.25) - 1.5*IQR(x))
}


most_common_val_w_percent <- function(value_vector = c("a", "a", "b")){
  mc_vals <- tabyl(value_vector) %>%
    arrange(desc(n)) %>%
    top_n(1, n) 
  paste0(round(100*mc_vals$percent), "% ", mc_vals$value_vector)
}

# StatBin2 allows depiction of empty bins as blank instead of a horizontal line:
# https://stackoverflow.com/questions/57128090/remove-baseline-color-for-geom-histogram
StatBin2 <- ggproto(
  "StatBin2", 
  StatBin,
  compute_group = function (data, scales, binwidth = NULL, bins = NULL, 
                            center = NULL, boundary = NULL, 
                            closed = c("right", "left"), pad = FALSE, 
                            breaks = NULL, origin = NULL, right = NULL, 
                            drop = NULL, width = NULL) {
    if (!is.null(breaks)) {
      if (!scales$x$is_discrete()) {
        breaks <- scales$x$transform(breaks)
      }
      bins <- ggplot2:::bin_breaks(breaks, closed)
    }
    else if (!is.null(binwidth)) {
      if (is.function(binwidth)) {
        binwidth <- binwidth(data$x)
      }
      bins <- ggplot2:::bin_breaks_width(scales$x$dimension(), binwidth, 
                                         center = center, boundary = boundary, 
                                         closed = closed)
    }
    else {
      bins <- ggplot2:::bin_breaks_bins(scales$x$dimension(), bins, 
                                        center = center, boundary = boundary, 
                                        closed = closed)
    }
    res <- ggplot2:::bin_vector(data$x, bins, weight = data$weight, pad = pad)

    # drop 0-count bins completely before returning the dataframe
    res <- res[res$count > 0, ] 

    res
  })

```

# Import data
```{r rna-seq-survey-MS-v20-5 }

col_spec <- cols(
  .default = col_double(),
  dataset_id = col_character(),
  disease = col_character(),
  pedaya = col_character(),
  race = col_character(),
  ethnicity = col_character(),
  gender = col_character(),
  cohort_name = col_character(),
  cohort_code = col_character(),
  doi_link = col_character(),
  source_repository = col_character(),
  repository_cohort_accession = col_character(),
  repository_dataset_accession = col_character()
)

counts_and_meta_raw <- read_tsv(file.path(inputs_dir, data_file), col_types = col_spec)

counts_and_meta <- counts_and_meta_raw %>%
  group_by(cohort_code) %>%
  mutate(n_in_cohort = n())

write_csv(counts_and_meta, "results/Table_S1.csv")

```


# Colors
```{r rna-seq-survey-MS-v20-6 }

# brewer.pal(8, "Paired")

these_colors <- c("#A6CEE3", "#1F78B4", "#B2DF8A", "#33A02C", "#FB9A99", "#E31A1C", 
"#FDBF6F", "#FF7F00")

# light blue, dark blue, etc: green, red, orange (then purple and brown)

these_colors <- c("#1F78B4", "#1F78B4", "#33A02C", "#33A02C", "#E31A1C", "#E31A1C", 
"#FF7F00", "#FF7F00")

names(these_colors) <-  c("not assigned", "Total reads", "Non-exonic reads",  "Exonic reads", "Duplicate reads",  "Non-duplicate reads",  "Unmapped reads", "Mapped reads")

these_colors_with_MEND = c(these_colors, "MEND reads"="#000000")

```


# Datasets analyzed
```{r rna-seq-survey-MS-v20-7 }

nrow(counts_and_meta)

```

# Diseases in dataset

```{r rna-seq-survey-MS-v20-8 }
counts_and_meta %>% tabyl(disease)

disease_freq_table <- counts_and_meta$disease %>%
  str_to_sentence %>%
  str_replace(" nos$", " NOS") %>%
  as.factor %>%
  fct_infreq %>%
  fct_lump(prop=0.01) %>%
  tabyl (var1 = "Disease") %>%
  adorn_pct_formatting(digits = 1)

colnames(disease_freq_table)[1]="Disease"

write_tsv(disease_freq_table, "results/disease_freq_table.tsv")
write_csv(disease_freq_table, "results/disease_freq_table.csv")

disease_freq_table  %>%
  kable %>% 
  kable_styling(bootstrap_options = "striped", full_width = F)

```

# Cohorts in dataset

```{r rna-seq-survey-MS-v20-9, fig.height=6, fig.width=5}

old_cohort_size_histogram <- ggplot(counts_and_meta) + 
  geom_histogram(aes(n_in_cohort, group=cohort_code), 
                 fill = "lightgrey", color = "black", stat = StatBin2) +
  theme_minimal(base_size = 20) +
  theme(axis.line.x = element_line(color="darkgrey", size = 0.5),
        axis.line.y = element_line(color="darkgrey", size = 0.5)) +  xlab("Cohort size") +
  ylab("Cohorts") + 
  theme(legend.position="none")


# max(table(counts_and_meta$n_in_cohort)/as.numeric(names(table(counts_and_meta$n_in_cohort))))
cohort_hist_bin_size <- 10


cohort_size_histogram_raw <- ggplot(counts_and_meta %>%
  select(cohort_code, n_in_cohort) %>%
  distinct ) + 
  geom_histogram(aes(n_in_cohort), 
                 fill = "lightgrey", color = "black", stat = StatBin2,
                 breaks = seq(0,500, by=10)) +
  theme_minimal(base_size = 20) +
  theme(axis.line.x = element_line(color="darkgrey", size = 0.5),
        axis.line.y = element_line(color="darkgrey", size = 0.5)) +  
  xlab("Cohort size") +
  theme(legend.position="none")

max_y_val <- max(ggplot_build(cohort_size_histogram_raw)$data[[1]][,"y"])


cohort_size_histogram <- cohort_size_histogram_raw + 
  scale_y_continuous("Cohorts", breaks = seq(0, ceiling(max_y_val/2) * 2, 2)) 
  

cohort_size_histogram

n_distinct(counts_and_meta$cohort_code)
sort(unique(counts_and_meta$n_in_cohort))

counts_and_meta %>%
  select(cohort_code, n_in_cohort) %>%
  distinct %>%
  arrange(desc(n_in_cohort)) %>%
  kable

counts_and_meta %>%
  select(cohort_code, n_in_cohort) %>%
  distinct %>%
  pull(n_in_cohort) %>%
  median

```

# Table S2 - repositories
```{r rna-seq-survey-MS-v20-10}

repo_info <- read_tsv("inputs/repos.tsv")

repo_table <- repo_info %>%
  filter(repo_abbrev %in% counts_and_meta$source_repository) %>%
  rename(`Name`=repo_name, `Abbreviation`=repo_abbrev, URL = repo_url) %>%
  select(-source_accession_pattern) 

write_tsv(repo_table, "results/Table_S2_repo_table.tsv")
write_csv(repo_table, "results/Table_S2_repo_table.csv")

repo_table %>%
  kable %>% 
  kable_styling(bootstrap_options = "striped", full_width = F)

```

# Table S3 - data access table
```{r rna-seq-survey-MS-v20-11}
cohort_table <- counts_and_meta %>% 
  select(cohort_code, cohort_name, source_repository, repository_cohort_accession, n_in_cohort) %>% # removed DOI link
#  group_by(across(c(-disease)))
  distinct %>%
  arrange(desc(n_in_cohort)) 

colnames(cohort_table) <- to_any_case(colnames(cohort_table), case = "sentence")

write_tsv(cohort_table, "results/Table_S3_cohort_table.tsv")
write_csv(cohort_table, "results/Table_S3_cohort_table.csv")




cohort_table %>% 
  kable %>% 
  kable_styling(bootstrap_options = "striped", full_width = F)
```




# Race and gender in dataset

```{r rna-seq-survey-MS-v20-12 }

counts_and_meta %>%
  tabyl(race)


counts_and_meta %>%
  select(race) %>%
  na.omit() %>%
  tabyl(race) %>%
  adorn_totals("row")


counts_and_meta %>%
    select(ethnicity) %>%
  na.omit() %>%
tabyl(ethnicity) %>%
  adorn_totals("row")


counts_and_meta %>%
  filter(ethnicity == "Hispanic or Latino") %>%
  tabyl(race)


counts_and_meta %>%
  filter(ethnicity == "Hispanic or Latino") %>%
  select(race) %>%
  na.omit() %>%
tabyl(race) %>%
  adorn_totals("row")
```


# Genders in dataset

```{r rna-seq-survey-MS-v20-13 }


table(counts_and_meta$gender, useNA = "always")

# samples with reported genders
counts_and_meta %>% 
  filter(! is.na(gender)) %>%
  filter(! gender %in% c("not reported", "unknown")) %>%
  tabyl(gender) %>%
  adorn_totals


```



# Ages in dataset

```{r rna-seq-survey-MS-v20-14 }


table(counts_and_meta$pedaya, useNA = "always")

# at least one target sample is >30

n_ped_aya <- sum(counts_and_meta$pedaya=="Yes, age < 30 years" | counts_and_meta$age_at_dx<30, na.rm = TRUE)
n_ped_aya
n_ped_aya/nrow(counts_and_meta)
# more than 95% of which are pediatric <30; of the Datasets with specific ages at diagnosis reported, the median was 7

median(counts_and_meta$age_at_dx, na.rm = TRUE)

table(is.na(counts_and_meta$age_at_dx))

```
# Read lengths in dataset
```{r rna-seq-survey-MS-v20-15 }

seq_len_round_25 <- round(counts_and_meta$seq_length / 25) * 25
table(seq_len_round_25, useNA = "always")
tabyl(seq_len_round_25)


table(is.na(counts_and_meta$seq_length))

table(counts_and_meta$seq_length)
summary(counts_and_meta$seq_length)
IQR(counts_and_meta$seq_length, na.rm = TRUE)

 old_read_length_histogram <-  ggplot(counts_and_meta) + 
  geom_histogram(aes(seq_length, group=cohort_code), 
                 fill = "lightgrey", color = "black", stat = StatBin2) +
      theme_minimal(base_size = 20) +
    theme(axis.line.x = element_line(color="darkgrey", size = 0.5),
        axis.line.y = element_line(color="darkgrey", size = 0.5)) +
  xlab("Sequence length") +
    ylab("Datasets")
 
 read_length_histogram <-  ggplot(counts_and_meta) + 
  geom_histogram(aes(seq_length), 
                 fill = "lightgrey", color = "black", stat = StatBin2,
                 binwidth = 1) +
      theme_minimal(base_size = 20) +
    theme(axis.line.x = element_line(color="darkgrey", size = 0.5),
        axis.line.y = element_line(color="darkgrey", size = 0.5)) +
  xlab("Sequence length") +
    ylab("Datasets") +
   expand_limits(x=0)
  
read_length_histogram

```
# Figure 1
```{r rna-seq-survey-MS-v20-16, fig.height=6, fig.width=5}



figure_name <- "f1"



f1 <- plot_grid(cohort_size_histogram,
read_length_histogram,
nrow=2,
labels = c("A", "B")) +
   draw_label(paste(figure_name, Sys.time()),
              x = 0, y = 0.02, hjust = 0, size = 6) 

f1
ggsave("results/Figure_1.png", f1, height=6, width=5)
ggsave("results/Figure_1.pdf", f1, height=6, width=5)

```



# Check for counting errors
```{r rna-seq-survey-MS-v20-17 }

if(nrow(subset(counts_and_meta, Mapped > Total_reads)) > 0) stop("Counting error")
if(nrow(subset(counts_and_meta, MND > Mapped)) > 0) stop("Counting error")
if(nrow(subset(counts_and_meta, MEND > MND)) > 0) stop("Counting error")

```


# Summarize total read counts
```{r rna-seq-survey-MS-v20-18 }

min(counts_and_meta$Total_reads)

Total_read_counts <- counts_and_meta %>% 
  ungroup %>%
  mutate(dataset_value = Total_reads/1e6) %>%
  summarize(variance= var(dataset_value),
            min = min(dataset_value),
            max = max(dataset_value),
            median = median(dataset_value),
            p25 = quantile(dataset_value, 0.25),
            p75 = quantile(dataset_value, 0.75),
            iqr = IQR(dataset_value))

kable(Total_read_counts %>%
        select(-variance), digits = 1)

kable(Total_read_counts %>%
        select(-variance), digits = 0)


total_read_count_histogram <-  ggplot(counts_and_meta) + 
  geom_histogram(aes(Total_reads/1E6), 
                 fill = "lightgrey", color = "black", stat = StatBin2) +
      theme_minimal(base_size = 20) +
    theme(axis.line.x = element_line(color="darkgrey", size = 0.5),
        axis.line.y = element_line(color="darkgrey", size = 0.5)) +
  xlab("Total reads") +
    ylab("Datasets")

total_read_count_histogram
```

# Calculate percentages


```{r rna-seq-survey-MS-v20-19 }

 read_counts_with_read_type_fractions <- counts_and_meta  %>%
  arrange(desc(Total_reads)) %>%
  mutate(
    pct_unmapped_of_total = (Total_reads - Mapped) / Total_reads,
    pct_dupe_of_mapped = (Mapped - MND) / Mapped,
    pct_non_exonic_of_non_dupe = (MND - MEND) / MND
  )

read_type_fractions_long <- read_counts_with_read_type_fractions %>%
  pivot_longer(cols = starts_with("pct_"), names_to = "read_type_fraction_name", values_to = "dataset_value") 


```

# Summarize read type fractions
```{r rna-seq-survey-MS-v20-20 }

comparison_of_distributions <- read_type_fractions_long %>% 
  group_by(read_type_fraction_name) %>%
  summarize(variance= var(dataset_value),
            min = min(dataset_value),
            max = max(dataset_value),
            median = median(dataset_value),
            p25 = quantile(dataset_value, 0.25),
            p75 = quantile(dataset_value, 0.75),
            iqr = IQR(dataset_value))

kable(comparison_of_distributions %>%
      mutate_if(is.double, ~100*.), digits = 3) %>% 
  kable_styling(bootstrap_options = "striped", full_width = T)


kable(comparison_of_distributions %>%
        select(-variance, -p25, -p75) %>%
      mutate_if(is.double, ~100*.), digits = 0) %>% 
  kable_styling(bootstrap_options = "striped", full_width = T)


kable(comparison_of_distributions %>%
        select(-variance) %>%
      mutate_if(is.double, ~100*.), digits = 0) %>% 
  kable_styling(bootstrap_options = "striped", full_width = T)


```

# statements about read type fractions
See comparison_of_distributions table above

## unmapped
```{r rna-seq-survey-MS-v20-21 }
# in 77 datasets, more than 25% of reads are unmapped
table(read_counts_with_read_type_fractions$pct_unmapped_of_total>0.25)
```


## duplicate
```{r rna-seq-survey-MS-v20-22 }

# 426 datasets have more than 50% duplicates (Fig 2A).
read_counts_with_read_type_fractions %>%
  mutate(above_50 = pct_dupe_of_mapped>0.50) %>%
  tabyl(above_50)

fiftypct <- read_counts_with_read_type_fractions %>%
  group_by(cohort_code) %>%
  summarize(has_a_dataset_with_more_than_50pct_dupes = any(pct_dupe_of_mapped>0.50),
            median_pct_dupe_of_mapped = median(pct_dupe_of_mapped)) %>%
  mutate(median_lt_50 = median_pct_dupe_of_mapped < 0.5,
         median_lt_50_and_has_a_dataset = median_lt_50 & has_a_dataset_with_more_than_50pct_dupes)


fiftypct %>% tabyl(median_lt_50)
fiftypct %>% tabyl(has_a_dataset_with_more_than_50pct_dupes)
fiftypct %>% tabyl(median_lt_50_and_has_a_dataset)
```

## exonic
```{r rna-seq-survey-MS-v20-23 }
read_counts_with_read_type_fractions %>%
  mutate(above_50 = pct_non_exonic_of_non_dupe>0.50) %>%
  tabyl(above_50)

fiftypcte <- read_counts_with_read_type_fractions %>%
  group_by(cohort_code) %>%
  summarize(has_a_dataset_with_more_than_50pct_ne = any(pct_non_exonic_of_non_dupe>0.50),
            median_pct_dupe_of_mapped = median(pct_non_exonic_of_non_dupe)) %>%
  mutate(median_lt_50 = median_pct_dupe_of_mapped < 0.5,
         median_lt_50_and_has_a_dataset = median_lt_50 & has_a_dataset_with_more_than_50pct_ne)


fiftypcte %>% tabyl(median_lt_50)
fiftypcte %>% tabyl(has_a_dataset_with_more_than_50pct_ne)
fiftypcte %>% tabyl(median_lt_50_and_has_a_dataset)
```


## Correlation of depth of sequencing to duplicate read fraction
### Does the depth of sequencing correlate to the number of duplicates? How much variance does that explain?
```{r rna-seq-survey-MS-v20-24 }
cor_tot_reads_and_dupes <- cor(read_counts_with_read_type_fractions$Total_reads, 
    read_counts_with_read_type_fractions$pct_dupe_of_mapped,
    method = "spearman")
round(cor_tot_reads_and_dupes, 2)
round(cor_tot_reads_and_dupes^2, 2)


```
# Figure 4
## Plot correlation of depth of sequencing to duplicate read fraction
```{r rna-seq-survey-MS-v20-25 }

this_data <- read_counts_with_read_type_fractions %>%
  select(Total_reads, pct_dupe_of_mapped)

this_plot_title <- "% Duplicates vs Total_reads"

t2 <- cor.test(this_data$Total_reads, this_data$pct_dupe_of_mapped, 
    alternative = "greater", # "greater" corresponds to positive association, 
    method = "spearman")

label_data <- tibble(
  label=paste0("r=", round(t2$estimate,2), "; p=", signif(t2$p.value, digits=1)),
  x=max(this_data$Total_reads/1E6),
  y=max(this_data$pct_dupe_of_mapped))

f4 <- ggplot(this_data,
       aes(x=Total_reads/1E6, y=pct_dupe_of_mapped)) + 
  geom_point(alpha = 0.5, shape = 20, size =0.1) +
  # geom_smooth(method = "lm", color = "black", se = FALSE) +
  geom_label(data=label_data,
             aes(label=label,
             x=x,
             y=y),
             hjust = 1, vjust = 1
             ) +
    theme_minimal(base_size = 20) +
    theme(axis.line.x = element_line(color="darkgrey", size = 0.5),
        axis.line.y = element_line(color="darkgrey", size = 0.5)) +
  ggtitle(this_plot_title)  + 
  theme(legend.position="none", plot.title=element_text(size=22)) +
  xlab("Read counts (million)") +
  ylab("% Duplicates")

f4

ggsave("results/Figure_4.png", f4 + ggtitle(""))
ggsave("results/Figure_4.pdf", f4 + ggtitle(""))

```



# Generate schematic defining read types
```{r rna-seq-survey-MS-v20-26, fig.height=2, fig.width=6}

colors_for_read_types <- these_colors_with_MEND
names(colors_for_read_types) <- gsub(" reads", "", names(these_colors_with_MEND))

colors_for_read_types <- c(colors_for_read_types, 
                           "Genome" = "darkblue",
                           "Exon" = "darkblue",
                           "Mapped reads" = "darkgrey")

plot1_colnames <- c("label", "x1", "x2", "y1", "y2")

read_type_schematic_data_as_vector <- c("Duplicate", 3, 4, 5, 6, 
                       "MEND", 3, 4, 4, 5, 
                       "MEND", 2.5, 3.5, 3, 4, 
                       "Non-exonic", 4.5, 5.5, 3, 4,
                       "Unmapped", 6, 7, 3, 4,
                       "Exon", 2.5, 4, 1, 2,
                       "Mapped reads", 2.8, 2.8, 6, 7)

rows_in_schematic <- length(read_type_schematic_data_as_vector)/length(plot1_colnames)

read_type_schematic_data <- matrix(read_type_schematic_data_as_vector, 
                     ncol = length(plot1_colnames), byrow = TRUE,
                     dimnames = list(c(1:rows_in_schematic), plot1_colnames)) %>%
  as_tibble %>%
  type_convert() %>%
  mutate(
    base_color = colors_for_read_types[match(label, names(colors_for_read_types))],
    border_color = case_when(
      label %in% c("Mapped reads", "Genome", "MEND", "Exon") ~ NA_character_,
      TRUE ~ base_color),
    fill_color = ifelse(label %in% c("MEND", "Exon"),  base_color, "white"),
    text_color = case_when(
      label %in% c("MEND", "Exon")~  "white",
      TRUE ~ base_color),
    mid_bar_x = map2_dbl(x1, x2,  function(x, y) mean(c(x,y))),
    mid_bar_y = map2_dbl(y1, y2,  function(x, y) mean(c(x,y)))
  )



padding_size = 0.35
end_of_box <- sort(read_type_schematic_data$x2,partial=length(read_type_schematic_data$x2)-1)[length(read_type_schematic_data$x2)-1] + padding_size
read_type_schematic <- ggplot(read_type_schematic_data, 
                              aes(xmin=x1, xmax=x2, ymin=y1, ymax = y2, 
                                  fill = fill_color, color = border_color)) +
  geom_rect() +
  
  geom_segment(x=min(read_type_schematic_data$x1-padding_size), y=1.5, xend = 5.75, yend = 1.5, color = "darkblue") +   
  geom_text(aes(label = label, x = mid_bar_x, y=mid_bar_y,
                color = text_color), 
            size = 4,
            fontface = "bold") +
  geom_rect(aes(xmin = min(x1 - padding_size), 
                ymin = 2.5, 
                xmax = end_of_box,
                ymax=max(y2 + padding_size)), 
            fill = NA, color = "darkgrey", size = 0.25) +
  scale_fill_identity() +
  scale_color_identity() +
  theme_void() + 
  theme(
    panel.grid.major.x = element_line(color="lightgrey", size=0.2),
    ) + 
 scale_x_continuous(
   breaks = seq(min(read_type_schematic_data$x1) - padding_size + 0.25, end_of_box , 0.25))

read_type_schematic

```

# Generate schematic defining read type fractions
```{r rna-seq-survey-MS-v20-27 }
stat_names <- c("pct_unmapped_of_total",  "pct_dupe_of_mapped", "pct_non_exonic_of_non_dupe")
stat_labels <- c("Unmapped",  "Duplicate", "Non-exonic") 
complement_stat_names <- c("Mapped reads",  "Non-duplicate reads", "Exonic reads") 
divisor_name = c("total", "mapped", "non-duplicate")

comparison_of_remaining <- read_type_fractions_long %>% 
  mutate(read_type_fraction_name = complement_stat_names[match(read_type_fraction_name, stat_names)]) %>%
  group_by(read_type_fraction_name) %>%
  summarize(min = round(1-max(dataset_value), 2),
    max = round(1-min(dataset_value), 2),
    median = round(1-median(dataset_value), 2),
    p25 = round(1-quantile(dataset_value, 0.25), 2),
    p75 = round(1-quantile(dataset_value, 0.75), 2)) %>%
#      mutate(bar_label = paste0 (read_type_fraction_name, "\n(", 100*min, "-", 100*max, "% of ", divisor_name[match(read_type_fraction_name, complement_stat_names)], "; x͂=",100* median, "%)"))
mutate(bar_label = paste0 (read_type_fraction_name, "\n(", 100*min, "-", 100*max, "% of ", divisor_name[match(read_type_fraction_name, complement_stat_names)], "; med.=", 100* median, "%)"))

positive_bars <- comparison_of_remaining %>%
  select(-min, -max, -p25, -p75) %>%
  mutate(read_type_fraction_label = read_type_fraction_name, 
         position = "1") %>%
  rename(abs_median = median)

negative_bars <- tibble(read_type_fraction_name = positive_bars$read_type_fraction_name,
                        read_type_fraction_label = stat_labels[match(read_type_fraction_name, complement_stat_names)],
                        bar_label = read_type_fraction_label,
                        abs_median = 1-positive_bars$abs_median,
                        position = "2")

total_bar <- tibble(read_type_fraction_name = "Total reads",
                    read_type_fraction_label = read_type_fraction_name,
                    bar_label = "Total reads\n(100% of reads)",
                    abs_median = 1,
                    position = "1"
                    )


MEND_stats_of_total <- counts_and_meta  %>%
  arrange(desc(Total_reads)) %>%
  mutate(
    pct_MEND_of_total = (Total_reads - MEND) / Total_reads
  ) %>%
  ungroup %>% 
  summarize(min = round(1-max(pct_MEND_of_total), 2),
    max = round(1-min(pct_MEND_of_total), 2),
    median = round(1-median(pct_MEND_of_total), 2))

MEND_bar <- tibble(read_type_fraction_name = "MEND reads",
                   read_type_fraction_label = read_type_fraction_name,
                   bar_label = with(MEND_stats_of_total, paste0 ("MEND reads\n(", 100*min, "-", 100*max, "% of total reads; med.=",100* median, "%)")),
                   abs_median = 1,
                   position = "1"
)

bar_names_in_order <- c("Total reads", "Mapped reads", "Non-duplicate reads", "Exonic reads", "MEND reads")
  
all_bars <- bind_rows(positive_bars, negative_bars, total_bar, MEND_bar) %>%
  mutate(read_type_fraction_name = factor(read_type_fraction_name, 
                                          levels = bar_names_in_order))

cat_space <- all_bars %>% filter(position == 1) %>%
  arrange(read_type_fraction_name) %>%
  select(read_type_fraction_name, abs_median) %>%
  ungroup %>%
  mutate(lagged1_abs_median = lag(abs_median, default =1),
         lagged2_abs_median = lag(abs_median, n=2, default =1),
         category_space = lagged1_abs_median*lagged2_abs_median
         ) %>%
  select(read_type_fraction_name, category_space)

these_colors_with_MEND = c(these_colors, "MEND reads"="#000000")


all_bars_rel <- left_join(all_bars, cat_space, by = "read_type_fraction_name") %>%
  mutate(rel_median=ifelse(read_type_fraction_name == "MEND reads",
                           abs_median[read_type_fraction_label=="Exonic reads"]*
                             category_space[read_type_fraction_label=="Exonic reads"], 
                           abs_median*category_space),
         read_type_fraction_name = factor(read_type_fraction_name, levels = rev(bar_names_in_order)),
         position = factor(position, levels = rev(sort(unique(position)))),
         bar_color = these_colors_with_MEND[match(read_type_fraction_name, names(these_colors_with_MEND))],
         fill_val = ifelse(position == 1, bar_color, NA),
         border_color_val = bar_color,
         text_color = ifelse(position ==1, "white", "black"),
font_size = ifelse(abs_median<0.1, 3,3.5),
text_angle = ifelse(abs_median<0.1, 90,0)
         ) 

```

# read_type_fractions_schematic
```{r rna-seq-survey-MS-v20-28, fig.height=4, fig.width=8}
read_type_fractions_schematic <- ggplot(all_bars_rel, aes(y=read_type_fraction_name,
                x=rel_median, group = position)) + 
   geom_col(aes(fill = fill_val,
                color = border_color_val
                )) +
   geom_text(aes(label=bar_label,
                 color = text_color,
                 size = font_size,
                 angle= text_angle), 
             position = position_stack(vjust = .5),
              fontface = "bold") +
   scale_fill_identity() +
   scale_color_identity()  +
  scale_size_identity() + 
   theme_void() +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank()) 


ggplot(all_bars_rel, aes(y=read_type_fraction_name,
                x=rel_median, group = position)) + 
   geom_col(aes(fill = fill_val,
                color = border_color_val
                )) +
   geom_text(aes(label=read_type_fraction_label,
                 color = text_color,
                 size = font_size,
                 angle = text_angle), 
             position = position_stack(vjust = .5),
              fontface = "bold") +
   scale_fill_identity() +
   scale_color_identity()  +
  scale_size_identity() + 
   theme_void() +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank())
 

read_type_fractions_schematic

```

# Figure 2
```{r rna-seq-survey-MS-v20-29, fig.height=6, fig.width=5}

figure_name <- "f2"

f2 <- plot_grid(read_type_schematic,
read_type_fractions_schematic,
nrow=2,
rel_heights = c(1, 2),
labels = c("A", "B")) +
   draw_label(paste(figure_name, Sys.time()),
              x = 0, y = 0.01, hjust = 0, size = 6) 
# 
# f2 <- plot_grid(left_side,
# right_side,
# ncol = 2) +
#    draw_label(paste(figure_name, Sys.time()),
#               x = 0, y = 0.02, hjust = 0, size = 6) 

f2
ggsave("results/Figure_2.png", f2, height=6, width=5)
ggsave("results/Figure_2.pdf", f2, height=6, width=5)

```


```{r rna-seq-survey-MS-v20-30 }

stat_names <- c("pct_unmapped_of_total",  "pct_dupe_of_mapped", "pct_non_exonic_of_non_dupe")
stat_label <- c("% unmapped", "% duplicate \n(of mapped)", "% non-exonic \n(of non-duplicate)")
names(stat_label) <- stat_names


```

# Plot read fractions
```{r rna-seq-survey-MS-v20-31, fig.height=4.5, fig.width=8.25}

colors_for_stat_names <- these_colors[c("Mapped reads", "Non-duplicate reads","Exonic reads")]
names(colors_for_stat_names) <- stat_names

these_colors <- c("#A6CEE3", "#1F78B4", "#B2DF8A", "#33A02C", "#FB9A99", "#E31A1C", 
"#FDBF6F", "#FF7F00")
# light blue, dark blue, etc: green, red, orange (then purple and brown)
names(these_colors) <-  c("not assigned", "Total reads", "Non-exonic reads",  "Exonic reads", "Duplicate reads",  "Non-duplicate reads",  "Unmapped reads", "Mapped reads")

# 6 x 11 was a good ratio, but the text was too small

read_type_fractions_long_for_histogram <- read_type_fractions_long %>%
  mutate(read_type_fraction_name = factor(read_type_fraction_name,
                                          levels = stat_names))

read_type_fractions_histogram <-  ggplot(read_type_fractions_long_for_histogram) +
  geom_histogram(aes(x = dataset_value*100, color = read_type_fraction_name), 
                 fill = "white", stat = StatBin2)  +
  scale_color_manual(values = colors_for_stat_names) +
  facet_wrap(~read_type_fraction_name, 
             nrow = 1, 
             scales = "free_y",
             labeller = labeller(
               read_type_fraction_name = stat_label
             )
  ) +
  ylab("Datasets") +
  xlab(paste0("Read type percentages in ", length(unique(read_type_fractions_long_for_histogram$dataset_id)),  " datasets")) +
#  scale_x_continuous(breaks = seq(0, 1, by = 0.5)) +
    scale_x_continuous(breaks = seq(0, 100, by = 50)) +
  theme_minimal() +
  theme(axis.line.x = element_line(color="darkgrey", size = 0.5),
        axis.line.y = element_line(color="darkgrey", size = 0.5)) +
  theme(strip.text.x = element_text(size = 12, face = "bold")) +
  theme(legend.position = "none")

read_type_fractions_histogram

```


# Calculate MEND reads as a fraction of total
```{r rna-seq-survey-MS-v20-32 }

read_counts_with_MEND_fractions <- counts_and_meta  %>%
  arrange(desc(Total_reads)) %>%
  mutate(pct_MEND_of_total = MEND /Total_reads)

ggplot(read_counts_with_MEND_fractions) + 
  geom_histogram(aes(pct_MEND_of_total), stat = StatBin2) +
    theme_minimal(base_size = 20) +
    theme(axis.line.x = element_line(color="darkgrey", size = 0.5),
        axis.line.y = element_line(color="darkgrey", size = 0.5)) 

MEND_pct <- read_counts_with_MEND_fractions %>% 
  ungroup %>%
  mutate(dataset_value = pct_MEND_of_total) %>%
  summarize(variance= var(dataset_value),
            min = min(dataset_value),
            max = max(dataset_value),
            median = median(dataset_value),
            p25 = quantile(dataset_value, 0.25),
            p75 = quantile(dataset_value, 0.75),
            iqr = IQR(dataset_value))

kable(MEND_pct %>%
        select(-variance) %>%
      mutate_if(is.double, ~100*.), digits = 0) %>% 
  kable_styling(bootstrap_options = "striped", full_width = F)

head(sort(read_counts_with_MEND_fractions$pct_MEND_of_total))

sum(read_counts_with_MEND_fractions$pct_MEND_of_total<0.25)

read_counts_with_MEND_fractions %>%
  ungroup %>%
  mutate(pct_under_25 = pct_MEND_of_total<0.25) %>%
  tabyl(pct_under_25) %>%
  adorn_totals()

read_counts_with_MEND_fractions %>%
  ungroup %>%
  filter(Mapped > 30E6) %>%
  mutate(pct_under_25 = pct_MEND_of_total<0.25) %>%
  tabyl(pct_under_25) %>%
  adorn_totals()



```


# Show per-cohort distribution of read_type_fractions 
```{r rna-seq-survey-MS-v20-33, fig.height=6, fig.width=8}

figure_name <- "fracs_by_group"

read_type_fractions_long_for_boxplot <- read_type_fractions_long_for_histogram %>%
  ungroup %>%
  mutate(boxplot_label = fct_reorder(paste0(as.character(cohort_code), ", n=", n_in_cohort), n_in_cohort))


read_type_fractions_per_cohort_boxplot <- ggplot(read_type_fractions_long_for_boxplot) +
  geom_boxplot(aes(x = dataset_value*100, y=boxplot_label), outlier.size = 0.5)  +
  facet_wrap(~read_type_fraction_name, 
             nrow = 1, 
             labeller = labeller(
               read_type_fraction_name = stat_label
             )
  ) +
  ylab("Cohorts") +
  xlab(paste0("Read type percentages in ", length(unique(read_type_fractions_long_for_boxplot$dataset_id)),  " datasets")) +
  #  scale_x_continuous(breaks = seq(0, 1, by = 0.5)) +
  scale_x_continuous(breaks = seq(0, 100, by = 50)) +
  
  theme_minimal() +
  theme(axis.line.x = element_line(color="darkgrey", size = 0.5),
        axis.line.y = element_line(color="darkgrey", size = 0.5)) +
  theme(strip.text.x = element_text(size = 12, face = "bold")) +
  theme(legend.position = "none")


read_type_fractions_per_cohort_boxplot

```


# Characterize Cohort 4

```{r rna-seq-survey-MS-v20-34 }

read_type_names <- c("Total_reads", "Mapped",  "MND", "MEND")

```

```{r rna-seq-survey-MS-v20-35 }
this_cohort_name <- "EGAD00001003279"

this_cohort_code <- subset(read_counts_with_read_type_fractions, cohort_name == this_cohort_name)$cohort_code[1]

this_cohort_code
```

## in MS
```{r rna-seq-survey-MS-v20-36 }

 read_counts_with_read_type_fractions %>%
  filter(cohort_code == this_cohort_code) %>%
   mutate(gt_98 = pct_dupe_of_mapped > 0.98) %>%
   tabyl(gt_98) %>%
   # add_totals_row
  adorn_totals("row")

 read_counts_with_read_type_fractions %>%
  filter(cohort_code == this_cohort_code) %>%
  filter(pct_dupe_of_mapped > 0.98) %>%
  pull(Total_reads) %>% min

 read_counts_with_read_type_fractions %>%
  filter(cohort_code == this_cohort_code) %>%
  filter(pct_dupe_of_mapped > 0.98) %>%
   pull(genes_with_at_least_1_count) %>% max
 
```

## other
```{r rna-seq-survey-MS-v20-37, fig.height=6, fig.width=8}


 counts_and_meta_long <- counts_and_meta  %>%
  pivot_longer(cols = c(Total_reads, Mapped, MEND, MND), names_to = "read_type_name", values_to = "dataset_value")  %>%
   ungroup %>%
   mutate(read_type_name = factor(read_type_name, levels = read_type_names))

counts_and_meta_long_one_proj <-  counts_and_meta_long %>%
  filter(cohort_code == this_cohort_code)

table(counts_and_meta_long_one_proj$read_type_name)
 
 
# 78 datasets in the cohort have fewer than 0.2 million MEND reads. 
sum((counts_and_meta_long_one_proj %>% filter(read_type_name=="MEND"))$dataset_value < 200000)
those_datasets <- counts_and_meta_long_one_proj %>% filter(read_type_name=="MEND", dataset_value < 200000) %>% pull(dataset_id)


summary((counts_and_meta_long_one_proj %>% filter(read_type_name=="Total_reads", dataset_id %in% those_datasets))$dataset_value)


max_plots_to_include <- 10
n_datasets_to_plot <- ifelse(n_distinct(counts_and_meta_long_one_proj$dataset_id)>max_plots_to_include, 
                         max_plots_to_include, 
                         n_distinct(counts_and_meta_long_one_proj$dataset_id))
plot_word <- ifelse(n_datasets_to_plot==max_plots_to_include, max_plots_to_include, paste("all", n_distinct(counts_and_meta_long_one_proj$dataset_id)))

set.seed(100)
datasets_to_plot <- sample(unique(counts_and_meta_long_one_proj$dataset_id), n_datasets_to_plot)

ggplot(subset(counts_and_meta_long_one_proj, dataset_id %in% datasets_to_plot)) + 
  geom_col(aes(x=read_type_name, y=dataset_value/1e6, fill = read_type_name)) +
  facet_wrap(~dataset_id, ncol = 1, strip.position="right")  +
    theme_minimal(base_size = 20) +
    theme(axis.line.x = element_line(color="darkgrey", size = 0.5),
        axis.line.y = element_line(color="darkgrey", size = 0.5)) +
  theme(strip.text.y = element_text(angle = 0)) +
  scale_fill_brewer(palette = "Set1")  +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  ggtitle(paste("Read counts for", plot_word, "datasets from", this_cohort_code)) +
  coord_flip()


ggplot(counts_and_meta_long_one_proj %>% filter(read_type_name == "MEND")) + 
  geom_histogram(aes(dataset_value/1e6), stat = StatBin2) +
    theme_minimal(base_size = 20) +
    theme(axis.line.x = element_line(color="darkgrey", size = 0.5),
        axis.line.y = element_line(color="darkgrey", size = 0.5)) +
  ggtitle(paste("Distribution of MEND counts in", this_cohort_code))

bl <- colorRampPalette(c("navy","royalblue","lightskyblue"))(200)                      
re <- colorRampPalette(c("mistyrose", "red2","darkred"))(200)


ggplot(counts_and_meta %>% 
         filter(cohort_code == this_cohort_code) %>%
         mutate(pct_dupes = (Mapped - MND)/Mapped)) + 
  geom_point(aes(x=MEND/1e6, y=pct_dupes, fill = Total_reads/1e6), shape=21, color = "black") +
    theme_minimal(base_size = 20) +
    theme(axis.line.x = element_line(color="darkgrey", size = 0.5),
        axis.line.y = element_line(color="darkgrey", size = 0.5)) +
  scale_fill_gradientn(colours = c(bl,"white", re)) +
  ggtitle(paste("Relationship of dupe fraction to MEND count in", this_cohort_code), "For detecting whether datasets with less info were sequenced more deeply. \nNormally, higher depth datasets have more MEND reads and more duplicate \nreads than the same dataset sequenced at a lower total depth.")

counts_and_meta_long_one_proj %>%
  select(seq_length, dataset_id) %>%
  distinct %>%
  tabyl(seq_length)

```


# 
```{r}
# what dataset outside this project has the highest fraction of duplicate reads

 read_counts_with_read_type_fractions %>%
  filter(cohort_code != this_cohort_code) %>% 
  arrange(desc(pct_dupe_of_mapped))

# samples with greater than 90% duplicates
 read_counts_with_read_type_fractions %>%
  filter(cohort_code != this_cohort_code) %>% 
   filter(pct_dupe_of_mapped > 0.9) %>%
   select(dataset_id, n_in_cohort, Total_reads, pct_dupe_of_mapped, genes_with_at_least_1_count) %>% 
  arrange(desc(pct_dupe_of_mapped))
 
 
  
 
 summary(read_counts_with_read_type_fractions$genes_with_at_least_1_count)
 
 


```


# Expressed genes

## What is the distribution of measured genes?
```{r rna-seq-survey-MS-v20-38 }

gene_counts_long <- counts_and_meta %>%
  pivot_longer(cols = starts_with("genes_"), names_to = "Expression_threshold")  %>%
  mutate(Expression_threshold = paste(str_replace(Expression_threshold, "genes_expressed_above_", "genes > "), "TPM"))  


ggplot(gene_counts_long) + 
  geom_histogram(aes(x=value/1e3, fill = Expression_threshold), stat = StatBin) +
  facet_wrap(~Expression_threshold)  +
    theme_minimal() +
  scale_fill_brewer(palette = "Set1") +
    theme(axis.line.x = element_line(color="darkgrey", size = 0.5),
        axis.line.y = element_line(color="darkgrey", size = 0.5)) +
    theme(legend.position = "none") +
  ylab("Datasets") +
  xlab("Measured genes") 
  
```

# What are the exact lowest gene counts?
```{r rna-seq-survey-MS-v20-39 }

sort(counts_and_meta$genes_with_at_least_1_count) [1:100]

```

```{r}
n_datasets_above_31_Total_reads <- counts_and_meta %>%
  filter(Total_reads>31E6) %>%
  nrow

n_datasets_above_155_MEND_read_and_above_31_Total_reads <- counts_and_meta %>%
  filter(MEND<=15.5e6 & Total_reads>31E6) %>%
  nrow

n_datasets_above_31_Total_reads
n_datasets_above_155_MEND_read_and_above_31_Total_reads


n_datasets_above_155_MEND_read_and_above_31_Total_reads/n_datasets_above_31_Total_reads

# sampels that would be excluded


# samples with fewer than 10% mend reads
n_datasets_above_31_Total_reads_and_fewer_than_10_pct_MEND <- counts_and_meta %>%
  filter(MEND/Total_reads < 0.1 & Total_reads>31E6) %>%
  nrow

n_datasets_above_31_Total_reads
n_datasets_above_31_Total_reads_and_fewer_than_10_pct_MEND


n_datasets_above_31_Total_reads_and_fewer_than_10_pct_MEND/n_datasets_above_31_Total_reads




```

# Checking distribution of values
```{r}

# 
# ggplot(this_data,
#        aes(x=Read_count/1E6,fill = Read_type)) + 
#   geom_histogram(alpha = 0.5, shape = 20, size =0.1) +
# #  geom_smooth(method = "lm", color = "black") +
# # geom_label(data = these_stats, 
#              # aes(label=paste0("r=", r_corr, "\np=", format(r_p_value, digits = 1))), 
#              # x=Inf, #max(this_data$Read_count/1E6),
#              # y=Inf, #max(this_data$Gene_count/1e3),
#              # fill = NA,
#              # hjust = 1, vjust = 1, size = 3
#              # ) +
#     theme_minimal() +
#     theme(axis.line.x = element_line(color="darkgrey", size = 0.5),
#         axis.line.y = element_line(color="darkgrey", size = 0.5),
#         strip.text = element_text(face = "bold")) +
#   #scale_fill_manual(values = colors_for_correlation_plot) +
#   facet_grid(Read_type~Expression_threshold) +
#   theme(legend.position="none")
# 
# 
# 
# ggplot(this_data,
#        aes(x=Gene_count/1e3,fill = Read_type)) + 
#   geom_histogram(alpha = 0.5, shape = 20, size =0.1) +
#     theme_minimal() +
#     theme(axis.line.x = element_line(color="darkgrey", size = 0.5),
#         axis.line.y = element_line(color="darkgrey", size = 0.5),
#         strip.text = element_text(face = "bold")) +
#   scale_fill_manual(values = colors_for_correlation_plot) +
#   facet_grid(Read_type~Expression_threshold) +
#   theme(legend.position="none") 
# 
# this_data %>%
#   group_by(Read_type, Expression_threshold) %>%
#   summarize(Gene_count_st = list(shapiro.test(Gene_count)),
#          Read_count_st = list(shapiro.test(Read_count))) %>%
#   select(Gene_count_st, Read_count_st)

```

# Pipeline duration
```{r}

# The 382 RNA-Seq datasets recorded were obtained from SRA and EGA and had reads that were 100 bases in length

sum(! is.na(counts_and_meta$duration_of_expression))
sum(! is.na(counts_and_meta$duration_of_mend))

# For datasets within 10 million total reads of the median total read size in our survey of 2179 datasets (61 million total reads), the median duration is 290 minutes (4.8 hours) for the expression pipeline and 143 minutes (2.4 hours) for the MEND pipeline. 

median_total_reads <- median(counts_and_meta$Total_reads)

total_read_range_expansion <- 10E6

typical_size_samples <- counts_and_meta %>%
  ungroup %>%
  filter(Total_reads > median_total_reads - total_read_range_expansion, 
         Total_reads < median_total_reads + total_read_range_expansion)

medians <- typical_size_samples %>%
summarize(median_expr_duration = round(median(duration_of_expression, na.rm = TRUE)),
          median_mend_duration = round(median(duration_of_mend, na.rm = TRUE)))

medians
round(medians/60, 1)

```


# Figure 3
```{r rna-seq-survey-MS-v20-48, fig.height=8, fig.width=8}
# 
f3a <- read_type_fractions_histogram +
            theme(strip.text.x = element_text(size = 10, face = "bold"),
                  axis.title.x=element_blank(),
                  axis.text.x=element_blank())

f3b <- read_type_fractions_per_cohort_boxplot +
            theme(strip.text.x = element_blank(),
                  strip.background = element_blank(),
                  strip.text = element_blank())

f3ab <- plot_grid(f3a, 
                  f3b,
                  ncol = 1,
                  rel_heights = c(1.5, 6),
                  align = "v",
                  labels = "AUTO")
# f3c <- plot_grid(cor_plots_excluding_high_read_depth_and_low_gene_count_datasets,
#           labels = "C")


# f3 <- plot_grid(f3ab,
#                  f3c,
#                  nrow=1,
#                  rel_widths = c(4,2.3))
# 
figure_name <- "f3"
f3_labeled <- f3ab +   draw_label(paste(figure_name, Sys.time()),
             x = 0, y = 0.02, hjust = 0, size = 6)


f3_labeled
ggsave("results/Figure_3.png", f3_labeled, height=8, width=10)
ggsave("results/Figure_3.pdf", f3_labeled, height=8, width=10)


```


# Plot sequence length against pct duplicates
```{r rna-seq-survey-MS-v20-49 }

ggplot(read_counts_with_read_type_fractions) + 
  geom_point(aes(x=seq_length, y=pct_dupe_of_mapped))


ggplot(read_counts_with_read_type_fractions) + 
  geom_boxplot(aes(x=seq_length, y=pct_dupe_of_mapped, group = seq_length))

```

# Generate bibtex formatted citations for packages
```{r rna-seq-survey-MS-v20-50}

used_pkgs <- c("base", "tidyverse", "janitor", "knitr", "cowplot", "RColorBrewer", "pander", "kableExtra", "snakecase", "corrr")
write_bib(used_pkgs)
```

## Does the depth of sequencing correlate to the number of duplicates? 
```{r rna-seq-survey-MS-v20-51 }

v1 <- read_counts_with_read_type_fractions$Total_reads
v2 <- read_counts_with_read_type_fractions$pct_dupe_of_mapped
# 
# # cor_tot_reads_and_dupes <- cor(v1, 
# #     v2,
# #     method = c("pearson"))
# # 
# # round(cor_tot_reads_and_dupes, 2)
# # round(cor_tot_reads_and_dupes^2, 2)
# # 
# # c1 <- cor(v1, 
# #     v2,
# #     method = c("kendall"))
# # c1
# # c1^2
# 
# c2 <- cor(v1, 
#     v2,
#     method = "spearman")
# c2
# c2^2
# 
# # 
# # cor.test(v1, v2, 
# #     alternative = c("two.sided"),
# #     method = c("pearson"))
# # 
# # cor.test(v1, v2, 
# #     alternative = c("two.sided"),
# #     method = c("kendall"))
# 
# cor.test(v1, v2, 
#     alternative = c("two.sided"),
#     method = "spearman")
# 

cor.test(v1, v2, 
    alternative = "greater", # "greater" corresponds to positive association, 
    method = "spearman")
t.test(v1, v2, 
       alternative = "greater",
       var.equal=FALSE)$p.value

# 
# # test normality
# shapiro.test(v1)
# shapiro.test(v2)
# 
# ggplot(read_counts_with_read_type_fractions) +
#   geom_histogram(aes(Total_reads)) +
#   theme_minimal()
# 
# ggplot(read_counts_with_read_type_fractions) +
#   geom_histogram(aes(pct_dupe_of_mapped))  +
#   theme_minimal()
# 
# ggpubr::ggqqplot(v1)
# ggpubr::ggqqplot(v2)
```

# SessionInfo
```{r rna-seq-survey-MS-v20-52 }
pander(sessionInfo(), compact = FALSE)
```