TMRC3 Comprehensive Data Analysis: 202104

1 Introduction

This document is intended to provide an overview of TMRC3 samples which have been sequenced. It includes some plots and analyses showing the relationships among the samples as well as some differential analyses when possible.

2 Annotation

We take the annotation data from ensembl’s biomart instance. The genome which was used to map the data was hg38 revision 100. My default when using biomart is to load the data from 1 year before the current date.

hs_annot <- sm(load_biomart_annotations(year = "2020"))
hs_annot <- hs_annot[["annotation"]]
hs_annot[["transcript"]] <- paste0(rownames(hs_annot), ".", hs_annot[["version"]])
rownames(hs_annot) <- make.names(hs_annot[["ensembl_gene_id"]], unique = TRUE)
tx_gene_map <- hs_annot[, c("transcript", "ensembl_gene_id")]

summary(hs_annot)

##  ensembl_transcript_id ensembl_gene_id       version     transcript_version
##  Length:227921         Length:227921      Min.   : 1.0   Min.   : 1.00     
##  Class :character      Class :character   1st Qu.: 6.0   1st Qu.: 1.00     
##  Mode  :character      Mode  :character   Median :12.0   Median : 1.00     
##                                           Mean   :10.7   Mean   : 3.08     
##                                           3rd Qu.:16.0   3rd Qu.: 5.00     
##                                           Max.   :29.0   Max.   :17.00     
##                                                                            
##  hgnc_symbol        description        gene_biotype         cds_length    
##  Length:227921      Length:227921      Length:227921      Min.   :     3  
##  Class :character   Class :character   Class :character   1st Qu.:   357  
##  Mode  :character   Mode  :character   Mode  :character   Median :   694  
##                                                           Mean   :  1139  
##                                                           3rd Qu.:  1446  
##                                                           Max.   :107976  
##                                                           NA's   :127343  
##  chromosome_name       strand          start_position      end_position     
##  Length:227921      Length:227921      Min.   :5.77e+02   Min.   :6.47e+02  
##  Class :character   Class :character   1st Qu.:3.11e+07   1st Qu.:3.12e+07  
##  Mode  :character   Mode  :character   Median :6.04e+07   Median :6.06e+07  
##                                        Mean   :7.41e+07   Mean   :7.42e+07  
##                                        3rd Qu.:1.09e+08   3rd Qu.:1.09e+08  
##                                        Max.   :2.49e+08   Max.   :2.49e+08  
##                                                                             
##   transcript       
##  Length:227921     
##  Class :character  
##  Mode  :character  
##                    
##                    
##                    
## 

hs_go <- sm(load_biomart_go()[["go"]])
hs_length <- hs_annot[, c("ensembl_gene_id", "cds_length")]
colnames(hs_length) <- c("ID", "length")

3 Sample Estimation

3.1 Generate expressionsets

The sample sheet is copied from our shared online sheet and updated with each release of sequencing data.

samplesheet <- "sample_sheets/tmrc3_samples_20210512.xlsx"

3.1.1 Hisat2 expressionsets

The first thing to note is the large range in coverage. There are multiple samples with coverage which is too low to use. These will be removed shortly.

In the following block I immediately exclude any non-coding reads as well.

## Create the expressionset and immediately pass it to a filter
## removing the non protein coding genes.
hs_expt <- create_expt(samplesheet,
                       file_column = "hg38100hisatfile",
                       savefile = glue::glue("rda/hs_expt_all-v{ver}.rda"),
                       gene_info = hs_annot) %>%
  exclude_genes_expt(column = "gene_biotype", method = "keep", pattern = "protein_coding")

## Reading the sample metadata.

## Dropped 113 rows from the sample metadata because they were blank.

## The sample definitions comprises: 131 rows(samples) and 79 columns(metadata fields).

## Reading count tables.

## Reading count files with read.table().

## /fs01/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/TMRC30001/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows.

## preprocessing/TMRC30002/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30003/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30004/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30005/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30006/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30007/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30008/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30009/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30010/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30015/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30011/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30012/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30013/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30016/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30017/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30050/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30052/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30071/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30056/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30058/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30105/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30094/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30080/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30103/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30018/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30107/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30083/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30019/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30082/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30093/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30113/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30096/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30014/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30021/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30029/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30020/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30038/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30039/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30023/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30025/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30022/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30046/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30047/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30048/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30026/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30030/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30031/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30032/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30024/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30040/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30033/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30049/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30053/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30054/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30115/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30037/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30027/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30028/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30034/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30035/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30036/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30044/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30055/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30068/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30070/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30041/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30042/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30043/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30045/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30059/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30060/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30061/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30062/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30063/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30051/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30064/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30065/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30066/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30067/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30057/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30069/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30116/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30074/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30072/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30076/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30077/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30078/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30088/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30079/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30097/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30075/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30085/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30086/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30087/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30101/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30089/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30090/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30081/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30092/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30104/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30106/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30114/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30095/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30108/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30130/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30124/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30131/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30109/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30084/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30098/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30099/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30100/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30110/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30111/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30102/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30091/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30112/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## preprocessing/TMRC30073/outputs/hisat2_hg38_100/r1_trimmed.count_hg38_100_sno_gene_gene_id.count.xz contains 21486 rows and merges to 21486 rows.

## Finished reading count data.

## Warning in create_expt(samplesheet, file_column = "hg38100hisatfile", savefile =
## glue::glue("rda/hs_expt_all-v{ver}.rda"), : Some samples were removed when cross
## referencing the samples against the count data.

## Matched 21452 annotations and counts.

## Bringing together the count matrix and gene information.

## Some annotations were lost in merging, setting them to 'undefined'.

## The final expressionset has 21481 rows and 119 columns.

## Before removal, there were 21481 genes, now there are 19941.

## There are 13 samples which kept less than 90 percent counts.

## TMRC30015 TMRC30017 TMRC30019 TMRC30044 TMRC30045 TMRC30097 TMRC30075 TMRC30087 
##     79.24     85.72     89.75     80.34     73.33     89.90     86.97     83.63 
## TMRC30101 TMRC30104 TMRC30114 TMRC30131 TMRC30073 
##     88.41     80.29     87.62     86.82     89.26

3.1.1.1 Initial metrics

Once the data was loaded, there are a couple of metrics which may be plotted immediately.

nonzero <- plot_nonzero(hs_expt)
nonzero$plot

## Warning: ggrepel: 89 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

Najib doesn’t want this plot, but I am using it to check new samples, so will hide it from general use.

libsize <- plot_libsize(hs_expt)
libsize$plot

3.2 Minimum coverage sample filtering

I arbitrarily chose 11,000 non-zero genes as a minimum. We may want this to be higher.

hs_valid <- subset_expt(hs_expt, nonzero = 11000)

## The samples (and read coverage) removed when filtering 11000 non-zero genes are:

## TMRC30010 TMRC30050 TMRC30052 
##     52471    808149   3087347

## There were 119, now there are 116 samples.

valid_write <- sm(write_expt(hs_valid, excel = glue("excel/hs_valid-v{ver}.xlsx")))

4 Project Aims

The project seeks to determine the relationship of the innate immune response and inflammatory signaling to the clinical outcome of antileishmanial drug treatment. We will test the hypothesis that the profile of innate immune cell activation and their dynamics through the course of treatment differ between CL patients with prospectively determined therapeutic cure or failure.

This will be achieved through the characterization of the in vivo dynamics of blood-derived monocyte, neutrophil and eosinophil transcriptome before, during and at the end of treatment in CL patients. Cell-type specific transcriptomes, composite signatures and time-response expression profiles will be contrasted among patients with therapeutic cure or failure.

4.1 Preparation

To address these, I added to the end of the sample sheet columns named ‘condition’, ‘batch’, ‘donor’, and ‘time’. These are filled in with shorthand values according to the above.

4.2 Global view

Before addressing the questions explicitly by subsetting the data, I want to get a look at the samples as they are.

hs_valid <- hs_valid %>%
  set_expt_batches(fact = "cellssource") %>%
  set_expt_conditions(fact = "typeofcells") %>%
  set_expt_samplenames(newnames = pData(hs_valid)[["samplename"]])

all_norm <- sm(normalize_expt(hs_valid, transform = "log2", norm = "quant",
                              convert = "cpm", filter = TRUE))

all_pca <- plot_pca(all_norm, plot_labels = FALSE, plot_title = "PCA - Cell type")
pp(file = glue("images/tmrc3_pca_nolabels-v{ver}.png"), image = all_pca$plot)

write.csv(all_pca$table, file = "coords/hs_donor_pca_coords.csv")
plot_corheat(all_norm, plot_title = "Heirarchical clustering - cell types")$plot

4.3 Examine samples relevant to clinical outcome

Now let us consider only the samples for which we have a clinical outcome. These fall primarily into either ‘cured’ or ‘failed’, but some people have not yet returned to the clinic after the first or second visit. These are deemed ‘lost’.

chosen_colors <- c("#D95F02", "#7570B3", "#1B9E77", "#FF0000")
names(chosen_colors) <- c("cure", "failure", "lost", "null")
newnames <- pData(hs_valid)[["samplename"]]

hs_clinical <- hs_valid %>%
  set_expt_conditions(fact = "clinicaloutcome") %>%
  set_expt_batches(fact = "typeofcells") %>%
  set_expt_colors(colors = chosen_colors) %>%
  set_expt_samplenames(newnames = newnames) %>%
  subset_expt(subset = "typeofcells!='PBMCs'&typeofcells!='Macrophages'")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'sampleNames': 'names' attribute [116] must be the same length as the vector [87]

hs_clinical_norm <- sm(normalize_expt(hs_clinical, filter = TRUE, transform = "log2",
                                      convert = "cpm", norm = "quant"))

## Error in normalize_expt(hs_clinical, filter = TRUE, transform = "log2", : object 'hs_clinical' not found

clinical_pca <- plot_pca(hs_clinical_norm, plot_labels = FALSE, cis = NULL,
                         plot_title = "PCA - clinical samples")

## Error in plot_pca(hs_clinical_norm, plot_labels = FALSE, cis = NULL, plot_title = "PCA - clinical samples"): object 'hs_clinical_norm' not found

pp(file = glue("images/all_clinical_nobatch_pca-v{ver}.png"), image = clinical_pca$plot,
   height = 8, width = 20)

## Error in pp(file = glue("images/all_clinical_nobatch_pca-v{ver}.png"), : object 'clinical_pca' not found

4.3.1 Repeat without the biopsy samples

hs_clinical_nobiop <- hs_clinical %>%
  subset_expt(subset = "typeofcells!='Biopsy'")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'sampleNames': object 'hs_clinical' not found

hs_clinical_nobiop_norm <- sm(normalize_expt(hs_clinical_nobiop, filter = TRUE, transform = "log2",
                                             convert = "cpm", norm = "quant"))

## Error in normalize_expt(hs_clinical_nobiop, filter = TRUE, transform = "log2", : object 'hs_clinical_nobiop' not found

clinical_nobiop_pca <- plot_pca(hs_clinical_nobiop_norm, plot_labels = FALSE, cis = NULL,
                                plot_title = "PCA - clinical samples without biopsies")

## Error in plot_pca(hs_clinical_nobiop_norm, plot_labels = FALSE, cis = NULL, : object 'hs_clinical_nobiop_norm' not found

pp(file = glue("images/all_clinical_nobiop_nobatch_pca-v{ver}.png"),
   image = clinical_nobiop_pca$plot)

## Error in pp(file = glue("images/all_clinical_nobiop_nobatch_pca-v{ver}.png"), : object 'clinical_nobiop_pca' not found

4.3.2 Attempt to correct for the surrogate variables

At this time we have two primary data structures of interest: hs_clinical and hs_clinical_nobiop

hs_clinical_nb <- normalize_expt(hs_clinical, filter = TRUE, batch = "svaseq",
                                 transform = "log2", convert = "cpm")

## Error in normalize_expt(hs_clinical, filter = TRUE, batch = "svaseq", : object 'hs_clinical' not found

clinical_batch_pca <- plot_pca(hs_clinical_nb, plot_labels = FALSE, cis = NULL,
                               plot_title = "PCA - clinical samples")

## Error in plot_pca(hs_clinical_nb, plot_labels = FALSE, cis = NULL, plot_title = "PCA - clinical samples"): object 'hs_clinical_nb' not found

clinical_batch_pca$plot

## Error in eval(expr, envir, enclos): object 'clinical_batch_pca' not found

hs_clinical_nobiop_nb <- sm(normalize_expt(hs_clinical_nobiop, filter = TRUE, batch = "svaseq",
                                           transform = "log2", convert = "cpm"))

## Error in normalize_expt(hs_clinical_nobiop, filter = TRUE, batch = "svaseq", : object 'hs_clinical_nobiop' not found

clinical_nobiop_batch_pca <- plot_pca(hs_clinical_nobiop_nb,
                                      plot_title = "PCA - clinical samples without biopsies",
                                      plot_labels = FALSE)

## Error in plot_pca(hs_clinical_nobiop_nb, plot_title = "PCA - clinical samples without biopsies", : object 'hs_clinical_nobiop_nb' not found

pp(file = "images/clinical_batch.png", image = clinical_nobiop_batch_pca$plot)

## Error in pp(file = "images/clinical_batch.png", image = clinical_nobiop_batch_pca$plot): object 'clinical_nobiop_batch_pca' not found

clinical_nobiop_batch_tsne <- plot_tsne(hs_clinical_nobiop_nb,
                                        plot_title = "tSNE - clinical samples without biopsies",
                                        plot_labels = FALSE)

## Error in plot_pca(..., pc_method = "tsne"): object 'hs_clinical_nobiop_nb' not found

clinical_nobiop_batch_tsne$plot

## Error in eval(expr, envir, enclos): object 'clinical_nobiop_batch_tsne' not found

test <- simple_varpart(hs_clinical_nobiop)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'hs_clinical_nobiop' not found

test$partition_plot

## Error in eval(expr, envir, enclos): object 'test' not found

4.4 Perform DE of the clinical samples cure vs. fail

individual_celltypes <- subset_expt(hs_clinical_nobiop, subset="condition!='lost'")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'sampleNames': object 'hs_clinical_nobiop' not found

hs_clinic_de <- sm(all_pairwise(individual_celltypes, model_batch = "svaseq", filter = TRUE))

## Error in normalize_expt(input, filter = filter): object 'individual_celltypes' not found

hs_clinic_table <- sm(combine_de_tables(
    hs_clinic_de,
    excel = glue::glue("excel/individual_celltypes_table-v{ver}.xlsx")))

## Error in combine_de_tables(hs_clinic_de, excel = glue::glue("excel/individual_celltypes_table-v{ver}.xlsx")): object 'hs_clinic_de' not found

hs_clinic_sig <- sm(extract_significant_genes(
    hs_clinic_table,
    excel = glue::glue("excel/individual_celltypes_sig-v{ver}.xlsx")))

## Error in extract_significant_genes(hs_clinic_table, excel = glue::glue("excel/individual_celltypes_sig-v{ver}.xlsx")): object 'hs_clinic_table' not found

hs_clinic_sig[["summary_df"]]

## Error in eval(expr, envir, enclos): object 'hs_clinic_sig' not found

4.4.1 Look at only the differential genes

A good suggestion from Theresa was to examine only the most variant genes from failure vs. cure and see how they change the clustering/etc results. This is my attempt to address this query.

hs_clinic_topn <- sm(extract_significant_genes(hs_clinic_table, n = 100))

## Error in extract_significant_genes(hs_clinic_table, n = 100): object 'hs_clinic_table' not found

table <- "failure_vs_cure"
wanted <- rbind(hs_clinic_topn[["deseq"]][["ups"]][[table]],
                hs_clinic_topn[["deseq"]][["downs"]][[table]])

## Error in eval(quote(list(...)), env): object 'hs_clinic_topn' not found

small_expt <- exclude_genes_expt(hs_clinical_nobiop, ids = rownames(wanted), method = "keep")

## Error in exclude_genes_expt(hs_clinical_nobiop, ids = rownames(wanted), : object 'hs_clinical_nobiop' not found

small_norm <- sm(normalize_expt(small_expt, transform = "log2", convert = "cpm",
                                norm = "quant", filter = TRUE))

## Error in normalize_expt(small_expt, transform = "log2", convert = "cpm", : object 'small_expt' not found

plot_pca(small_norm)$plot

## Error in plot_pca(small_norm): object 'small_norm' not found

small_nb <- normalize_expt(small_expt, transform = "log2", convert = "cpm",
                           batch = "svaseq", norm = "quant", filter = TRUE)

## Error in normalize_expt(small_expt, transform = "log2", convert = "cpm", : object 'small_expt' not found

plot_pca(small_nb)$plot

## Error in plot_pca(small_nb): object 'small_nb' not found

## DESeq2 MA plot of failure / cure
hs_clinic_table[["plots"]][["failure_vs_cure"]][["deseq_ma_plots"]]$plot

## Error in eval(expr, envir, enclos): object 'hs_clinic_table' not found

## DESeq2 Volcano plot of failure / cure
hs_clinic_table[["plots"]][["failure_vs_cure"]][["deseq_vol_plots"]]$plot

## Error in eval(expr, envir, enclos): object 'hs_clinic_table' not found

4.4.2 g:Profiler results using the significant up and down genes

ups <- hs_clinic_sig[["deseq"]][["ups"]][[1]]

## Error in eval(expr, envir, enclos): object 'hs_clinic_sig' not found

downs <- hs_clinic_sig[["deseq"]][["downs"]][[1]]

## Error in eval(expr, envir, enclos): object 'hs_clinic_sig' not found

hs_clinic_gprofiler_ups <- simple_gprofiler(ups)

## Error in simple_gprofiler(ups): object 'ups' not found

hs_clinic_gprofiler_ups[["pvalue_plots"]][["bpp_plot_over"]]

## Error in eval(expr, envir, enclos): object 'hs_clinic_gprofiler_ups' not found

hs_clinic_gprofiler_ups[["pvalue_plots"]][["mfp_plot_over"]]

## Error in eval(expr, envir, enclos): object 'hs_clinic_gprofiler_ups' not found

hs_clinic_gprofiler_ups[["pvalue_plots"]][["reactome_plot_over"]]

## Error in eval(expr, envir, enclos): object 'hs_clinic_gprofiler_ups' not found

##hs_try2 <- simple_gprofiler2(ups)

hs_clinic_gprofiler_downs <- simple_gprofiler(downs)

## Error in simple_gprofiler(downs): object 'downs' not found

hs_clinic_gprofiler_downs[["pvalue_plots"]][["bpp_plot_over"]]

## Error in eval(expr, envir, enclos): object 'hs_clinic_gprofiler_downs' not found

hs_clinic_gprofiler_downs[["pvalue_plots"]][["mfp_plot_over"]]

## Error in eval(expr, envir, enclos): object 'hs_clinic_gprofiler_downs' not found

hs_clinic_gprofiler_downs[["pvalue_plots"]][["reactome_plot_over"]]

## Error in eval(expr, envir, enclos): object 'hs_clinic_gprofiler_downs' not found

4.5 Perform GSVA on the clinical samples

hs_celltype_gsva_c2 <- sm(simple_gsva(individual_celltypes))

## Error in simple_gsva(individual_celltypes): object 'individual_celltypes' not found

hs_celltype_gsva_c2_sig <- sm(get_sig_gsva_categories(
    hs_celltype_gsva_c2,
    excel = "excel/individual_celltypes_gsva_c2.xlsx"))

## Error in get_sig_gsva_categories(hs_celltype_gsva_c2, excel = "excel/individual_celltypes_gsva_c2.xlsx"): object 'hs_celltype_gsva_c2' not found

broad_c7 <- GSEABase::getGmt("reference/msigdb/c7.all.v7.2.entrez.gmt",
                             collectionType = GSEABase::BroadCollection(category = "c7"),
                             geneIdType = GSEABase::EntrezIdentifier())
hs_celltype_gsva_c7 <- sm(simple_gsva(individual_celltypes, signatures = broad_c7,
                                      msig_xml = "reference/msigdb_v7.2.xml", cores = 10))

## Error in simple_gsva(individual_celltypes, signatures = broad_c7, msig_xml = "reference/msigdb_v7.2.xml", : object 'individual_celltypes' not found

hs_celltype_gsva_c7_sig <- sm(get_sig_gsva_categories(
    hs_celltype_gsva_c7,
    excel = "excel/individual_celltypes_gsva_c7.xlsx"))

## Error in get_sig_gsva_categories(hs_celltype_gsva_c7, excel = "excel/individual_celltypes_gsva_c7.xlsx"): object 'hs_celltype_gsva_c7' not found

4.5.1 Print some plots of the GSVA outputs

## The raw heatmap of the C2 values
hs_celltype_gsva_c2_sig[["raw_plot"]]

## Error in eval(expr, envir, enclos): object 'hs_celltype_gsva_c2_sig' not found

## The 'significance' scores of the C2 values
hs_celltype_gsva_c2_sig[["score_plot"]]

## Error in eval(expr, envir, enclos): object 'hs_celltype_gsva_c2_sig' not found

## The subset of scores for categories deemed significantly different.
hs_celltype_gsva_c2_sig[["subset_plot"]]

## Error in eval(expr, envir, enclos): object 'hs_celltype_gsva_c2_sig' not found

## The raw heatmap of the C7 values
hs_celltype_gsva_c7_sig[["raw_plot"]]

## Error in eval(expr, envir, enclos): object 'hs_celltype_gsva_c7_sig' not found

## The 'significance' scores of the C7 values
hs_celltype_gsva_c7_sig[["score_plot"]]

## Error in eval(expr, envir, enclos): object 'hs_celltype_gsva_c7_sig' not found

## The subset of scores for categories deemed significantly different.
hs_celltype_gsva_c7_sig[["subset_plot"]]

## Error in eval(expr, envir, enclos): object 'hs_celltype_gsva_c7_sig' not found

5 Individual Cell types

The following blocks split the samples into a few groups by sample type and look at the distributions between them.

5.1 Implementation details

Get top/bottom n genes for each cell type, using clinical outcome as the factor of interest. For the moment, use sva for the DE analysis. Provide cpms for the top/bottom n genes.

Start with top/bottom 200. Perform default logFC and p-value as well.

5.1.1 Shared contrasts

Here is the contrast we will use throughput, I am leaving open the option to add more.

keepers <- list(
  "fail_vs_cure" = c("failure", "cure"))

5.2 Monocytes

5.2.1 Evaluate Monocyte samples

mono <- subset_expt(hs_valid, subset = "typeofcells=='Monocytes'") %>%
  set_expt_conditions(fact = "clinicaloutcome") %>%
  set_expt_batches(fact = "donor") %>%
  set_expt_colors(colors = chosen_colors)

## There were 116, now there are 22 samples.

## FIXME set_expt_colors should speak up if there are mismatches here!!!

save_result <- save(mono, file = "rda/monocyte_expt.rda")
mono_norm <- normalize_expt(mono, convert = "cpm", filter = TRUE,
                            transform = "log2", norm = "quant")

## Removing 8966 low-count genes (10975 remaining).

## transform_counts: Found 9 values equal to 0, adding 1 to the matrix.

plt <- plot_pca(mono_norm, plot_labels = FALSE)$plot
pp(file = glue("images/mono_pca_normalized-v{ver}.pdf"), image = plt)

mono_nb <- normalize_expt(mono, convert = "cpm", filter = TRUE,
                          transform = "log2", batch = "svaseq")

## Removing 8966 low-count genes (10975 remaining).

## batch_counts: Before batch/surrogate estimation, 1202 entries are x==0: 0%.

## batch_counts: Before batch/surrogate estimation, 13643 entries are 0<x<1: 6%.

## Setting 382 low elements to zero.

## transform_counts: Found 382 values equal to 0, adding 1 to the matrix.

plt <- plot_pca(mono_nb, plot_labels = FALSE)$plot
pp(file = glue("images/mono_pca_normalized_batch-v{ver}.pdf"), image = plt)

5.2.2 DE of Monocyte samples

5.2.2.1 Without sva

mono_de <- sm(all_pairwise(mono, model_batch = FALSE, filter = TRUE))

mono_tables <- sm(combine_de_tables(
    mono_de, keepers = keepers,
    excel = glue::glue("excel/monocyte_clinical_all_tables-v{ver}.xlsx")))
written <- write_xlsx(data = mono_tables[["data"]][[1]],
                      excel = glue::glue("excel/monocyte_clinical_table-v{ver}.xlsx"))
mono_sig <- sm(extract_significant_genes(mono_tables, according_to = "deseq"))
written <- write_xlsx(data = mono_sig[["deseq"]][["ups"]][[1]],
                      excel = glue::glue("excel/monocyte_clinical_sigup-v{ver}.xlsx"))
written <- write_xlsx(data = mono_sig[["deseq"]][["downs"]][[1]],
                      excel = glue::glue("excel/monocyte_clinical_sigdown-v{ver}.xlsx"))

mono_pct_sig <- sm(extract_significant_genes(mono_tables, n = 200,
                                             lfc = NULL, p = NULL, according_to = "deseq"))
written <- write_xlsx(data = mono_pct_sig[["deseq"]][["ups"]][[1]],
                      excel = glue::glue("excel/monocyte_clinical_sigup_pct-v{ver}.xlsx"))
written <- write_xlsx(data = mono_pct_sig[["deseq"]][["downs"]][[1]],
                      excel = glue::glue("excel/monocyte_clinical_sigdown_pct-v{ver}.xlsx"))
mono_sig$summary_df

## data frame with 0 columns and 1 row

## Print out a table of the cpm values for other explorations.
mono_cpm <- sm(normalize_expt(mono, convert = "cpm"))
written <- write_xlsx(data = exprs(mono_cpm),
                      excel = glue::glue("excel/monocyte_cpm_before_batch-v{ver}.xlsx"))
mono_bcpm <- sm(normalize_expt(mono, filter = TRUE, convert = "cpm", batch = "svaseq"))
written <- write_xlsx(data = exprs(mono_bcpm),
                      excel = glue::glue("excel/monocyte_cpm_after_batch-v{ver}.xlsx"))

5.2.2.2 With sva

mono_de_sva <- sm(all_pairwise(mono, model_batch = "svaseq", filter = TRUE))

mono_tables_sva <- sm(combine_de_tables(
    mono_de_sva, keepers = keepers,
    excel = glue::glue("excel/monocyte_clinical_all_tables_sva-v{ver}.xlsx")))
mono_sig_sva <- sm(extract_significant_genes(
    mono_tables_sva,
    excel = glue::glue("excel/monocyte_clinical_sig_tables_sva-v{ver}.xlsx"),
    according_to = "deseq"))

5.2.2.3 Monocyte DE plots

First print out the DE plots without and then with sva estimates.

## DESeq2 MA plot of failure / cure
mono_tables[["plots"]][["fail_vs_cure"]][["deseq_ma_plots"]]$plot

## DESeq2 Volcano plot of failure / cure
mono_tables[["plots"]][["fail_vs_cure"]][["deseq_vol_plots"]]$plot

## DESeq2 MA plot of failure / cure with svaseq
mono_tables_sva[["plots"]][["fail_vs_cure"]][["deseq_ma_plots"]]$plot

## DESeq2 Volcano plot of failure / cure with svaseq
mono_tables_sva[["plots"]][["fail_vs_cure"]][["deseq_vol_plots"]]$plot

5.2.2.4 Monocyte ontology search

ups <- mono_sig[["deseq"]][["ups"]][["fail_vs_cure"]]
downs <- mono_sig[["deseq"]][["downs"]][["fail_vs_cure"]]

mono_up_gp <- simple_gprofiler(sig_genes = ups)

## Performing gProfiler GO search of 150 genes against hsapiens.

## GO search found 37 hits.

## Performing gProfiler KEGG search of 150 genes against hsapiens.

## KEGG search found 3 hits.

## Performing gProfiler REAC search of 150 genes against hsapiens.

## REAC search found 0 hits.

## Performing gProfiler MI search of 150 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 150 genes against hsapiens.

## TF search found 5 hits.

## Performing gProfiler CORUM search of 150 genes against hsapiens.

## CORUM search found 1 hits.

## Performing gProfiler HP search of 150 genes against hsapiens.

## HP search found 0 hits.

mono_up_gp[["pvalue_plots"]][["bpp_plot_over"]]

mono_up_gp[["pvalue_plots"]][["mfp_plot_over"]]

mono_up_gp[["pvalue_plots"]][["reactome_plot_over"]]

## NULL

mono_down_gp <- simple_gprofiler(sig_genes = downs)

## Performing gProfiler GO search of 296 genes against hsapiens.

## GO search found 57 hits.

## Performing gProfiler KEGG search of 296 genes against hsapiens.

## KEGG search found 2 hits.

## Performing gProfiler REAC search of 296 genes against hsapiens.

## REAC search found 2 hits.

## Performing gProfiler MI search of 296 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 296 genes against hsapiens.

## TF search found 29 hits.

## Performing gProfiler CORUM search of 296 genes against hsapiens.

## CORUM search found 1 hits.

## Performing gProfiler HP search of 296 genes against hsapiens.

## HP search found 0 hits.

mono_down_gp[["pvalue_plots"]][["bpp_plot_over"]]

mono_down_gp[["pvalue_plots"]][["mfp_plot_over"]]

mono_down_gp[["pvalue_plots"]][["reactome_plot_over"]]

5.2.2.5 Monocyte MSigDB query

broad_c7 <- GSEABase::getGmt("reference/msigdb/c7.all.v7.2.entrez.gmt",
                             collectionType = GSEABase::BroadCollection(category = "c7"),
                             geneIdType = GSEABase::EntrezIdentifier())

mono_up_goseq_msig <- goseq_msigdb(sig_genes = ups, signatures = broad_c7,
                                   signature_category = "c7", length_db = hs_length)

## Before conversion: 150, after conversion: 151.

## Warning: non-unique value when setting 'row.names': '51326'

## Error in `.rowNamesDF<-`(x, value = value): duplicate 'row.names' are not allowed

mono_down_goseq_msig <- goseq_msigdb(sig_genes = downs, signatures = broad_c7,
                                     signature_category = "c7", length_db = hs_length)

## Before conversion: 296, after conversion: 293.

## Before conversion: 227921, after conversion: 35341.

## Found 281 go_db genes and 293 length_db genes out of 293.

5.2.2.6 Plot of similar experiments

## Monocyte genes with increased expression in the failed samples
## share genes with the following experiments
mono_up_goseq_msig[["pvalue_plots"]][["mfp_plot_over"]]

## Error in eval(expr, envir, enclos): object 'mono_up_goseq_msig' not found

## Monocyte genes with increased expression in the cured samples
## share genes with the following experiments
mono_down_goseq_msig[["pvalue_plots"]][["mfp_plot_over"]]

5.2.3 Evaluate Neutrophil samples

neut <- subset_expt(hs_valid, subset = "typeofcells=='Neutrophils'") %>%
  set_expt_conditions(fact = "clinicaloutcome") %>%
  set_expt_batches(fact = "donor") %>%
  set_expt_colors(colors = chosen_colors)

## There were 116, now there are 21 samples.

save_result <- save(neut, file = "rda/neutrophil_expt.rda")
neut_norm <- sm(normalize_expt(neut, convert = "cpm", filter = TRUE, transform = "log2"))
plt <- plot_pca(neut_norm, plot_labels = FALSE)$plot
pp(file = glue("images/neut_pca_normalized-v{ver}.pdf"), image = plt)

neut_nb <- sm(normalize_expt(neut, convert = "cpm", filter = TRUE,
                             transform = "log2", batch = "svaseq"))
plt <- plot_pca(neut_nb, plot_labels = FALSE)$plot
pp(file = glue("images/neut_pca_normalized_svaseq-v{ver}.pdf"), image = plt)

5.2.4 DE of Netrophil samples

5.2.4.1 Without sva

neut_de <- sm(all_pairwise(neut, model_batch = FALSE, filter = TRUE))

neut_tables <- sm(combine_de_tables(
    neut_de, keepers = keepers,
    excel = glue::glue("excel/neutrophil_clinical_all_tables-v{ver}.xlsx")))
written <- write_xlsx(data = neut_tables[["data"]][[1]],
                      excel = glue::glue("excel/neutrophil_clinical_table-v{ver}.xlsx"))
neut_sig <- sm(extract_significant_genes(neut_tables, according_to = "deseq"))
written <- write_xlsx(data = neut_sig[["deseq"]][["ups"]][[1]],
                      excel = glue::glue("excel/neutrophil_clinical_sigup-v{ver}.xlsx"))
written <- write_xlsx(data = neut_sig[["deseq"]][["downs"]][[1]],
                      excel = glue::glue("excel/neutrophil_clinical_sigdown-v{ver}.xlsx"))

neut_pct_sig <- sm(extract_significant_genes(neut_tables, n = 200, lfc = NULL,
                                             p = NULL, according_to = "deseq"))
written <- write_xlsx(data = neut_sig[["deseq"]][["ups"]][[1]],
                      excel = glue::glue("excel/neutrophil_clinical_sigup_pct-v{ver}.xlsx"))
written <- write_xlsx(data = neut_sig[["deseq"]][["downs"]][[1]],
                      excel = glue::glue("excel/neutrophil_clinical_sigdown_pct-v{ver}.xlsx"))
neut_cpm <- sm(normalize_expt(neut, convert = "cpm"))
written <- write_xlsx(data = exprs(neut_cpm),
                      excel = glue::glue("excel/neutrophil_cpm_before_batch-v{ver}.xlsx"))
neut_bcpm <- sm(normalize_expt(neut, filter = TRUE, batch = "svaseq", convert = "cpm"))
written <- write_xlsx(data = exprs(neut_bcpm),
                      excel = glue::glue("excel/neutrophil_cpm_after_batch-v{ver}.xlsx"))

5.2.4.2 With sva

neut_de_sva <- sm(all_pairwise(neut, model_batch = "svaseq", filter = TRUE))

neut_tables_sva <- sm(combine_de_tables(
    neut_de_sva, keepers = keepers,
    excel = glue::glue("excel/neutrophil_clinical_all_tables_sva-v{ver}.xlsx")))
neut_sig_sva <- sm(extract_significant_genes(
    neut_tables_sva,
    excel = glue::glue("excel/neutrophil_clinical_sig_tables_sva-v{ver}.xlsx"),
    according_to = "deseq"))

5.2.4.3 Neutrophil DE plots

## DESeq2 MA plot of failure / cure
neut_tables[["plots"]][["fail_vs_cure"]][["deseq_ma_plots"]]$plot

## DESeq2 Volcano plot of failure / cure
neut_tables[["plots"]][["fail_vs_cure"]][["deseq_vol_plots"]]$plot

## DESeq2 MA plot of failure / cure with sva
neut_tables_sva[["plots"]][["fail_vs_cure"]][["deseq_ma_plots"]]$plot

## DESeq2 Volcano plot of failure / cure with sva
neut_tables_sva[["plots"]][["fail_vs_cure"]][["deseq_vol_plots"]]$plot

5.2.4.4 Neutrophil ontology search

ups <- neut_sig[["deseq"]][["ups"]][["fail_vs_cure"]]
downs <- neut_sig[["deseq"]][["downs"]][["fail_vs_cure"]]

neut_up_gp <- simple_gprofiler(sig_genes = ups)

## Performing gProfiler GO search of 337 genes against hsapiens.

## GO search found 33 hits.

## Performing gProfiler KEGG search of 337 genes against hsapiens.

## KEGG search found 6 hits.

## Performing gProfiler REAC search of 337 genes against hsapiens.

## REAC search found 62 hits.

## Performing gProfiler MI search of 337 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 337 genes against hsapiens.

## TF search found 52 hits.

## Performing gProfiler CORUM search of 337 genes against hsapiens.

## CORUM search found 9 hits.

## Performing gProfiler HP search of 337 genes against hsapiens.

## HP search found 0 hits.

neut_up_gp[["pvalue_plots"]][["bpp_plot_over"]]

neut_up_gp[["pvalue_plots"]][["mfp_plot_over"]]

neut_up_gp[["pvalue_plots"]][["reactome_plot_over"]]

neut_down_gp <- simple_gprofiler(downs)

## Performing gProfiler GO search of 296 genes against hsapiens.

## GO search found 4 hits.

## Performing gProfiler KEGG search of 296 genes against hsapiens.

## KEGG search found 1 hits.

## Performing gProfiler REAC search of 296 genes against hsapiens.

## REAC search found 0 hits.

## Performing gProfiler MI search of 296 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 296 genes against hsapiens.

## TF search found 73 hits.

## Performing gProfiler CORUM search of 296 genes against hsapiens.

## CORUM search found 0 hits.

## Performing gProfiler HP search of 296 genes against hsapiens.

## HP search found 3 hits.

neut_down_gp[["pvalue_plots"]][["bpp_plot_over"]]

## NULL

neut_down_gp[["pvalue_plots"]][["mfp_plot_over"]]

neut_down_gp[["pvalue_plots"]][["reactome_plot_over"]]

## NULL

5.2.4.5 Neutrophil GSVA query

neut_up_goseq_msig <- goseq_msigdb(sig_genes = ups, signatures = broad_c7,
                                   signature_category = "c7", length_db = hs_length)

## Before conversion: 337, after conversion: 332.

## Before conversion: 227921, after conversion: 35341.

## Found 315 go_db genes and 332 length_db genes out of 332.

neut_down_goseq_msig <- goseq_msigdb(sig_genes = downs, signatures = broad_c7,
                                     signature_category = "c7", length_db = hs_length)

## Before conversion: 296, after conversion: 289.

## Before conversion: 227921, after conversion: 35341.

## Found 282 go_db genes and 289 length_db genes out of 289.

5.2.4.6 Plot of similar experiments

## Neutrophil genes with increased expression in the failed samples
## share genes with the following experiments
neut_up_goseq_msig[["pvalue_plots"]][["mfp_plot_over"]]

## Neutrophil genes with increased expression in the cured samples
## share genes with the following experiments
neut_down_goseq_msig[["pvalue_plots"]][["mfp_plot_over"]]

5.3 Eosinophils

5.3.1 Evaluate Eosinophil samples

eo <- subset_expt(hs_valid, subset = "typeofcells=='Eosinophils'") %>%
  set_expt_conditions(fact = "clinicaloutcome") %>%
  set_expt_batches(fact = "donor") %>%
  set_expt_colors(colors = chosen_colors)

## There were 116, now there are 15 samples.

save_result <- save(eo, file = "rda/eosinophil_expt.rda")
eo_norm <- sm(normalize_expt(eo, convert = "cpm", transform = "log2",
                             norm = "quant", filter = TRUE))
plt <- plot_pca(eo_norm, plot_labels = FALSE)$plot
pp(file = glue("images/eo_pca_normalized-v{ver}.pdf"), image = plt)

eo_nb <- sm(normalize_expt(eo, convert = "cpm", transform = "log2",
                           filter = TRUE, batch = "svaseq"))
plot_pca(eo_nb)$plot

5.3.2 DE of Eosinophil samples

5.3.2.1 Withouth sva

eo_de <- sm(all_pairwise(eo, model_batch = FALSE, filter = TRUE))

eo_tables <- sm(combine_de_tables(
    eo_de, keepers = keepers,
    excel = glue::glue("excel/eosinophil_clinical_all_tables-v{ver}.xlsx")))
written <- write_xlsx(data = eo_tables[["data"]][[1]],
                      excel = glue::glue("excel/eosinophil_clinical_table-v{ver}.xlsx"))
eo_sig <- sm(extract_significant_genes(eo_tables, according_to = "deseq"))
written <- write_xlsx(data = eo_sig[["deseq"]][["ups"]][[1]],
                      excel = glue::glue("excel/eosinophil_clinical_sigup-v{ver}.xlsx"))
written <- write_xlsx(data = eo_sig[["deseq"]][["downs"]][[1]],
                      excel = glue::glue("excel/eosinophil_clinical_sigdown-v{ver}.xlsx"))

eo_pct_sig <- sm(extract_significant_genes(eo_tables, n = 200,
                                           lfc = NULL, p = NULL, according_to = "deseq"))
written <- write_xlsx(data = eo_pct_sig[["deseq"]][["ups"]][[1]],
                      excel = glue::glue("excel/eosinophil_clinical_sigup_pct-v{ver}.xlsx"))
written <- write_xlsx(data = eo_pct_sig[["deseq"]][["downs"]][[1]],
                      excel = glue::glue("excel/eosinophil_clinical_sigdown_pct-v{ver}.xlsx"))

eo_cpm <- sm(normalize_expt(eo, convert = "cpm"))
written <- write_xlsx(data = exprs(eo_cpm),
                      excel = glue::glue("excel/eosinophil_cpm_before_batch-v{ver}.xlsx"))
eo_bcpm <- sm(normalize_expt(eo, filter = TRUE, batch = "svaseq", convert = "cpm"))
written <- write_xlsx(data = exprs(eo_bcpm),
                      excel = glue::glue("excel/eosinophil_cpm_after_batch-v{ver}.xlsx"))

5.3.2.2 With sva

eo_de_sva <- sm(all_pairwise(eo, model_batch = "svaseq", filter = TRUE))

eo_tables_sva <- sm(combine_de_tables(
    eo_de_sva, keepers = keepers,
    excel = glue::glue("excel/eosinophil_clinical_all_tables_sva-v{ver}.xlsx")))
eo_sig_sva <- sm(extract_significant_genes(
    eo_tables_sva,
    excel = glue::glue("excel/eosinophil_clinical_sig_tables_sva-v{ver}.xlsx"),
    according_to = "deseq"))

5.3.2.3 Eosinophil DE plots

## DESeq2 MA plot of failure / cure
eo_tables[["plots"]][["fail_vs_cure"]][["deseq_ma_plots"]]$plot

## DESeq2 Volcano plot of failure / cure
eo_tables[["plots"]][["fail_vs_cure"]][["deseq_vol_plots"]]$plot

## DESeq2 MA plot of failure / cure with sva
eo_tables_sva[["plots"]][["fail_vs_cure"]][["deseq_ma_plots"]]$plot

## DESeq2 Volcano plot of failure / cure with sva
eo_tables_sva[["plots"]][["fail_vs_cure"]][["deseq_vol_plots"]]$plot

5.3.2.4 Eosinophil ontology search

ups <- eo_sig[["deseq"]][["ups"]][["fail_vs_cure"]]
downs <- eo_sig[["deseq"]][["downs"]][["fail_vs_cure"]]

eo_up_gp <- simple_gprofiler(sig_genes = ups)

## Performing gProfiler GO search of 240 genes against hsapiens.

## GO search found 70 hits.

## Performing gProfiler KEGG search of 240 genes against hsapiens.

## KEGG search found 7 hits.

## Performing gProfiler REAC search of 240 genes against hsapiens.

## REAC search found 10 hits.

## Performing gProfiler MI search of 240 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 240 genes against hsapiens.

## TF search found 48 hits.

## Performing gProfiler CORUM search of 240 genes against hsapiens.

## CORUM search found 1 hits.

## Performing gProfiler HP search of 240 genes against hsapiens.

## HP search found 0 hits.

eo_up_gp[["pvalue_plots"]][["bpp_plot_over"]]

eo_up_gp[["pvalue_plots"]][["mfp_plot_over"]]

eo_up_gp[["pvalue_plots"]][["reactome_plot_over"]]

eo_down_gp <- simple_gprofiler(downs)

## Performing gProfiler GO search of 208 genes against hsapiens.

## GO search found 16 hits.

## Performing gProfiler KEGG search of 208 genes against hsapiens.

## KEGG search found 2 hits.

## Performing gProfiler REAC search of 208 genes against hsapiens.

## REAC search found 0 hits.

## Performing gProfiler MI search of 208 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 208 genes against hsapiens.

## TF search found 9 hits.

## Performing gProfiler CORUM search of 208 genes against hsapiens.

## CORUM search found 0 hits.

## Performing gProfiler HP search of 208 genes against hsapiens.

## HP search found 0 hits.

eo_down_gp[["pvalue_plots"]][["bpp_plot_over"]]

eo_down_gp[["pvalue_plots"]][["mfp_plot_over"]]

## NULL

eo_down_gp[["pvalue_plots"]][["reactome_plot_over"]]

## NULL

5.3.2.5 Eosinophil MSigDB query

eo_up_goseq_msig <- goseq_msigdb(sig_genes = ups, signatures = broad_c7,
                                 signature_category = "c7", length_db = hs_length)

## Before conversion: 240, after conversion: 237.

## Before conversion: 227921, after conversion: 35341.

## Found 222 go_db genes and 237 length_db genes out of 237.

eo_down_goseq_msig <- goseq_msigdb(sig_genes = downs, signatures = broad_c7,
                                   signature_category = "c7", length_db = hs_length)

## Before conversion: 208, after conversion: 207.

## Before conversion: 227921, after conversion: 35341.

## Found 196 go_db genes and 207 length_db genes out of 207.

5.3.2.6 Plot of similar experiments

## Eosinophil genes with increased expression in the failed samples
## share genes with the following experiments
eo_up_goseq_msig[["pvalue_plots"]][["mfp_plot_over"]]

## Eosinophil genes with increased expression in the cured samples
## share genes with the following experiments
eo_down_goseq_msig[["pvalue_plots"]][["mfp_plot_over"]]

5.4 Biopsies

5.4.1 Evaluate Biopsy samples

biop <- subset_expt(hs_valid, subset = "typeofcells=='Biopsy'") %>%
  set_expt_conditions(fact = "clinicaloutcome") %>%
  set_expt_batches(fact = "donor") %>%
  set_expt_colors(colors = chosen_colors)

## There were 116, now there are 40 samples.

save_result <- save(biop, file = "rda/biopsy_expt.rda")
biop_norm <- normalize_expt(biop, filter = TRUE, convert = "cpm",
                            transform = "log2", norm = "quant")

## Removing 5816 low-count genes (14125 remaining).

## transform_counts: Found 5 values equal to 0, adding 1 to the matrix.

plt <- plot_pca(biop_norm, plot_labels = FALSE)$plot

## Error in data.frame(sampleid = as.character(design[["sampleid"]]), condition = design[[cond_column]], : arguments imply differing number of rows: 40, 11

pp(file = glue("images/biop_pca_normalized-v{ver}.pdf"), image = plt)

biop_nb <- sm(normalize_expt(biop, convert = "cpm", filter = TRUE,
                             transform = "log2", batch = "svaseq"))
plt <- plot_pca(biop_nb, plot_labels = FALSE)$plot

## Error in data.frame(sampleid = as.character(design[["sampleid"]]), condition = design[[cond_column]], : arguments imply differing number of rows: 40, 11

pp(file = glue("images/biop_pca_normalized_svaseq-v{ver}.pdf"), image = plt)

5.4.2 DE of Biopsy samples

5.4.2.1 Without sva

biop_de <- sm(all_pairwise(biop, model_batch = FALSE, filter = TRUE))

## Error in data.frame(sampleid = as.character(design[["sampleid"]]), condition = design[[cond_column]], : arguments imply differing number of rows: 40, 11

biop_tables <- combine_de_tables(biop_de, keepers = keepers,
                                 excel = glue::glue("excel/biopsy_clinical_all_tables-v{ver}.xlsx"))

## Deleting the file excel/biopsy_clinical_all_tables-v202104.xlsx before writing the tables.

## Error in combine_de_tables(biop_de, keepers = keepers, excel = glue::glue("excel/biopsy_clinical_all_tables-v{ver}.xlsx")): object 'biop_de' not found

written <- write_xlsx(data = biop_tables[["data"]][[1]],
                      excel = glue::glue("excel/biopsy_clinical_table-v{ver}.xlsx"))

## Error in write_xlsx(data = biop_tables[["data"]][[1]], excel = glue::glue("excel/biopsy_clinical_table-v{ver}.xlsx")): object 'biop_tables' not found

biop_sig <- extract_significant_genes(biop_tables, according_to = "deseq")

## Error in extract_significant_genes(biop_tables, according_to = "deseq"): object 'biop_tables' not found

##written <- write_xlsx(data = biop_sig[["deseq"]][["ups"]][[1]],
##                      excel = glue::glue("excel/biopsy_clinical_sigup-v{ver}.xlsx"))
written <- write_xlsx(data = biop_sig[["deseq"]][["downs"]][[1]],
                      excel = glue::glue("excel/biopsy_clinical_sigdown-v{ver}.xlsx"))

## Error in write_xlsx(data = biop_sig[["deseq"]][["downs"]][[1]], excel = glue::glue("excel/biopsy_clinical_sigdown-v{ver}.xlsx")): object 'biop_sig' not found

biop_pct_sig <- extract_significant_genes(biop_tables, n = 200, lfc = NULL, p = NULL, according_to = "deseq")

## Error in extract_significant_genes(biop_tables, n = 200, lfc = NULL, p = NULL, : object 'biop_tables' not found

written <- write_xlsx(data = biop_pct_sig[["deseq"]][["ups"]][[1]],
                      excel = glue::glue("excel/biopsy_clinical_sigup_pct-v{ver}.xlsx"))

## Error in write_xlsx(data = biop_pct_sig[["deseq"]][["ups"]][[1]], excel = glue::glue("excel/biopsy_clinical_sigup_pct-v{ver}.xlsx")): object 'biop_pct_sig' not found

written <- write_xlsx(data = biop_pct_sig[["deseq"]][["downs"]][[1]],
                      excel = glue::glue("excel/biopsy_clinical_sigdown_pct-v{ver}.xlsx"))

## Error in write_xlsx(data = biop_pct_sig[["deseq"]][["downs"]][[1]], excel = glue::glue("excel/biopsy_clinical_sigdown_pct-v{ver}.xlsx")): object 'biop_pct_sig' not found

biop_cpm <- sm(normalize_expt(biop, convert = "cpm"))
written <- write_xlsx(data = exprs(biop_cpm),
                      excel = glue::glue("excel/biopsy_cpm_before_batch-v{ver}.xlsx"))
biop_bcpm <- sm(normalize_expt(biop, filter = TRUE, batch = "svaseq", convert = "cpm"))
written <- write_xlsx(data = exprs(biop_bcpm),
                      excel = glue::glue("excel/biopsy_cpm_after_batch-v{ver}.xlsx"))

5.4.2.2 with sva

biop_de_sva <- sm(all_pairwise(biop, model_batch = "svaseq", filter = TRUE))

## Error in data.frame(sample = samples, factor = as.integer(as.factor(expt[["design"]][[batch_column]])), : arguments imply differing number of rows: 40, 11, 1

biop_tables_sva <- sm(combine_de_tables(
    biop_de_sva, keepers = keepers,
    excel = glue::glue("excel/biopsy_clinical_all_tables_sva-v{ver}.xlsx")))

## Error in combine_de_tables(biop_de_sva, keepers = keepers, excel = glue::glue("excel/biopsy_clinical_all_tables_sva-v{ver}.xlsx")): object 'biop_de_sva' not found

biop_sig_sva <- sm(extract_significant_genes(
    biop_tables_sva,
    excel = glue::glue("excel/biopsy_clinical_sig_tables_sva-v{ver}.xlsx"),
    according_to = "deseq"))

## Error in extract_significant_genes(biop_tables_sva, excel = glue::glue("excel/biopsy_clinical_sig_tables_sva-v{ver}.xlsx"), : object 'biop_tables_sva' not found

5.4.2.3 Biopsy DE plots

## DESeq2 MA plot of failure / cure
biop_tables[["plots"]][["fail_vs_cure"]][["deseq_ma_plots"]]$plot

## Error in eval(expr, envir, enclos): object 'biop_tables' not found

## DESeq2 Volcano plot of failure / cure
biop_tables[["plots"]][["fail_vs_cure"]][["deseq_vol_plots"]]$plot

## Error in eval(expr, envir, enclos): object 'biop_tables' not found

## DESeq2 MA plot of failure / cure
biop_tables_sva[["plots"]][["fail_vs_cure"]][["deseq_ma_plots"]]$plot

## Error in eval(expr, envir, enclos): object 'biop_tables_sva' not found

## DESeq2 Volcano plot of failure / cure
biop_tables_sva[["plots"]][["fail_vs_cure"]][["deseq_vol_plots"]]$plot

## Error in eval(expr, envir, enclos): object 'biop_tables_sva' not found

6 Monocytes by visit

1. Can you please share with us a PCA (SVA and non-SVA) of the monocytes of the TMRC3 project, but labeling them based on the visit (V1, V2, V3)?
2. Can you please share DE lists of V1 vs V2, V1 vs V3, V1 vs. V2+V3 and V2 vs V3?

visit_colors <- chosen_colors <- c("#D95F02", "#7570B3", "#1B9E77")
names(visit_colors) <- c(1, 2, 3)
mono_visit <- subset_expt(hs_valid, subset = "typeofcells=='Monocytes'") %>%
  set_expt_conditions(fact = "visitnumber") %>%
  set_expt_batches(fact = "clinicaloutcome") %>%
  set_expt_colors(colors = chosen_colors)

## There were 116, now there are 22 samples.

mono_visit_norm <- normalize_expt(mono_visit, filter = TRUE, norm = "quant", convert = "cpm",
                                  transform = "log2")

## Removing 8966 low-count genes (10975 remaining).

## transform_counts: Found 9 values equal to 0, adding 1 to the matrix.

mono_visit_pca <- plot_pca(mono_visit_norm)

## Potentially check over the experimental design, there appear to be missing values.

## Warning in plot_pca(mono_visit_norm): There are NA values in the component data.
## This can lead to weird plotting errors.

pp(file = "images/monocyte_by_visit.png", image = mono_visit_pca$plot)

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

mono_visit_nb <- normalize_expt(mono_visit, filter = TRUE, convert = "cpm",
                                batch = "svaseq", transform = "log2")

## Removing 8966 low-count genes (10975 remaining).

## batch_counts: Before batch/surrogate estimation, 1202 entries are x==0: 0%.

## batch_counts: Before batch/surrogate estimation, 13643 entries are 0<x<1: 6%.

## Setting 333 low elements to zero.

## transform_counts: Found 333 values equal to 0, adding 1 to the matrix.

mono_visit_nb_pca <- plot_pca(mono_visit_nb)

## Potentially check over the experimental design, there appear to be missing values.

## Warning in plot_pca(mono_visit_nb): There are NA values in the component data.
## This can lead to weird plotting errors.

pp(file = "images/monocyte_by_visit_nb.png", image = mono_visit_nb_pca$plot)

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

table(pData(mono_visit_norm)$batch)

## 
##    cure failure    lost 
##       6      10       5

keepers <- list(
    "second_vs_first" = c("c2", "c1"),
    "third_vs_second" = c("c3", "c2"),
    "third_vs_first" = c("c3", "c1"))
mono_visit_de <- all_pairwise(mono_visit, model_batch = "svaseq", filter = TRUE)

## batch_counts: Before batch/surrogate estimation, 1202 entries are x==0: 0%.

## Plotting a PCA before surrogate/batch inclusion.

## Potentially check over the experimental design, there appear to be missing values.

## Warning in plot_pca(pre_batch, plot_labels = FALSE, ...): There are NA values in
## the component data. This can lead to weird plotting errors.

## Using svaseq to visualize before/after batch inclusion.

## Performing a test normalization with: raw

## Removing 0 low-count genes (10975 remaining).

## batch_counts: Before batch/surrogate estimation, 1202 entries are x==0: 0%.

## batch_counts: Before batch/surrogate estimation, 13643 entries are 0<x<1: 6%.

## Setting 333 low elements to zero.

## transform_counts: Found 333 values equal to 0, adding 1 to the matrix.

## Potentially check over the experimental design, there appear to be missing values.

## Warning in plot_pca(post_batch, plot_labels = FALSE, ...): There are NA values
## in the component data. This can lead to weird plotting errors.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

mono_visit_tables <- combine_de_tables(
    mono_visit_de,
    keepers = keepers,
    excel = glue::glue("excel/mono_visit_tables-v{ver}.xlsx"))

## Deleting the file excel/mono_visit_tables-v202104.xlsx before writing the tables.

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

new_factor <- pData(mono_visit)[["condition"]]
not_one_idx <- new_factor != 1
new_factor[not_one_idx] <- "not_1"
mono_one_vs <- set_expt_conditions(mono_visit, new_factor)

mono_one_vs_de <- all_pairwise(mono_one_vs, model_batch = "svaseq", filter = TRUE)

## batch_counts: Before batch/surrogate estimation, 1202 entries are x==0: 0%.

## Plotting a PCA before surrogate/batch inclusion.

## Potentially check over the experimental design, there appear to be missing values.

## Warning in plot_pca(pre_batch, plot_labels = FALSE, ...): There are NA values in
## the component data. This can lead to weird plotting errors.

## Using svaseq to visualize before/after batch inclusion.

## Performing a test normalization with: raw

## Removing 0 low-count genes (10975 remaining).

## batch_counts: Before batch/surrogate estimation, 1202 entries are x==0: 0%.

## batch_counts: Before batch/surrogate estimation, 13643 entries are 0<x<1: 6%.

## Setting 333 low elements to zero.

## transform_counts: Found 333 values equal to 0, adding 1 to the matrix.

## Potentially check over the experimental design, there appear to be missing values.

## Warning in plot_pca(post_batch, plot_labels = FALSE, ...): There are NA values
## in the component data. This can lead to weird plotting errors.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

mono_one_vs_tables <- combine_de_tables(
    mono_one_vs_de,
    excel = glue::glue("excel/mono_one_vs_tables-v{ver}.xlsx"))

## Deleting the file excel/mono_one_vs_tables-v202104.xlsx before writing the tables.

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

7 Test TSP

In writing the following, I quickly realized that tspair was not joking when it said it is intended for small numbers of genes. For a full expressionset of human data it is struggling. I like the idea, it may prove worth while to spend some time optimizing the package so that it is more usable.

expt <- hs_clinical_nobiop

simple_tsp <- function(expt, column = "condition") {
  facts <- levels(as.factor(pData(expt)[[column]]))
  retlist <- list()
  if (length(facts) < 2) {
    stop("This requires factors with at least 2 levels.")
  } else if (length(facts) == 2) {
    retlist <- simple_tsp_pair(expt, column = column)
  } else {
    for (first in 1:(length(facts) - 1)) {
      for (second in 2:(length(facts))) {
        if (first < second) {
          name <- glue::glue("{facts[first]}_vs_{facts[second]}")
          message("Starting ", name, ".")
          substring <- glue::glue("{column}=='{facts[first]}'|{column}=='{facts[second]}'")
          subby <- subset_expt(expt, subset=as.character(substring))
          retlist[[name]] <- simple_tsp_pair(subby, column = column)
        }
      }
    }
  }
}

simple_tsp_pair <- function(subby, column = "condition", repetitions = 50) {
  tsp_input <- subby[["expressionset"]]
  tsp_output <- tspcalc(tsp_input, column)
  tsp_scores <- tspsig(tsp_input, column, B = repetitions)
}

tsp1 <- tspcalc(tsp_input, "condition")

if (!isTRUE(get0("skip_load"))) {
  pander::pander(sessionInfo())
  message(paste0("This is hpgltools commit: ", get_git_commit()))
  message(paste0("Saving to ", savefile))
  tmp <- sm(saveme(filename = savefile))
}

## If you wish to reproduce this exact build of hpgltools, invoke the following:

## > git clone http://github.com/abelew/hpgltools.git

## > git reset 53433f808ad055552025c90161db331405085a9e

## This is hpgltools commit: Tue May 4 12:44:03 2021 -0400: 53433f808ad055552025c90161db331405085a9e

## Saving to tmrc3_02sample_estimation_v202104.rda.xz

tmp <- loadme(filename = savefile)