L. panamensis 202205: Reorganizing TMRC3 Analyses

1 Introduction

The set of analyses performed in tmrc3 sample estimation has become unorganized and difficult to follow. This is intended to start at the level of broad organization before stepping into the analyses.

As of 20220224, this document is attempting to synchronize with the google doc named ‘TMRC3_Aug18_2021’. I am hoping therefore to have blocks named according to the figures/etc in that document.

1.1 Goals

These samples are from patients who either successfully cleared a Leishmania panamensis infection following treatment, or did not. They include biopsies from each patient along with purifications for Monocytes, Neutrophils, and Eosinophils. When possible, this process was repeated over three visits; but some patients did not return for the second or third visit.

The over-arching goal is to look for attributes(most likely genes) which distinguish patients who do and do not cure the infection after treatment. If possible, these will be apparent on the first visit.

1.2 Relevant Metadata

The metadata factors (in no particular order) of the experiment are:

1. Visit. 1, 2, or 3.
2. Patient.
3. Clinical outcome: cure or fail. This includes the population of patients who did not return and are labeled ‘lost’. This is, by definition, confounded with patient.
4. Drug treatment. Most of the patients were treated with an antimonial, but a few were treated with miltefosine.
5. Cell type or biopsy.

Metadata was also collected for the patients and there are many factors which have potential to affect the outcome. These analyses will generally remain agnostic for these factors, which include (among other things):

1. sex
2. ethnicity
3. age
4. ulcer/lesion attributes (by visit)
5. amount of time spent with infection

1.3 Libraries, sequencing, mapping, quantification

The samples used in these analyses were all collected, purified, and libraries generated by the scientists/doctors at CIDEIM. The sequencing libraries were generated via the TruSeq non-stranded library kit and sequenced either at JHU or UMD; earlier samples were single-ended, but most were paired.

All samples were trimmed with trimomatic using the same set of parameters. All mapping was performed with hisat2 version 2.2.1. Quantifications were also performed with salmon version 1.2.0. The reference genome used was hg38 100 (released 202004). When samples were mapped against L.panamensis, the TriTrypDB version 36 reference was used. The set of annotations were therefore limited to ensembl’s 2020 release. The early samples were actually first mapped with hg38 91 and later redone.

1.4 Types of analyses included

This document will limit itself to a set of canonical RNAseq analyses. It must therefore create files containing the raw data to facilitate sharing the data. It will plot metrics of the data to demonstrate the sequencing quality and clustering of the samples under the various conditions examined and normalizations employed. It will perform differential expression analyses for the metadata factors of interest alongside likelihood ratio tests for factors like celltype and time. Given the sets of over/under expressed genes observed in the various DE methods, it will perform the likely gene set tests for over represented gene ontology groups, reactome, etc. Simultaneously, the raw data will be passed to gene set variance analyses to see if there are groups of genes overrepresented in other experiments.

1.5 Metadata collection

There are two metadata sources:

1. The online sample sheet, which I periodically update and download into the ‘sample_sheets/’ directory.
2. The crf metadata, describing the individual patients.

samplesheet <- "sample_sheets/tmrc3_samples_202205.xlsx"
crf_metadata <- "sample_sheets/20210825_EXP_ESPECIAL_TMRC3_VERSION_2.xlsx"

2 Annotation Collection

The primary annotation sources are:

1. Ensembl’s biomart archive from 2020 for human annotations.
2. The TriTrypDB release 36 for parasite annotations.

Both provide GO data. They also provide helpful links to other data sources. For the moment, we are focusing on the human annotations.

2.1 Gene annotations

These analyses have focused on gene-level abundances/differences. Thus, when htseq-count was invoked against the hisat2-based mappings, parameters were chosen to count genes rather than transcripts. Similarly, when salmon counts were used via tximport, a mapping of genes to transcripts was used to collapse the matrix to gene-level abundances. This decision may be revisited.

hs_annot <- load_biomart_annotations(year="2020", month="jan")

## The biomart annotations file already exists, loading from it.

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:227784         Length:227784      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.:15.0   3rd Qu.: 5.00     
##                                           Max.   :28.0   Max.   :17.00     
##                                                                            
##  hgnc_symbol        description        gene_biotype         cds_length    
##  Length:227784      Length:227784      Length:227784      Min.   :     3  
##  Class :character   Class :character   Class :character   1st Qu.:   357  
##  Mode  :character   Mode  :character   Mode  :character   Median :   694  
##                                                           Mean   :  1140  
##                                                           3rd Qu.:  1446  
##                                                           Max.   :107976  
##                                                           NA's   :127222  
##  chromosome_name       strand          start_position      end_position     
##  Length:227784      Length:227784      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:227784     
##  Class :character  
##  Mode  :character  
##                    
##                    
##                    
## 

2.2 Gene ontology data

The set of GO categories has not been limited to the 2020 data at the time of this writing.

hs_go <- load_biomart_go()[["go"]]

## The biomart annotations file already exists, loading from it.

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

3 The full dataset

Before we do any of the following subsets/analyses of the data, we need to collect it all in one place. Let’s do that here and name it ‘hs_valid’ meaning it is the set of valid data for Homo sapiens.

sanitize_columns <- c("visitnumber", "clinicaloutcome", "donor",
                      "typeofcells", "clinicalpresentation", "drug",
                      "condition", "batch")
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", meta_column="ncrna_lost") %>%
  sanitize_expt_metadata(columns=sanitize_columns) %>%
  set_expt_factors(columns=sanitize_columns, class="factor") %>%
  set_expt_conditions(fact="clinicaloutcome") %>%
  set_expt_batches(fact="visitnumber")

## Reading the sample metadata.

## Dropped 69 rows from the sample metadata because the sample ID is blank.

## The sample definitions comprises: 254 rows(samples) and 84 columns(metadata fields).

## 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 21447 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 features and 246 samples.

## remove_genes_expt(), before removal, there were 21481 genes, now there are 19923.

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

## TMRC30015 TMRC30017 TMRC30019 TMRC30044 TMRC30045 TMRC30154 TMRC30269 TMRC30241 
##     79.19     85.66     89.69     80.29     73.28     83.16     89.18     89.37

## The following should make visit 1 the largest if one uses that column as a size factor when plotting.
meta <- pData(hs_expt) %>%
  mutate(visitnumber = fct_relevel(visitnumber, c("notapplicable", "3", "2", "1")))
pData(hs_expt) <- meta

The above block does the following:

1. Creates an expressionset using the ‘hg38100hisatfile’ column from the most recently downloaded sample sheet (the column’s name has any puctuation/spaces/capitals/etc removed) and the set of human annotations downloaded above.
2. This expressionset is passed to a filter which pulls out only the protein_coding genes and uses the information from that process to add a new metadata column called ‘ncrna_lost’. Thus it keeps a tally of the number of reads lost in the filter and adds it to the sample sheet.
3. It is passed again to a function which sanitizes the metadata (there were a couple of entries which said ‘cure’ instead of ‘Cure’ or vice versa) and similarly removes problematic characters.
4. Passed to a function which sets some columns explicitly to factors instead of characters.
5. Sets the experimental ‘condition’ to the factor of cure vs. fail.
6. Sets the experimental ‘batch’ to visit number.
7. Resets the levels of the visit number so that the samples which were ‘notapplicable’ are logically before visit 1 which is before 2 before 3.

3.1 Add the CRF patient metadata

Let us also merge in the clinician’s metadata. I worry a little that this might not be allowed for dbGap data, but if it is a problem I suspect we can just remove the bad columns from it. Also note that I rarely use the join function, but it is somewhat required here because I do not want to risk shuffling the metadata when I add the new metadata, which comes from a spreadsheet sorted by patient, not sample. In doing this I therefore just created a new column ‘join’ which contains the shared information, e.g. the patient ID from the existing metadata and the same ID from the CRF file which has been coerced into lowercase.

hs_pd <- pData(hs_expt)
start <- rownames(hs_pd)
hs_crf <- openxlsx::read.xlsx(crf_metadata)
hs_crf[["join"]] <- tolower(hs_crf[["codigo_paciente"]])
hs_pd[["join"]] <- hs_pd[["tubelabelorigin"]]
test <- plyr::join(hs_pd, hs_crf, by="join")
test[["join"]] <- NULL
rownames(test) <- rownames(hs_pd)
na_idx <- is.na(test)
test[na_idx] <- "undefined"
pData(hs_expt) <- test

The above block does the following:

1. Extracts the metadata and makes a character vector of the rownames (e.g. the sample IDs).
2. Creates a table of the clinical metadata from the file I downloaded.
3. Creates a ‘join’ column from both the existing metadata and the clinical metadata which is the patient code, using whatever the column names were in each sheet.
4. Joins the two tables using this information. The join function is explicitly used in this context to make sure that the row-order of the entries does not change.
5. Gets rid of any entries in the new table which are NA.
6. Replaces the metadata of the expressionset with this new, larger table.

4 Sanity check, clinical outcome vs. CRF data

In our shared online sample sheet there is a clinical outcome column which should match the final CRF column ‘ef_lc_estado_final_estudio’ with the caveat that the CRF data is numeric:

0: curacion definitiva 1: fall terapeutica 2: perdida druante el seguimiento 3: excludio durante el estudio

two_columns <- pData(hs_expt)[, c("clinicaloutcome", "ef_lc_estado_final_estudio")]
undef_idx <- two_columns[[2]] == "undefined"
two_columns <- two_columns[!undef_idx, ]
two_columns[["rewritten"]] <- "undef"
cure_idx <- two_columns[[2]] == 0
two_columns[cure_idx, "rewritten"] <- "cure"
fail_idx <- two_columns[[2]] == 1
two_columns[fail_idx, "rewritten"] <- "failure"
lost_idx <- two_columns[[2]] == 2
two_columns[lost_idx, "rewritten"] <- "lost"
same_idx <- two_columns[["clinicaloutcome"]] == two_columns[["rewritten"]]
broken <- two_columns[!same_idx, ]
broken

##           clinicaloutcome ef_lc_estado_final_estudio rewritten
## TMRC30148            cure                          1   failure
## TMRC30149            cure                          1   failure
## TMRC30138            cure                          1   failure
## TMRC30150            cure                          1   failure
## TMRC30140            cure                          1   failure
## TMRC30151            cure                          1   failure
## TMRC30176            cure                          1   failure
## TMRC30153            cure                          1   failure
## TMRC30269            cure                          1   failure

4.1 Collect sample numbers before filtering

dim(pData(hs_expt))

## [1] 246 144

table(pData(hs_expt)$drug)

## 
##    antimony miltefosine        none 
##         216           8          22

table(pData(hs_expt)$clinic)

## 
##      Cali    Tumaco undefined 
##        75       143        28

table(pData(hs_expt)$typeofcells)

## 
##      biopsy eosinophils macrophages   monocytes neutrophils       pbmcs 
##          21          45          28          74          72           6

table(pData(hs_expt)$visit)

## 
## notapplicable             3             2             1 
##            28            54            56           108

summary(as.numeric(pData(hs_expt)$eb_lc_tiempo_evolucion))

## Warning in summary(as.numeric(pData(hs_expt)$eb_lc_tiempo_evolucion)): NAs
## introduced by coercion

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    2.00    4.00    6.00    7.37    8.00   21.00      60

summary(as.numeric(pData(hs_expt)$eb_lc_tto_mcto_glucan_dosis))

## Warning in summary(as.numeric(pData(hs_expt)$eb_lc_tto_mcto_glucan_dosis)): NAs
## introduced by coercion

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    13.0    14.0    19.0    59.5    20.0   999.0      60

summary(as.numeric(pData(hs_expt)$v3_lc_ejey_lesion_mm_1))

## Warning in summary(as.numeric(pData(hs_expt)$v3_lc_ejey_lesion_mm_1)): NAs
## introduced by coercion

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     0.0     7.2    32.4   300.0   999.0   999.0      60

summary(as.numeric(pData(hs_expt)$v3_lc_lesion_area_1))

## Warning in summary(as.numeric(pData(hs_expt)$v3_lc_lesion_area_1)): NAs
## introduced by coercion

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##       0     223     999    2133    3016   11487      60

summary(as.numeric(pData(hs_expt)$v3_lc_ejex_ulcera_mm_1))

## Warning in summary(as.numeric(pData(hs_expt)$v3_lc_ejex_ulcera_mm_1)): NAs
## introduced by coercion

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     0.0     0.0    12.5   292.8   999.0   999.0      60

table(pData(hs_expt)$eb_lc_sexo)

## 
##         1         2 undefined 
##       155        31        60

table(pData(hs_expt)$eb_lc_etnia)

## 
##         1         2         3 undefined 
##        99        53        34        60

summary(as.numeric(pData(hs_expt)$edad))

## Warning in summary(as.numeric(pData(hs_expt)$edad)): NAs introduced by coercion

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    18.0    24.0    28.0    29.9    36.0    51.0      60

table(pData(hs_expt)$eb_lc_peso)

## 
##     100.8      53.9      55.9      57.9      58.1      58.3      58.6        59 
##         2         9         2         2         7        10        11        14 
##      59.6        62        63        67      69.4      74.7      75.6      76.5 
##         1         8         6         6        10         3         3         3 
##        77        78      79.2        82      83.3      83.4      86.4        87 
##        18        10        10         9         4        10         9         3 
##        89      93.3 undefined 
##         9         7        60

table(pData(hs_expt)$eb_lc_estatura)

## 
##       145       152       154       158       159       160       162       163 
##         6         1        10        17         2         6        10         9 
##       164       165       166       167       172       173       174       175 
##        15        12        19         3        10         4        30         2 
##       176       177       182       183 undefined 
##         1        10         9        10        60

table(pData(hs_expt)$clinic)

## 
##      Cali    Tumaco undefined 
##        75       143        28

table(pData(hs_expt)$ef_lc_estado_final_estudio)

## 
##         0         1         2 undefined 
##        75        69        16        86

table(pData(hs_expt)$clinicaloutcome)

## 
##          cure       failure          lost notapplicable 
##           132            62            16            36

4.2 Define desired colors for the various subsets

There are lots of ways which we will categorize the data, here are some potential color choices for them.

cf_colors <- list(
    "cure" = "#998EC3",
    "failure" = "#F1A340")
type_visit_colors <- list(
    "monocytes_v1" = "#DD1C77",
    "monocytes_v2" = "#C994C7",
    "monocytes_v3" = "#E7E1EF",
    "eosinophils_v1" = "#31A354",
    "eosinophils_v2" = "#ADDD8E",
    "eosinophils_v3" = "#F7FCD9",
    "neutrophils_v1" = "#3182BD",
    "neutrophils_v2" = "#9ECAE1",
    "neutrophils_v3" = "#DEEBF7",
    "biopsy_v1" = "#D95F0E")
type_colors <- list(
    "monocytes" = "#DD1C77",
    "eosinophils" = "#31A354",
    "neutrophils" = "#3182BD",
    "biopsy" = "#D95F0E")
visit_colors <- list()
cf_type_colors <- list(
    cure_biopsy = "#D95F0E",
    failure_biopsy = "#FEC44F",
    cure_monocytes = "#DD1C77",
    failure_monocytes = "#C994C7",
    cure_eosinophils = "#31A354",
    failure_eosinophils = "#ADDD8E",
    cure_neutrophils = "#3182BD",
    failure_neutrophils = "#9ECAE1")

Note from Maria Adelaida: Some chemokines are suggestive of Eosinophil recruitment.

4.3 Define the starting data

The following block contains the primary dataset which is the parent of everything which follows.

4.3.1 Set our initial coverage goal

There exists a baseline coverage below which we do not wish to fall. One likely way to approach it heuristically is to assume we should observe some number of genes. With that in mind, I arbitrarily chose 11,000 non-zero genes as the minimum.

With this in mind, here is a non-zero plot before a cutoff of 11,000 genes.

all_nz <- plot_nonzero(hs_expt)
all_nz$plot

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

To my eyes, there are 3 or 4 samples which are likely candidates for removal.

hs_valid <- subset_expt(hs_expt, nonzero=11000) %>%
  subset_expt(subset="clinicaloutcome!='lost'") %>%
  subset_expt(subset="clinicaloutcome!='notapplicable'") %>%
  subset_expt(subset="clinicaloutcome!='null'") %>%
  set_expt_colors(cf_colors)

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

## TMRC30010 TMRC30050 TMRC30052 
##     52429    807571   3086349

## subset_expt(): There were 246, now there are 243 samples.

## subset_expt(): There were 243, now there are 228 samples.

## subset_expt(): There were 228, now there are 192 samples.

## subset_expt(): There were 192, now there are 192 samples.

4.3.2 After removal

Here is the distribution of samples after dropping the less than 11,000 gene-samples.

nz_post <- plot_nonzero(hs_valid)
nz_post$plot

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

4.4 Count up sample types

I recently convinced myself that there is a difference between the data when it does and does not include the miltefosine treated patients.

However, when I actually counted up the samples by a few criteria I immediately relized this is spurious.

table(pData(hs_valid)$drug)

## 
##    antimony miltefosine 
##         184           8

table(pData(hs_valid)$clinic)

## 
##   Cali Tumaco 
##     61    131

table(pData(hs_valid)$clinicaloutcome)

## 
##    cure failure 
##     132      60

table(pData(hs_valid)$typeofcells)

## 
##      biopsy eosinophils   monocytes neutrophils 
##          19          42          66          65

table(pData(hs_valid)$visit)

## 
##  3  2  1 
## 51 50 91

summary(as.numeric(pData(hs_valid)$eb_lc_tiempo_evolucion))

## Warning in summary(as.numeric(pData(hs_valid)$eb_lc_tiempo_evolucion)): NAs
## introduced by coercion

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    2.00    4.00    6.00    7.14    8.00   21.00      24

summary(as.numeric(pData(hs_valid)$eb_lc_tto_mcto_glucan_dosis))

## Warning in summary(as.numeric(pData(hs_valid)$eb_lc_tto_mcto_glucan_dosis)): NAs
## introduced by coercion

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    13.0    15.0    19.0    64.3    20.0   999.0      24

summary(as.numeric(pData(hs_valid)$v3_lc_ejey_lesion_mm_1))

## Warning in summary(as.numeric(pData(hs_valid)$v3_lc_ejey_lesion_mm_1)): NAs
## introduced by coercion

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     0.0     7.2    32.7   328.6   999.0   999.0      24

summary(as.numeric(pData(hs_valid)$v3_lc_lesion_area_1))

## Warning in summary(as.numeric(pData(hs_valid)$v3_lc_lesion_area_1)): NAs
## introduced by coercion

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##       0     222     999    1975    2448   11487      24

summary(as.numeric(pData(hs_valid)$v3_lc_ejex_ulcera_mm_1))

## Warning in summary(as.numeric(pData(hs_valid)$v3_lc_ejex_ulcera_mm_1)): NAs
## introduced by coercion

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     0.0     0.0    16.9   324.2   999.0   999.0      24

table(pData(hs_valid)$eb_lc_sexo)

## 
##         1         2 undefined 
##       143        25        24

table(pData(hs_valid)$eb_lc_etnia)

## 
##         1         2         3 undefined 
##        97        37        34        24

summary(as.numeric(pData(hs_valid)$edad))

## Warning in summary(as.numeric(pData(hs_valid)$edad)): NAs introduced by coercion

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    18.0    24.0    28.0    29.6    35.0    51.0      24

table(pData(hs_valid)$eb_lc_peso)

## 
##     100.8      53.9      55.9      57.9      58.1      58.3      58.6        59 
##         2         9         2         2         7        10         1         8 
##      59.6        62        63        67      69.4      74.7      75.6      76.5 
##         1         6         6         6        10         3         3         3 
##        77        78      79.2        82      83.3      83.4      86.4        87 
##        18        10        10         9         4        10         9         3 
##        89      93.3 undefined 
##         9         7        24

table(pData(hs_valid)$eb_lc_estatura)

## 
##       152       154       158       159       160       163       164       165 
##         1        10        15         2         6         9        15        12 
##       166       167       172       173       174       175       176       177 
##        19         3        10         4        30         2         1        10 
##       182       183 undefined 
##         9        10        24

table(pData(hs_valid)$clinic)

## 
##   Cali Tumaco 
##     61    131

The above block does what it says on the tin, subsets the data to exclude any sample with less than 11,000 observed genes.

4.5 Set up initial data subsets of interest

One of the first global metrics I would like to provide is the set of library sizes. Unfortunately, we have too many samples to fit on a screen now. Therefore, I am going to do an early split of the data by cell type in the hopes that doing so will make it possible to visualize the library sizes/nonzero genes.

The initial factor for this is ‘typeofcells’.

table(pData(hs_valid)[["typeofcells"]])

## 
##      biopsy eosinophils   monocytes neutrophils 
##          19          42          66          65

biopsy_samples <- subset_expt(hs_valid, subset="typeofcells=='biopsy'")

## subset_expt(): There were 192, now there are 19 samples.

eosinophil_samples <- subset_expt(hs_valid, subset="typeofcells=='eosinophils'")

## subset_expt(): There were 192, now there are 42 samples.

monocyte_samples <- subset_expt(hs_valid, subset="typeofcells=='monocytes'")

## subset_expt(): There were 192, now there are 66 samples.

neutrophil_samples <- subset_expt(hs_valid, subset="typeofcells=='neutrophils'")

## subset_expt(): There were 192, now there are 65 samples.

## Currently, these are not used, but instead I pulled the samples from the hs_clinical
## which means the biopsies are not included.
v1_samples <- subset_expt(hs_valid, subset="visitnumber=='1'")

## subset_expt(): There were 192, now there are 91 samples.

v1_monocytes <- subset_expt(v1_samples, subset="typeofcells=='monocytes'")

## subset_expt(): There were 91, now there are 29 samples.

v1_neutrophils <- subset_expt(v1_samples, subset="typeofcells=='neutrophils'")

## subset_expt(): There were 91, now there are 28 samples.

v1_eosinophils <- subset_expt(v1_samples, subset="typeofcells=='eosinophils'")

## subset_expt(): There were 91, now there are 15 samples.

v2_samples <- subset_expt(hs_valid, subset="visitnumber=='2'")

## subset_expt(): There were 192, now there are 50 samples.

v2_monocytes <- subset_expt(v2_samples, subset="typeofcells=='monocytes'")

## subset_expt(): There were 50, now there are 18 samples.

v2_neutrophils <- subset_expt(v2_samples, subset="typeofcells=='neutrophils'")

## subset_expt(): There were 50, now there are 18 samples.

v2_eosinophils <- subset_expt(v2_samples, subset="typeofcells=='eosinophils'")

## subset_expt(): There were 50, now there are 14 samples.

v3_samples <- subset_expt(hs_valid, subset="visitnumber=='3'")

## subset_expt(): There were 192, now there are 51 samples.

v3_monocytes <- subset_expt(v3_samples, subset="typeofcells=='monocytes'")

## subset_expt(): There were 51, now there are 19 samples.

v3_neutrophils <- subset_expt(v3_samples, subset="typeofcells=='neutrophils'")

## subset_expt(): There were 51, now there are 19 samples.

v3_eosinophils <- subset_expt(v3_samples, subset="typeofcells=='eosinophils'")

## subset_expt(): There were 51, now there are 13 samples.

The above block does a lot of subset operations to create separate data structures on a per-celltype and per-visit basis. Ergo, our large data structure is now joined by ~21 new, smaller data structures.

5 Parasite reads

Let us see if we can make an expressionset of the parasite reads in the TMRC3 samples and distinguish between the faux and real reads. E.g: Are there samples which definitively contain parasites?

lp_expt <- create_expt(samplesheet,
                       file_column="lpanamensisv36hisatfile", gene_info = NULL) %>%
  subset_expt(coverage=1000) %>%
  set_expt_conditions(fact="typeofcells")

## Reading the sample metadata.

## Dropped 69 rows from the sample metadata because the sample ID is blank.

## The sample definitions comprises: 254 rows(samples) and 84 columns(metadata fields).

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30050/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30052/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30071/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30056/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30058/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30105/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30094/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30080/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30103/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30107/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30083/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30082/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30093/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30113/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30096/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30118/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30119/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30165/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30166/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30180/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30046/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30047/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30048/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30194/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30195/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30196/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30049/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30053/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30054/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30115/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30121/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30122/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30169/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30170/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30164/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30055/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30068/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30070/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30191/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30192/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30171/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30158/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30159/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30188/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30189/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30190/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30139/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30160/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30161/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30059/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30060/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30061/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30062/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30063/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30051/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30064/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30065/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30162/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30066/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30067/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30117/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30057/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30069/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30132/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30167/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30168/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30152/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30123/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30116/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30074/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30177/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30072/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30076/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30077/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30157/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30181/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30182/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30183/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30133/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30136/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30078/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30088/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30079/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30184/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30134/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30135/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30155/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30129/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30137/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30174/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30175/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30154/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30172/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30173/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30142/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30143/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30144/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30145/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30146/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30147/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30185/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30186/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30156/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30148/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30149/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30138/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30150/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30140/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30151/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30176/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30153/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30178/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30179/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30197/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30198/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30199/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30200/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30201/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30202/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30203/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30204/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30205/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30206/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30207/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30208/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30217/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30218/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30219/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30220/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30209/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30210/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30211/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30212/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30213/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30214/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30215/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30216/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30221/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30222/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30223/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30224/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30225/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30226/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30227/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30228/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30229/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30230/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30231/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30232/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30233/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30234/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30235/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30237/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30238/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30266/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30268/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30286/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30249/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30252/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30250/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30251/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30245/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30246/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30247/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30248/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30244/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30260/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30261/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30262/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30263/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30264/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30265/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30269/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30253/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30270/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30241/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30242/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30273/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30275/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30271/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30274/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30276/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30272/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30255/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30256/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30254/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30258/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30257/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30239/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30240/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30283/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30284/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30281/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30277/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30279/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30280/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30278/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30282/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a
## multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, :
## The file: /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/rnaseq/
## lpanamensis_tmrc_2019/preprocessing/TMRC30285/outputs/03hisat2_lpanamensis_v36/
## sno_gene_gene_id.count.xz has mismatched rownames.

## Warning in create_expt(samplesheet, file_column = "lpanamensisv36hisatfile", :
## Some samples were removed when cross referencing the samples against the count
## data.

## Matched 8778 annotations and counts.

## Bringing together the count matrix and gene information.

## Warning in create_expt(samplesheet, file_column =
## "lpanamensisv36hisatfile", : The following samples have no counts!
## TMRC30010TMRC30144TMRC30145TMRC30146TMRC30185TMRC30207TMRC30217TMRC30219TMRC30220TMRC30209TMRC30210TMRC30211TMRC30212TMRC30213TMRC30214TMRC30216TMRC30221TMRC30222TMRC30223TMRC30225TMRC30226TMRC30227TMRC30228TMRC30229TMRC30230TMRC30231TMRC30233TMRC30234TMRC30235TMRC30274TMRC30272TMRC30279TMRC30278TMRC30285

## Saving the expressionset to 'expt.rda'.

## The final expressionset has 8778 features and 244 samples.

## The samples removed (and read coverage) when filtering samples with less than 1000 reads are:

## TMRC30001 TMRC30002 TMRC30003 TMRC30004 TMRC30005 TMRC30006 TMRC30007 TMRC30008 
##        12         8         9        16        25        29         3        16 
## TMRC30009 TMRC30010 TMRC30011 TMRC30012 TMRC30013 TMRC30050 TMRC30018 TMRC30118 
##        16         0         5         9        13       345       110       747 
## TMRC30119 TMRC30014 TMRC30021 TMRC30038 TMRC30023 TMRC30025 TMRC30165 TMRC30166 
##       419         4        88       208       120       954       433       496 
## TMRC30030 TMRC30031 TMRC30032 TMRC30024 TMRC30040 TMRC30033 TMRC30194 TMRC30195 
##         3         8        22        49       896        36       198       458 
## TMRC30196 TMRC30164 TMRC30037 TMRC30027 TMRC30028 TMRC30034 TMRC30035 TMRC30036 
##       687       424         9       108        93        21       153        11 
## TMRC30192 TMRC30041 TMRC30042 TMRC30043 TMRC30045 TMRC30171 TMRC30158 TMRC30159 
##       384       264       545       571       334       370       284       405 
## TMRC30189 TMRC30190 TMRC30139 TMRC30160 TMRC30161 TMRC30152 TMRC30123 TMRC30181 
##       307       537       140       296       300       495       304       613 
## TMRC30182 TMRC30155 TMRC30129 TMRC30137 TMRC30174 TMRC30154 TMRC30172 TMRC30173 
##       468       508       353       389       644       398       383       898 
## TMRC30142 TMRC30143 TMRC30144 TMRC30145 TMRC30146 TMRC30147 TMRC30185 TMRC30186 
##         4         1         0         0         0         1         0       852 
## TMRC30148 TMRC30138 TMRC30150 TMRC30140 TMRC30151 TMRC30178 TMRC30179 TMRC30197 
##       706       156       470       153       480       420       515       202 
## TMRC30198 TMRC30200 TMRC30201 TMRC30202 TMRC30203 TMRC30205 TMRC30206 TMRC30207 
##       158        96       149       145       104       169       125         0 
## TMRC30208 TMRC30217 TMRC30218 TMRC30219 TMRC30220 TMRC30209 TMRC30210 TMRC30211 
##         1         0         2         0         0         0         0         0 
## TMRC30212 TMRC30213 TMRC30214 TMRC30215 TMRC30216 TMRC30221 TMRC30222 TMRC30223 
##         0         0         0         1         0         0         0         0 
## TMRC30224 TMRC30225 TMRC30226 TMRC30227 TMRC30228 TMRC30229 TMRC30230 TMRC30231 
##         2         0         0         0         0         0         0         0 
## TMRC30232 TMRC30233 TMRC30234 TMRC30235 TMRC30238 TMRC30266 TMRC30260 TMRC30261 
##         1         0         0         0       335       822       291       500 
## TMRC30262 TMRC30263 TMRC30264 TMRC30265 TMRC30269 TMRC30270 TMRC30273 TMRC30275 
##       905       431       309       596        58        21         2         2 
## TMRC30271 TMRC30274 TMRC30276 TMRC30272 TMRC30254 TMRC30257 TMRC30239 TMRC30240 
##         4         0         1         0       208       427       209       314 
## TMRC30283 TMRC30284 TMRC30281 TMRC30277 TMRC30279 TMRC30280 TMRC30278 TMRC30282 
##        18         4       548         6         0         3         0         1 
## TMRC30285 
##         0

## subset_expt(): There were 244, now there are 99 samples.

visit_fact <- pData(lp_expt)[["visitnumber"]]
batch_na <- is.na(visit_fact)
visit_fact[batch_na] <- "undefined"
lp_expt <- set_expt_batches(lp_expt, fact = visit_fact)

lp_norm <- normalize_expt(lp_expt, filter="simple", norm="quant",
                          convert="cpm", transform="log2")

## Removing 66 low-count genes (8712 remaining).

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

plotted <- plot_pca(lp_norm, plot_labels=FALSE)
plotted$plot

plotted_3d <- plot_3d_pca(plotted)

The above block is similar in concept to the previous expressionset creation. It uses a different column and currently ignores the gene annotations. Given that many of the samples have essentially 0 reads, I set a cutoff of only 20 observed genes. Finally, I did a quick and dirty PCA plot of this peculiar data structure in the hopes of being able to ‘see’ the difference between what are assumed to be ‘real’ samples with a significant number of ‘real’ parasite reads vs. those samples which just have a couple of potentially spurious reads.

6 Distributions/Visualizations of interest

The sets of samples used to visualize the data will also comprise the sets used when later performing the various differential expression analyses.

6.1 Global metrics

Start out with some initial metrics of all samples. The most obvious are plots of the numbers of non-zero genes observed, heatmaps showing the relative relationships among the samples, the relative library sizes, and some PCA. It might be smart to split the library sizes up across subsets of the data, because they have expanded too far to see well on a computer screen.

The most likely factors to query when considering the entire dataset are cure/fail, visit, and cell type. This is the level at which we will choose samples to exclude from future analyses.

plot_legend(biopsy_samples)$plot

plot_libsize(biopsy_samples)$plot

plot_nonzero(biopsy_samples)$plot

biopsy_prepost <- plot_libsize_prepost(biopsy_samples)
biopsy_prepost$count_plot

biopsy_prepost$lowgene_plot

## Warning: Using alpha for a discrete variable is not advised.

## Minimum number of biopsy genes: ~ 14,000

plot_libsize(eosinophil_samples)$plot

plot_nonzero(eosinophil_samples)$plot

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

eosinophil_prepost <- plot_libsize_prepost(eosinophil_samples)
eosinophil_prepost$count_plot

eosinophil_prepost$lowgene_plot

## Warning: Using alpha for a discrete variable is not advised.

## Minimum number of eosinophil genes: ~ 13,500

plot_libsize(monocyte_samples)$plot

plot_nonzero(monocyte_samples)$plot

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

monocyte_prepost <- plot_libsize_prepost(monocyte_samples)
monocyte_prepost$count_plot

monocyte_prepost$lowgene_plot

## Warning: Using alpha for a discrete variable is not advised.

## Minimum number of monocyte genes: ~ 7,500 before setting the minimum.

plot_libsize(neutrophil_samples)$plot

plot_nonzero(neutrophil_samples)$plot

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

neutrophil_prepost <- plot_libsize_prepost(neutrophil_samples)
neutrophil_prepost$count_plot

neutrophil_prepost$lowgene_plot

## Warning: Using alpha for a discrete variable is not advised.

## Minimum number of neutrophil genes: ~ 10,000 before setting minimum coverage.

The above block just repeats the same two plots on a per-celltype basis: the number of reads observed / sample and a plot of observed genes with respect to coverage. I made some comments with my observations about the number of genes.

6.2 Global PCA

Now that those ‘global’ metrics are out of the way, lets look at some global metrics of the data following normalization; the most likely plots are of course PCA but also a couple of heatmaps.

6.2.1 Figure 1

In the google doc TMRC3_Aug18_2021, there is an example of an image for the first figure:

“Transcriptomic profiles of primary innate cells of CL patients show unique transcriptional signatures - Remove PBMCs and M0, maybe biopsies as well (but Remove WT samples)”

While we were talking in a meeting however, it sounded like there was some desire to keep all cell types. Therefore the following block will have one image with everything and one following the above.

type_valid <- set_expt_conditions(hs_valid, fact="typeofcells") %>%
  set_expt_batches(fact="clinicaloutcome") %>%
  set_expt_colors(type_colors)

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

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

all_pca_nosize <- plot_pca(all_norm, plot_labels=FALSE)
all_pca_nosize$plot

write.csv(all_pca$table, file="coords/hs_donor_pca_coords.csv")
all_cf_norm <- set_expt_batches(all_norm,
                                fact="visitnumber")
all_cf_corheat <- plot_corheat(all_cf_norm, plot_title="Heirarchical clustering:
         cell types")

dev <- pp(file=glue("images/tmrc3_corheat_cf-v{ver}.png"))
all_cf_corheat$plot
closed <- dev.off()
all_cf_corheat$plot

all_cf_disheat <- plot_disheat(all_cf_norm, plot_title="Heirarchical clustering:
         cell types")
dev <- pp(file=glue("images/tmrc3_disheat_cf-v{ver}.png"))
all_cf_disheat$plot
closed <- dev.off()
all_cf_disheat$plot

6.3 Figure 1B: Transcriptomic profiles of primary innate

The biggest caveat for this is to ensure that there are no Wellcome Trust samples.

A potential figure legend for the following images might include:

The observed counts per gene for all of the clinical samples were filtered, log transformed, cpm converted, and quantile normalized. The colors were defined by cell types and shapes by patient visit. When the first two principle components were plotted, clustering was observed by cell type. The biopsy samples were significantly different from the innate immune cell types.

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

## Removing 5601 low-count genes (14322 remaining).

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

fig1v2_pca <- plot_pca(fig1v2_norm, cis=FALSE)

## plot labels was not set and there are more than 100 samples, disabling it.

dev <- pp(file=glue("images/tmrc3_fig1v2.png"))
fig1v2_pca$plot
closed <- dev.off()
fig1v2_pca$plot

fig1v3_samples <- subset_expt(type_valid, subset="condition!='biopsy'")

## subset_expt(): There were 192, now there are 173 samples.

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

## Removing 7749 low-count genes (12174 remaining).

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

fig1v3_pca <- plot_pca(fig1v3_norm, cis=FALSE)

## plot labels was not set and there are more than 100 samples, disabling it.

dev <- pp(file="images/tmrc3_fig1v3.png")
fig1v3_pca$plot
closed <- dev.off()
fig1v3_pca$plot

Continue looking, but switch the conditions/colors so that the clinical outcome becomes the focus and get rid of a few samples which are not actually a part of the TMRC3 focus (e.g. the PBMC and macrophage samples, which are all from the Wellcome Trust).

6.3.1 Clinically relevant samples

Included in this group will be the samples from patients who were lost.

6.3.1.1 Remove the lost samples

In our 20220218 meeting, it was decided that we would focus explicitly on the cure and fail samples, ignoring lost/NA/null samples.

Thus I am adding explicit filters right at the top to exclude them.

6.3.1.2 Samples by clinic

hs_clinical <- hs_valid %>%
  set_expt_conditions(fact="clinicaloutcome") %>%
  set_expt_batches(fact="typeofcells") %>%
  set_expt_colors(cf_colors)

hs_clinical_nobiop <- subset_expt(hs_clinical, subset="typeofcells!='biopsy'")

## subset_expt(): There were 192, now there are 173 samples.

6.3.1.3 View biopsies

clinic_biopsy <- hs_valid %>%
  set_expt_conditions(fact="clinic") %>%
  set_expt_batches(fact="clinicaloutcome") %>%
  subset_expt(subset="typeofcells=='biopsy'")

## subset_expt(): There were 192, now there are 19 samples.

clinic_cf <- paste0(pData(clinic_biopsy)$condition, "_",
                    pData(clinic_biopsy)$batch)
table(clinic_cf)

## clinic_cf
##      Cali_cure    Tumaco_cure Tumaco_failure 
##              4             10              5

clinic_biopsy <- set_expt_conditions(clinic_biopsy, fact=clinic_cf) %>%
  set_expt_batches(fact="visitnumber")

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

## Removing 6278 low-count genes (13645 remaining).

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

clinic_biopsy_pca <- plot_pca(clinic_biopsy_norm, plot_labels=FALSE)
dev <- pp(file="images/biopsy_place.png")
clinic_biopsy_pca$plot
closed <- dev.off()
clinic_biopsy_pca$plot

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

## Removing 6278 low-count genes (13645 remaining).

## Setting 320 low elements to zero.

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

clinic_biopsy_nb_pca <- plot_pca(clinic_biopsy_nb, plot_labels=FALSE)
dev <- pp(file="images/biopsy_place_nb.png")
clinic_biopsy_nb_pca$plot
closed <- dev.off()
clinic_biopsy_nb_pca$plot

6.3.1.4 Clinic Eosinophils

clinic_eosinophil <- hs_clinical %>%
  set_expt_conditions(fact="clinic") %>%
  set_expt_batches(fact="clinicaloutcome") %>%
  subset_expt(subset="typeofcells=='eosinophils'")

## subset_expt(): There were 192, now there are 42 samples.

clinic_cf <- paste0(pData(clinic_eosinophil)$condition, "_",
                    pData(clinic_eosinophil)$batch)
table(clinic_cf)

## clinic_cf
##      Cali_cure    Tumaco_cure Tumaco_failure 
##             15             18              9

clinic_eosinophil <- set_expt_conditions(clinic_eosinophil, fact=clinic_cf) %>%
    set_expt_batches(fact="visitnumber")

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

## Removing 9049 low-count genes (10874 remaining).

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

clinic_eosinophil_pca <- plot_pca(clinic_eosinophil_norm, plot_labels=FALSE)
dev <- pp(file="images/eosinophil_place.png")
clinic_eosinophil_pca$plot
closed <- dev.off()
clinic_eosinophil_pca$plot

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

## Removing 9049 low-count genes (10874 remaining).

## Setting 1033 low elements to zero.

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

clinic_eosinophil_nb_pca <- plot_pca(clinic_eosinophil_nb, plot_labels=FALSE)
dev <- pp(file="images/eosinophil_place_nb.png")
clinic_eosinophil_nb_pca$plot
closed <- dev.off()
clinic_eosinophil_nb_pca$plot

6.3.1.5 Clinic Monocyte

clinic_monocyte <- hs_clinical %>%
  set_expt_conditions(fact="clinic") %>%
  set_expt_batches(fact="clinicaloutcome") %>%
  subset_expt(subset="typeofcells=='monocytes'")

## subset_expt(): There were 192, now there are 66 samples.

clinic_cf <- paste0(pData(clinic_monocyte)$condition, "_",
                    pData(clinic_monocyte)$batch)
table(clinic_cf)

## clinic_cf
##      Cali_cure   Cali_failure    Tumaco_cure Tumaco_failure 
##             18              3             25             20

clinic_monocyte <- set_expt_conditions(clinic_monocyte, fact=clinic_cf) %>%
    set_expt_batches(fact="visitnumber")

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

## Removing 8797 low-count genes (11126 remaining).

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

clinic_monocyte_pca <- plot_pca(clinic_monocyte_norm, plot_labels=FALSE)
dev <- pp(file="images/monocytes_place.png")
clinic_monocyte_pca$plot
closed <- dev.off()
clinic_monocyte_pca$plot

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

## Removing 8797 low-count genes (11126 remaining).

## Setting 1626 low elements to zero.

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

clinic_monocyte_nb_pca <- plot_pca(clinic_monocyte_nb, plot_labels=FALSE)
dev <- pp(file="images/monocytes_place_nb.png")
clinic_monocyte_nb_pca$plot
closed <- dev.off()
clinic_monocyte_nb_pca$plot

6.4 Clinic Neutrophils

clinic_neutrophil <- hs_clinical %>%
  set_expt_conditions(fact="clinic") %>%
  set_expt_batches(fact="clinicaloutcome") %>%
  subset_expt(subset="typeofcells=='neutrophils'")

## subset_expt(): There were 192, now there are 65 samples.

clinic_cf <- paste0(pData(clinic_neutrophil)$condition, "_",
                    pData(clinic_neutrophil)$batch)
table(clinic_cf)

## clinic_cf
##      Cali_cure   Cali_failure    Tumaco_cure Tumaco_failure 
##             18              3             24             20

clinic_neutrophil <- set_expt_conditions(clinic_neutrophil, fact=clinic_cf) %>%
    set_expt_batches(fact="visitnumber")

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

## Removing 10634 low-count genes (9289 remaining).

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

clinic_neutrophil_pca <- plot_pca(clinic_neutrophil_norm, plot_labels=FALSE)
dev <- pp(file="images/neutrophil_place.png")
clinic_neutrophil_pca$plot
closed <- dev.off()
clinic_neutrophil_pca$plot

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

## Removing 10634 low-count genes (9289 remaining).

## Setting 1462 low elements to zero.

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

clinic_neutrophil_nb_pca <- plot_pca(clinic_neutrophil_nb, plot_labels=FALSE)
dev <- pp(file="images/neutrophil_place_nb.png")
clinic_neutrophil_nb_pca$plot
closed <- dev.off()
clinic_neutrophil_nb_pca$plot

6.4.1 Only Tumaco samples

A majority of samples came from tumaco, let us see what happens if we limit the data to them.

Let us therefore consider a couple of different view of this portion of the data.

tumaco_clinical <- hs_clinical %>%
  subset_expt(subset="clinic=='Tumaco'")

## subset_expt(): There were 192, now there are 131 samples.

tumaco_monocyte <- tumaco_clinical %>%
  subset_expt(subset="typeofcells=='monocytes'") %>%
  set_expt_conditions(fact="clinicaloutcome") %>%
  set_expt_batches(fact="visitnumber") %>%
  set_expt_samplenames("tubelabelorigin")

## subset_expt(): There were 131, now there are 45 samples.

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

## Removing 9025 low-count genes (10898 remaining).

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

tumaco_monocyte_pca <- plot_pca(tumaco_monocyte_norm, plot_labels=FALSE)
dev <- pp(file="images/monocytes_oneplace_norm.png")
tumaco_monocyte_pca$plot
closed <- dev.off()
tumaco_monocyte_pca$plot

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

## Removing 9025 low-count genes (10898 remaining).

## Setting 799 low elements to zero.

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

tumaco_monocyte_nb_pca <- plot_pca(tumaco_monocyte_nb, plot_labels=FALSE)
dev <- pp(file="images/monocytes_oneplace_norm_sva.png")
tumaco_monocyte_nb_pca$plot
closed <- dev.off()
tumaco_monocyte_nb_pca$plot

6.5 Compare clinics

Let us move a direct comparison of the two clinics right to the top.

hs_clinic <- hs_valid %>%
  set_expt_conditions(fact="clinic") %>%
  set_expt_batches(fact="typeofcells")

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

## Removing 5601 low-count genes (14322 remaining).

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

hs_clinic_pca <- plot_pca(hs_clinic_norm)

## plot labels was not set and there are more than 100 samples, disabling it.

hs_clinic_pca$plot

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

## Removing 5601 low-count genes (14322 remaining).

## Setting 33229 low elements to zero.

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

hs_clinic_nb_pca <- plot_pca(hs_clinic_nb)

## plot labels was not set and there are more than 100 samples, disabling it.

hs_clinic_nb_pca$plot

6.5.1 Compare clinics, all samples

clinic_contrasts <- list(
    "clinics" = c("Cali", "Tumaco"))
hs_clinic_de <- all_pairwise(hs_clinic, model_batch="svaseq", filter=TRUE)

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

## Setting 33229 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

hs_clinic_table <- combine_de_tables(
    hs_clinic_de, keepers=clinic_contrasts,
    excel=glue::glue("excel/hs_clinic_table-v{ver}.xlsx"))
hs_clinic_sig <- extract_significant_genes(
    hs_clinic_table,
    excel=glue::glue("excel/hs_clinic_sig-v{ver}.xlsx"))

6.5.2 Clinic biopsy comparison

clinic_biopsy_de <- all_pairwise(clinic_biopsy, model_batch="svaseq", filter=TRUE)

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

## Setting 320 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

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

## Deleting the file excel/clinic_biopsy_table-v202205.xlsx before writing the tables.

clinic_biopsy_sig <- extract_significant_genes(
    clinic_biopsy_table,
    excel=glue::glue("excel/clinic_biopsy_sig-v{ver}.xlsx"))

## Deleting the file excel/clinic_biopsy_sig-v202205.xlsx before writing the tables.

6.5.3 Clinic eosinophil comparison

clinic_eosinophil_de <- all_pairwise(clinic_eosinophil, model_batch=FALSE, filter=TRUE)

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

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

## Deleting the file excel/clinic_eosinophil_table-v202205.xlsx before writing the tables.

clinic_eosinophil_sig <- extract_significant_genes(
    clinic_eosinophil_table,
    excel=glue::glue("excel/clinic_eosinophil_sig-v{ver}.xlsx"))

## Deleting the file excel/clinic_eosinophil_sig-v202205.xlsx before writing the tables.

6.5.4 Clinic monocyte comparison

clinic_monocyte_de <- all_pairwise(clinic_monocyte, model_batch=FALSE, filter=TRUE)

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

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

## Deleting the file excel/clinic_monocyte_table-v202205.xlsx before writing the tables.

clinic_monocyte_sig <- extract_significant_genes(
    clinic_monocyte_table,
    excel=glue::glue("excel/clinic_monocyte_sig-v{ver}.xlsx"))

## Deleting the file excel/clinic_monocyte_sig-v202205.xlsx before writing the tables.

6.5.5 Clinic neutrophil comparison

clinic_neutrophil_de <- all_pairwise(clinic_neutrophil, model_batch=FALSE, filter=TRUE)

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

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

## Deleting the file excel/clinic_neutrophil_table-v202205.xlsx before writing the tables.

clinic_neutrophil_sig <- extract_significant_genes(
    clinic_neutrophil_table,
    excel=glue::glue("excel/clinic_neutrophil_sig-v{ver}.xlsx"))

## Deleting the file excel/clinic_neutrophil_sig-v202205.xlsx before writing the tables.

6.5.6 Compare clinic genes

clinic_sigenes_up <- rownames(hs_clinic_sig[["deseq"]][["ups"]][[1]])
clinic_sigenes_down <- rownames(hs_clinic_sig[["deseq"]][["downs"]][[1]])
clinic_sigenes <- c(clinic_sigenes_up, clinic_sigenes_down)

contrast <- "Tumacocure_vs_Calicure"
clinic_biopsy_sigenes <- c(rownames(clinic_biopsy_sig[["deseq"]][["ups"]][[contrast]]),
                           rownames(clinic_biopsy_sig[["deseq"]][["downs"]][[contrast]]))
clinic_eosinophil_sigenes_up <- rownames(clinic_eosinophil_sig[["deseq"]][["ups"]][[contrast]])
clinic_eosinophil_sigenes_down <- rownames(clinic_eosinophil_sig[["deseq"]][["downs"]][[contrast]])
clinic_monocyte_sigenes_up <- rownames(clinic_monocyte_sig[["deseq"]][["ups"]][[contrast]])
clinic_monocyte_sigenes_down <- rownames(clinic_monocyte_sig[["deseq"]][["downs"]][[contrast]])
clinic_neutrophil_sigenes_up <- rownames(clinic_neutrophil_sig[["deseq"]][["ups"]][[contrast]])
clinic_neutrophil_sigenes_down <- rownames(clinic_neutrophil_sig[["deseq"]][["downs"]][[contrast]])

clinic_eosinophil_sigenes <- c(clinic_eosinophil_sigenes_up,
                               clinic_eosinophil_sigenes_down)
clinic_monocyte_sigenes <- c(clinic_monocyte_sigenes_up,
                             clinic_monocyte_sigenes_down)
clinic_neutrophil_sigenes <- c(clinic_neutrophil_sigenes_up,
                               clinic_neutrophil_sigenes_down)

6.5.7 Place gprofiler

clinic_gp <- simple_gprofiler(clinic_sigenes)

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

## GO search found 167 hits.

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

## KEGG search found 12 hits.

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

## REAC search found 13 hits.

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

## MI search found 0 hits.

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

## TF search found 17 hits.

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

## CORUM search found 1 hits.

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

## HP search found 17 hits.

clinic_gp$pvalue_plots$kegg_plot_over

clinic_gp$pvalue_plots$reactome_plot_over

clinic_gp <- simple_gprofiler(clinic_sigenes)

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

## GO search found 167 hits.

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

## KEGG search found 12 hits.

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

## REAC search found 13 hits.

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

## MI search found 0 hits.

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

## TF search found 17 hits.

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

## CORUM search found 1 hits.

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

## HP search found 17 hits.

clinic_gp$pvalue_plots$kegg_plot_over

clinic_gp$pvalue_plots$reactome_plot_over

clinic_gp <- simple_gprofiler(clinic_sigenes)

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

## GO search found 167 hits.

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

## KEGG search found 12 hits.

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

## REAC search found 13 hits.

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

## MI search found 0 hits.

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

## TF search found 17 hits.

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

## CORUM search found 1 hits.

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

## HP search found 17 hits.

clinic_gp$pvalue_plots$kegg_plot_over

clinic_gp$pvalue_plots$reactome_plot_over

clinic_eosinophil_gp <- simple_gprofiler(clinic_eosinophil_sigenes)

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

## GO search found 329 hits.

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

## KEGG search found 32 hits.

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

## REAC search found 17 hits.

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

## MI search found 0 hits.

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

## TF search found 371 hits.

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

## CORUM search found 10 hits.

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

## HP search found 2 hits.

clinic_eosinophil_gp$pvalue_plots$kegg_plot_over

clinic_eosinophil_gp$pvalue_plots$reactome_plot_over

clinic_eosinophil_up_gp <- simple_gprofiler(clinic_eosinophil_sigenes_up)

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

## GO search found 234 hits.

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

## KEGG search found 28 hits.

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

## REAC search found 9 hits.

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

## MI search found 0 hits.

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

## TF search found 341 hits.

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

## CORUM search found 7 hits.

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

## HP search found 2 hits.

clinic_eosinophil_up_gp$pvalue_plots$kegg_plot_over

clinic_eosinophil_up_gp$pvalue_plots$reactome_plot_over

clinic_eosinophil_down_gp <- simple_gprofiler(clinic_eosinophil_sigenes_down)

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

## GO search found 184 hits.

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

## KEGG search found 4 hits.

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

## REAC search found 11 hits.

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

## MI search found 0 hits.

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

## TF search found 39 hits.

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

## CORUM search found 3 hits.

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

## HP search found 1 hits.

clinic_eosinophil_down_gp$pvalue_plots$kegg_plot_over

clinic_eosinophil_down_gp$pvalue_plots$reactome_plot_over

clinic_monocyte_gp <- simple_gprofiler(clinic_monocyte_sigenes)

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

## GO search found 417 hits.

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

## KEGG search found 21 hits.

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

## REAC search found 13 hits.

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

## MI search found 0 hits.

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

## TF search found 403 hits.

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

## CORUM search found 7 hits.

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

## HP search found 22 hits.

clinic_monocyte_gp$pvalue_plots$kegg_plot_over

clinic_monocyte_gp$pvalue_plots$reactome_plot_over

clinic_monocyte_up_gp <- simple_gprofiler(clinic_monocyte_sigenes_up)

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

## GO search found 382 hits.

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

## KEGG search found 19 hits.

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

## REAC search found 10 hits.

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

## MI search found 0 hits.

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

## TF search found 385 hits.

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

## CORUM search found 6 hits.

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

## HP search found 13 hits.

clinic_monocyte_up_gp$pvalue_plots$kegg_plot_over

clinic_monocyte_up_gp$pvalue_plots$reactome_plot_over

clinic_monocyte_down_gp <- simple_gprofiler(clinic_monocyte_sigenes_down)

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

## GO search found 60 hits.

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

## KEGG search found 3 hits.

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

## REAC search found 4 hits.

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

## MI search found 0 hits.

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

## TF search found 241 hits.

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

## CORUM search found 0 hits.

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

## HP search found 10 hits.

clinic_monocyte_down_gp$pvalue_plots$kegg_plot_over

clinic_monocyte_down_gp$pvalue_plots$reactome_plot_over

clinic_neutrophil_gp <- simple_gprofiler(clinic_neutrophil_sigenes)

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

## GO search found 153 hits.

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

## KEGG search found 10 hits.

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

## REAC search found 20 hits.

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

## MI search found 0 hits.

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

## TF search found 416 hits.

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

## CORUM search found 1 hits.

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

## HP search found 1 hits.

clinic_neutrophil_gp$pvalue_plots$kegg_plot_over

clinic_neutrophil_gp$pvalue_plots$reactome_plot_over

clinic_neutrophil_up_gp <- simple_gprofiler(clinic_neutrophil_sigenes_up)

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

## GO search found 134 hits.

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

## KEGG search found 5 hits.

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

## REAC search found 3 hits.

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

## MI search found 0 hits.

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

## TF search found 388 hits.

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

## CORUM search found 0 hits.

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

## HP search found 0 hits.

clinic_neutrophil_up_gp$pvalue_plots$kegg_plot_over

clinic_neutrophil_up_gp$pvalue_plots$reactome_plot_over

clinic_neutrophil_down_gp <- simple_gprofiler(clinic_neutrophil_sigenes_down)

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

## GO search found 62 hits.

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

## KEGG search found 12 hits.

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

## REAC search found 77 hits.

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

## MI search found 0 hits.

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

## TF search found 21 hits.

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

## CORUM search found 9 hits.

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

## HP search found 7 hits.

clinic_neutrophil_down_gp$pvalue_plots$kegg_plot_over

clinic_neutrophil_down_gp$pvalue_plots$reactome_plot_over

6.5.8 Plot the clinical samples

Let us recolor the same plot by cure/fail followed by a concatenation of the cell type and cure/fail. In the following block, the clinical samples are plotted once with the most common normalization (log,cpm,quant,filtered) method followed by a plot of the data using svaseq adjusted values and without quantile normalization.

Thus, the following block switches the colors of the groups to the clinical state (cure/fail) and shapes by cell type.

hs_clinical_norm <- sm(normalize_expt(hs_clinical, filter="simple", transform="log2",
                                      norm="quant", convert="cpm"))
clinical_pca <- plot_pca(hs_clinical_norm, plot_labels=FALSE,
                         cis=NULL,
                         plot_title="PCA - clinical samples")
dev <- pp(file=glue("images/all_clinical_nobatch_pca-v{ver}.png"), height=8, width=16)
clinical_pca$plot
closed <- dev.off()
clinical_pca$plot

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

## Removing 1865 low-count genes (18058 remaining).

## Setting 158794 low elements to zero.

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

hs_clinical_nb_pca <- plot_pca(hs_clinical_nb)

## plot labels was not set and there are more than 100 samples, disabling it.

dev <- pp(file=glue("images/all_clinical_svaseqbatch_pca-v{ver}.png"), height=6, width=8)
hs_clinical_nb_pca$plot
closed <- dev.off()
hs_clinical_nb_pca$plot

clinical_pca_info <- pca_information(
    hs_clinical_norm, plot_pcas=TRUE, num_components = 30,
    expt_factors=c("visitnumber", "typeofcells", "clinicaloutcome",
                   "clinic"))

## plot labels was not set and there are more than 100 samples, disabling it.

dev <- pp(file="images/clinical_samples_neglogp_pcs.png")
clinical_pca_info$anova_neglogp_heatmap
closed <- dev.off()
clinical_pca_info$anova_neglogp_heatmap

clinical_pca_info$pca_plots$PC20_PC29

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

clinical_scores <- pca_highscores(hs_clinical_norm)
clinical_scores[["highest"]][,"Comp.20"]

##  [1] "14.03:ENSG00000266302" "9.344:ENSG00000163993" "6.935:ENSG00000129824"
##  [4] "6.362:ENSG00000067048" "5.933:ENSG00000171860" "5.535:ENSG00000196526"
##  [7] "5.304:ENSG00000198692" "5.088:ENSG00000099725" "4.65:ENSG00000185897" 
## [10] "4.499:ENSG00000106565" "4.471:ENSG00000144681" "4.422:ENSG00000012817"
## [13] "4.381:ENSG00000176834" "4.354:ENSG00000118432" "4.285:ENSG00000073464"
## [16] "4.145:ENSG00000129295" "4.134:ENSG00000178538" "3.972:ENSG00000198178"
## [19] "3.944:ENSG00000134460" "3.93:ENSG00000152766"

clinical_scores[["highest"]][,"Comp.27"]

##  [1] "13.72:ENSG00000244734" "11.84:ENSG00000188536" "9.716:ENSG00000206172"
##  [4] "5.305:ENSG00000266302" "4.884:ENSG00000183570" "4.333:ENSG00000277632"
##  [7] "3.539:ENSG00000130766" "3.435:ENSG00000198848" "3.419:ENSG00000158578"
## [10] "3.401:ENSG00000136732" "3.168:ENSG00000134824" "3.142:ENSG00000101220"
## [13] "3.069:ENSG00000237541" "3.023:ENSG00000122877" "2.958:ENSG00000120738"
## [16] "2.957:ENSG00000123358" "2.948:ENSG00000223609" "2.937:ENSG00000142583"
## [19] "2.919:ENSG00000133048" "2.913:ENSG00000153283"

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

## Removing 5601 low-count genes (14322 remaining).

## Setting 206667 low elements to zero.

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

first_info <- pca_information(
    first, plot_pcas=TRUE, num_components = 30,
    expt_factors=c("visitnumber", "typeofcells", "clinicaloutcome",
                   "clinic"))

## plot labels was not set and there are more than 100 samples, disabling it.

first_info$anova_neglogp_heatmap

first_info$pca_plots[["PC1_PC2"]]

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

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

second <- normalize_expt(hs_clinical, transform="log2", convert="cpm",
                         filter = TRUE, batch="svaseq", surrogates=2) %>%
  set_expt_batches(fact="clinic")

## Removing 5601 low-count genes (14322 remaining).

## Setting 33691 low elements to zero.

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

second_info <- pca_information(
    second, plot_pcas=TRUE, num_components = 30,
    expt_factors=c("visitnumber", "typeofcells", "clinicaloutcome",
                   "clinic"))

## plot labels was not set and there are more than 100 samples, disabling it.

second_info$anova_neglogp_heatmap

third <- normalize_expt(hs_clinical, transform="log2", convert="cpm",
                        filter = TRUE, batch="svaseq", surrogates=3) %>%
  set_expt_batches(fact="clinic")

## Removing 5601 low-count genes (14322 remaining).

## Setting 28785 low elements to zero.

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

third_info <- pca_information(
    third, plot_pcas=TRUE, num_components = 30,
    expt_factors=c("visitnumber", "typeofcells", "clinicaloutcome",
                   "clinic"))

## plot labels was not set and there are more than 100 samples, disabling it.

third_info$anova_neglogp_heatmap

fourth <- normalize_expt(hs_clinical, transform="log2", convert="cpm",
                         filter = TRUE, batch="svaseq", surrogates=4) %>%
  set_expt_batches(fact="clinic")

## Removing 5601 low-count genes (14322 remaining).

## Setting 27785 low elements to zero.

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

fourth_info <- pca_information(
    fourth, plot_pcas=TRUE, num_components = 30,
    expt_factors=c("visitnumber", "typeofcells", "clinicaloutcome",
                   "clinic"))

## plot labels was not set and there are more than 100 samples, disabling it.

fourth_info$anova_neglogp_heatmap

fourth_info[["pca_plots"]][["PC1_PC2"]]

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

fifth <- normalize_expt(hs_clinical, transform="log2", convert="cpm",
                        filter = TRUE, batch="svaseq", surrogates=5) %>%
  set_expt_batches(fact="clinic")

## Removing 5601 low-count genes (14322 remaining).

## Setting 28885 low elements to zero.

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

fifth_info <- pca_information(
    fifth, plot_pcas=TRUE, num_components = 30,
    expt_factors=c("visitnumber", "typeofcells", "clinicaloutcome",
                   "clinic"))

## plot labels was not set and there are more than 100 samples, disabling it.

fifth_info$anova_neglogp_heatmap

fifth_info[["pca_plots"]][["PC1_PC12"]]

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

sixth <- normalize_expt(hs_clinical, transform="log2", convert="cpm",
                        filter = TRUE, batch="svaseq", surrogates=6) %>%
  set_expt_batches(fact="clinic")

## Removing 5601 low-count genes (14322 remaining).

## Setting 25457 low elements to zero.

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

sixth_info <- pca_information(
    sixth, plot_pcas=TRUE, num_components = 30,
    expt_factors=c("visitnumber", "typeofcells", "clinicaloutcome",
                   "clinic"))

## plot labels was not set and there are more than 100 samples, disabling it.

sixth_info$anova_neglogp_heatmap

seventh <- normalize_expt(hs_clinical, transform="log2", convert="cpm",
                          filter = TRUE, batch="svaseq", surrogates=7) %>%
  set_expt_batches(fact="clinic")

## Removing 5601 low-count genes (14322 remaining).

## Setting 25447 low elements to zero.

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

seventh_info <- pca_information(
    seventh, plot_pcas=TRUE, num_components = 30,
    expt_factors=c("visitnumber", "typeofcells", "clinicaloutcome",
                   "clinic"))

## plot labels was not set and there are more than 100 samples, disabling it.

seventh_info$anova_neglogp_heatmap

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

## Removing 5601 low-count genes (14322 remaining).

## Setting 25130 low elements to zero.

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

eighth_info <- pca_information(
    eighth, plot_pcas=TRUE, num_components = 30,
    expt_factors=c("visitnumber", "typeofcells", "clinicaloutcome",
                   "clinic"))

## plot labels was not set and there are more than 100 samples, disabling it.

eighth_info$anova_neglogp_heatmap

6.5.9 Concatenate cure/fail and cell type

In the following block the experimental condition was reset to the concatenation of clinical outcome and type of cells. There are an insufficient number of biopsy samples for them to be useful in this visualization, so they are ignored.

desired_levels <- c("cure_biopsy", "failure_biopsy", "cure_eosinophils", "failure_eosinophils",
                    "cure_monocytes", "failure_monocytes", "cure_neutrophils", "failure_neutrophils")
new_fact <- factor(
    paste0(pData(hs_clinical)[["condition"]], "_",
           pData(hs_clinical)[["batch"]]),
    levels=desired_levels)

hs_clinical_concat <- set_expt_conditions(hs_clinical, fact = new_fact) %>%
  set_expt_batches(fact = "visitnumber") %>%
  set_expt_colors(cf_type_colors) %>%
  subset_expt(subset="typeofcells!='biopsy'")

## subset_expt(): There were 192, now there are 173 samples.

## Try to ensure that the levels stay in the order I want
meta <- pData(hs_clinical_concat) %>%
  mutate(condition = fct_relevel(condition, desired_levels))

## Warning: Unknown levels in `f`: cure_biopsy, failure_biopsy

pData(hs_expt) <- meta

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

## Removing 7749 low-count genes (12174 remaining).

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

clinical_concat_norm_pca <- plot_pca(clinical_concat_norm)

## plot labels was not set and there are more than 100 samples, disabling it.

dev <- pp(file=glue("images/clinical_concatenated_normalized_pca-v{ver}.png"), height=6, width=10)
clinical_concat_norm_pca$plot
closed <- dev.off()
clinical_concat_norm_pca$plot

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

## Removing 7749 low-count genes (12174 remaining).

## Setting 37479 low elements to zero.

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

clinical_concat_nb_pca <- plot_pca(clinical_concat_nb)

## plot labels was not set and there are more than 100 samples, disabling it.

dev <- pp(file=glue("images/clinical_concatenated_svaseqbatch_pca-v{ver}.png"), height=6, width=12)
clinical_concat_nb_pca$plot
closed <- dev.off()
clinical_concat_nb_pca$plot

tumaco_concat <- subset_expt(hs_clinical_concat, subset="clinic=='Tumaco'")

## subset_expt(): There were 173, now there are 116 samples.

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

## Removing 7985 low-count genes (11938 remaining).

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

tumaco_concat_norm_pca <- plot_pca(tumaco_concat_norm)

## plot labels was not set and there are more than 100 samples, disabling it.

dev <- pp(file=glue("images/tumaco_concatenated_normalized_pca-v{ver}.png"), height=6, width=10)
tumaco_concat_norm_pca$plot
closed <- dev.off()
tumaco_concat_norm_pca$plot

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

## Removing 7985 low-count genes (11938 remaining).

## Setting 32996 low elements to zero.

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

tumaco_concat_nb_pca <- plot_pca(tumaco_concat_nb)

## plot labels was not set and there are more than 100 samples, disabling it.

dev <- pp(file=glue("images/clinical_concatenated_svaseqbatch_pca-v{ver}.png"), height=6, width=12)
tumaco_concat_nb_pca$plot
closed <- dev.off()
tumaco_concat_nb_pca$plot

library(ggplot2)
wanted_legend_order <- c("cure_eosinophils","cure_monocytes","cure_neutrophils",
                         "failure_eosinophils","failure_monocytes","failure_neutrophils")
test <- plot_pca(clinical_concat_norm)

## plot labels was not set and there are more than 100 samples, disabling it.

test$plot

test$plot +
  scale_color_discrete(breaks=wanted_legend_order) +
  guides(color=guide_legend(ncol=2))

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

test2 <- plot_legend(clinical_concat_norm)

## plot labels was not set and there are more than 100 samples, disabling it.

6.5.10 Samples separated by visit

Separate the samples by visit in order to more easily see what patterns emerge across cell type and clinical outcome.

The most likely data set for this is hs_clinical.

6.5.10.1 All visits together

Now let us shift the view slightly to focus on changes observed over time.

visit_expt <- set_expt_conditions(hs_clinical, fact = "visitnumber") %>%
  set_expt_batches(fact = "clinicaloutcome") %>%
  subset_expt(subset="typeofcells!='biopsy'")

## subset_expt(): There were 192, now there are 173 samples.

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

## Removing 7749 low-count genes (12174 remaining).

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

plot_pca(visit_norm)$plot

## plot labels was not set and there are more than 100 samples, disabling it.

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

## Removing 7749 low-count genes (12174 remaining).

## Setting 18454 low elements to zero.

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

visit_nb_pca <- plot_pca(visit_nb)

## plot labels was not set and there are more than 100 samples, disabling it.

dev <- pp(file=glue("images/visit_svaseqbatch_pca-v{ver}.png"), height=7, width=9)
visit_nb_pca$plot
closed <- dev.off()
visit_nb_pca$plot

6.6 Top 1000 genes

repeat with only the most expressed genes from one cell type.

rpkm_values <- normalize_expt(type_valid, transform="log2", convert="rpkm",
                              filter=TRUE, column="cds_length")

## Removing 5601 low-count genes (14322 remaining).

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

rpkm_median <- median_by_factor(rpkm_values)

## The factor biopsy has 19 rows.

## The factor eosinophils has 42 rows.

## The factor monocytes has 66 rows.

## The factor neutrophils has 65 rows.

rpkm_tumaco_monocyte_idx <- order(rpkm_median[["medians"]][["monocytes"]], decreasing=TRUE)
rpkm_tumaco_monocyte <- rpkm_median[["medians"]][rpkm_tumaco_monocyte_idx, ]
tumaco_monocyte_gene <- head(rownames(rpkm_tumaco_monocyte), n=1000)

6.6.0.1 All visits, separated by cell type

Now separate the visits by cell type and whether or not the drug used was the antimonial.

visit_monocyte <- subset_expt(visit_expt, subset = "typeofcells=='monocytes'")

## subset_expt(): There were 173, now there are 66 samples.

visit_neutrophil <- subset_expt(visit_expt, subset = "typeofcells=='neutrophils'")

## subset_expt(): There were 173, now there are 65 samples.

visit_eosinophil <- subset_expt(visit_expt, subset = "typeofcells=='eosinophils'")

## subset_expt(): There were 173, now there are 42 samples.

6.6.0.2 Visit 1 alone

v1_clinical <- subset_expt(visit_expt, subset = "visitnumber=='1'") %>%
  set_expt_conditions(fact = "clinicaloutcome") %>%
  set_expt_batches(fact = "typeofcells")

## subset_expt(): There were 173, now there are 72 samples.

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

## Removing 8111 low-count genes (11812 remaining).

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

plot_pca(v1_norm)$plot

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

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

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

## Removing 8111 low-count genes (11812 remaining).

## Setting 5286 low elements to zero.

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

plot_pca(v1_nb, plot_labels = FALSE)$plot

So, given the data on hand, there is slight ability to discern cure and fail among the visit 1 samples.

6.6.0.3 Visit 2 alone

v2_clinical <- subset_expt(visit_expt, subset="visitnumber=='2'") %>%
  set_expt_conditions(fact = "clinicaloutcome") %>%
  set_expt_batches(fact = "typeofcells")

## subset_expt(): There were 173, now there are 50 samples.

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

## Warning in normalize_expt(v2_clinical, transform = "log2", convert = "cpm", :
## Quantile normalization and sva do not always play well together.

## Removing 8184 low-count genes (11739 remaining).

## Setting 2960 low elements to zero.

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

plot_pca(v2_nb, plot_labels = FALSE)$plot

There might be some variance among the visit 2 samples which are cure/fail-ish.

6.6.0.4 Visit 3 alone

v3_clinical <- subset_expt(visit_expt, subset="visitnumber=='3'") %>%
  set_expt_conditions(fact = "clinicaloutcome") %>%
  set_expt_batches(fact = "typeofcells")

## subset_expt(): There were 173, now there are 51 samples.

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

## Warning in normalize_expt(v3_clinical, transform = "log2", convert = "cpm", :
## Quantile normalization and sva do not always play well together.

## Removing 8259 low-count genes (11664 remaining).

## Setting 2584 low elements to zero.

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

plot_pca(v3_nb, plot_labels = FALSE)$plot

In a similar fashion, it seems that the visit 3 samples are not particularly informative.

6.6.1 Samples separated by cell type

Separate the samples by cell type in order to more easily observe patterns with respect to visit and clinical outcome.

6.6.1.1 Biopsies

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

## Removing 6278 low-count genes (13645 remaining).

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

plot_pca(biopsy_norm, plot_labels = FALSE)$plot

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

## Removing 6278 low-count genes (13645 remaining).

## Setting 243 low elements to zero.

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

plot_pca(biopsy_nb, plot_labels = FALSE)$plot

6.6.1.2 Monocytes

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

## Removing 8797 low-count genes (11126 remaining).

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

monocyte_pca <- plot_pca(monocyte_norm, plot_labels = FALSE)
dev <- pp(file="images/monocytes_cf_norm_pca.png")
monocyte_pca$plot
closed <- dev.off()
monocyte_pca$plot

monocyte_disheat <- plot_disheat(monocyte_norm)
dev <- pp(file="images/monocytes_cf_norm_disheat.png")
monocyte_disheat$plot
closed <- dev.off()
monocyte_disheat$plot

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

## Removing 8797 low-count genes (11126 remaining).

## Setting 1403 low elements to zero.

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

monocyte_nb_pca <- plot_pca(monocyte_nb, plot_labels = FALSE)
dev <- pp(file="images/monocytes_cf_norm_sva_pca.png")
monocyte_nb_pca$plot
closed <- dev.off()
monocyte_nb_pca$plot

6.6.1.2.1 Monocytes, visit 1

monocyte_v1 <- subset_expt(monocyte_samples, subset = "visitnumber=='1'")

## subset_expt(): There were 66, now there are 29 samples.

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

## Removing 9053 low-count genes (10870 remaining).

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

monocyte_v1_pca <- plot_pca(monocyte_v1_norm, plot_labels = FALSE)
dev <- pp(file="images/monocytes_v1_cf_norm_pca.png")
monocyte_v1_pca$plot
closed <- dev.off()
monocyte_v1_pca$plot

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

## Removing 9053 low-count genes (10870 remaining).

## Setting 558 low elements to zero.

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

monocyte_v1_nb_pca <- plot_pca(monocyte_v1_nb, plot_labels = FALSE)
dev <- pp(file="images/monocytes_v1_cf_norm_sva_pca.png")
monocyte_v1_nb_pca$plot
closed <- dev.off()
monocyte_v1_nb_pca$plot

The monocytes, on the other hand, appear to have some real information lurking in them.

6.6.1.3 Neutrophils

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

## Removing 10634 low-count genes (9289 remaining).

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

plot_pca(neutrophil_norm, plot_labels = FALSE)$plot

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

## Removing 10634 low-count genes (9289 remaining).

## Setting 1703 low elements to zero.

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

plot_pca(neutrophil_nb, plot_labels = FALSE)$plot

This appears also to be the case for the neutrophil samples.

6.6.1.4 Eosinophils

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

## Removing 9049 low-count genes (10874 remaining).

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

plot_pca(eosinophil_norm, plot_labels = FALSE)$plot

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

## Removing 9049 low-count genes (10874 remaining).

## Setting 936 low elements to zero.

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

plot_pca(eosinophil_nb, plot_labels = FALSE)$plot

We have fewer eosinophil samples currently, but they appear to have some real differences.

6.6.1.5 Monocytes, Neutrophils, and Eosinophils

7 Contrasts and colors of interest

This might be a bit early to consider the contrasts, but I think we should consider this question immediately. The main thing we will be comparing is of course cure vs. fail; but we may also look at the visits and compare cell types.

cf_contrasts <- list(
    "fail_vs_cure" = c("failure", "cure"))
visit_contrasts <- list(
    "v2v1" = c("c2", "c1"),
    "v3v1" = c("c3", "c1"),
    "v3v2" = c("c3", "c2"))
type_contrasts <- list(
    "mono_biopsy" = c("monocytes", "biopsy"),
    "eosinophil_biopsy" = c("eosinophils", "biopsy"),
    "neutrophil_biopsy" = c("neutrophils", "biopsy"))
visit_cf_contrasts <- list(
    "v1fail_vs_cure" = c("v1failure", "v1cure"),
    "v2fail_vs_cure" = c("v2failure", "v2cure"),
    "v3fail_vs_cure" = c("v3failure", "v3cure"))

Sample IDs starting with 1: Cali IDs starting with 2: Tumaco

8 Differential expression analyses

The primary goal is to learn about cure vs. fail.

8.1 Cure/Fail, all samples

Now let us start performing the various differential expression analyses, starting with the set of all/most clinical samples.

cf_clinical_de <- all_pairwise(hs_clinical, model_batch = "svaseq", filter = TRUE)

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

## Setting 25457 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

cf_clinical_tables <- combine_de_tables(
    cf_clinical_de,
    keepers = cf_contrasts,
    excel = glue::glue("excel/cf_clinical_tables-v{ver}.xlsx"))

## Deleting the file excel/cf_clinical_tables-v202205.xlsx before writing the tables.

cf_clinical_sig <- extract_significant_genes(
    cf_clinical_tables,
    excel = glue::glue("excel/cf_clinical_sig-v{ver}.xlsx"))

## Deleting the file excel/cf_clinical_sig-v202205.xlsx before writing the tables.

8.1.1 By cell type

8.1.1.1 Cure/Fail, Biopsies

cf_biopsy_de <- all_pairwise(biopsy_samples, model_batch = "svaseq", filter = TRUE)

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

## Setting 243 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

cf_biopsy_tables <- combine_de_tables(
    cf_biopsy_de,
    keepers = cf_contrasts,
    excel = glue::glue("excel/cf_biopsy_tables-v{ver}.xlsx"))

## Deleting the file excel/cf_biopsy_tables-v202205.xlsx before writing the tables.

cf_biopsy_sig <- extract_significant_genes(
    cf_biopsy_tables,
    excel = glue::glue("excel/cf_biopsy_sig-v{ver}.xlsx"))

## Deleting the file excel/cf_biopsy_sig-v202205.xlsx before writing the tables.

8.1.1.2 Cure/Fail, Monocytes

cf_monocyte_sva_de <- all_pairwise(monocyte_samples, model_batch = "svaseq", filter = TRUE)

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

## Setting 1403 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

cf_monocyte_sva_tables <- combine_de_tables(
    cf_monocyte_sva_de,
    keepers = cf_contrasts,
    excel = glue::glue("excel/cf_monocyte_sva_tables-v{ver}.xlsx"))

## Deleting the file excel/cf_monocyte_sva_tables-v202205.xlsx before writing the tables.

cf_monocyte_sva_sig <- extract_significant_genes(
    cf_monocyte_sva_tables,
    excel = glue::glue("excel/cf_monocyte_sva_sig-v{ver}.xlsx"))

## Deleting the file excel/cf_monocyte_sva_sig-v202205.xlsx before writing the tables.

cf_monocyte_batchvisit_de <- all_pairwise(monocyte_samples, model_batch = TRUE, filter = TRUE)

## Using limma's removeBatchEffect to visualize with(out) batch inclusion.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

cf_monocyte_batchvisit_tables <- combine_de_tables(
    cf_monocyte_batchvisit_de,
    keepers = cf_contrasts,
    excel = glue::glue("excel/cf_monocyte_batchvisit_tables-v{ver}.xlsx"))

## Deleting the file excel/cf_monocyte_batchvisit_tables-v202205.xlsx before writing the tables.

cf_monocyte_batchvisit_sig <- extract_significant_genes(
    cf_monocyte_batchvisit_tables,
    excel = glue::glue("excel/cf_monocyte_batchvisit_sig-v{ver}.xlsx"))

## Deleting the file excel/cf_monocyte_batchvisit_sig-v202205.xlsx before writing the tables.

8.1.1.3 Cure/Fail, only the Tumaco samples

cf_tumaco_monocyte_sva_de <- all_pairwise(tumaco_monocyte, model_batch = "svaseq", filter = TRUE)

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

## Setting 799 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

cf_tumaco_monocyte_sva_tables <- combine_de_tables(
    cf_tumaco_monocyte_sva_de,
    keepers = cf_contrasts,
    excel = glue::glue("excel/cf_monocyte_tables_sva_tumaco-v{ver}.xlsx"))

## Deleting the file excel/cf_monocyte_tables_sva_tumaco-v202205.xlsx before writing the tables.

cf_tumaco_monocyte_sva_sig <- extract_significant_genes(
    cf_tumaco_monocyte_sva_tables,
    excel = glue::glue("excel/cf_monocyte_sva_tumaco_sig-v{ver}.xlsx"))

## Deleting the file excel/cf_monocyte_sva_tumaco_sig-v202205.xlsx before writing the tables.

cf_tumaco_monocyte_batchvisit_de <- all_pairwise(tumaco_monocyte, model_batch = TRUE, filter = TRUE)

## Using limma's removeBatchEffect to visualize with(out) batch inclusion.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

cf_tumaco_monocyte_batchvisit_tables <- combine_de_tables(
    cf_tumaco_monocyte_batchvisit_de,
    keepers = cf_contrasts,
    excel = glue::glue("excel/cf_monocyte_tables_batchvisit_tumaco-v{ver}.xlsx"))

## Deleting the file excel/cf_monocyte_tables_batchvisit_tumaco-v202205.xlsx before writing the tables.

cf_tumaco_monocyte_batchvisit_sig <- extract_significant_genes(
    cf_tumaco_monocyte_batchvisit_tables,
    excel = glue::glue("excel/cf_monocyte_batchvisit_tumaco_sig-v{ver}.xlsx"))

## Deleting the file excel/cf_monocyte_batchvisit_tumaco_sig-v202205.xlsx before writing the tables.

tumaco_sva_aucc <- calculate_aucc(cf_tumaco_monocyte_sva_tables[["data"]][["fail_vs_cure"]],
                            tbl2=cf_monocyte_sva_tables[["data"]][["fail_vs_cure"]],
                            py="deseq_adjp", ly="deseq_logfc")
tumaco_sva_aucc

## $aucc
## [1] 0.4024
## 
## $plot

shared_ids <- rownames(cf_tumaco_monocyte_sva_tables[["data"]][["fail_vs_cure"]]) %in%
  rownames(cf_monocyte_sva_tables[["data"]][["fail_vs_cure"]])
first <- cf_monocyte_sva_tables[["data"]][["fail_vs_cure"]][shared_ids, ]
second <- cf_tumaco_monocyte_sva_tables[["data"]][["fail_vs_cure"]][rownames(first), ]
cor.test(first[["deseq_logfc"]], second[["deseq_logfc"]])

## 
##  Pearson's product-moment correlation
## 
## data:  first[["deseq_logfc"]] and second[["deseq_logfc"]]
## t = 161, df = 10896, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8338 0.8449
## sample estimates:
##    cor 
## 0.8394

tumaco_batch_aucc <- calculate_aucc(cf_tumaco_monocyte_batchvisit_tables[["data"]][["fail_vs_cure"]],
                            tbl2=cf_monocyte_batchvisit_tables[["data"]][["fail_vs_cure"]],
                            py="deseq_adjp", ly="deseq_logfc")
tumaco_batch_aucc

## $aucc
## [1] 0.4004
## 
## $plot

shared_ids <- rownames(cf_tumaco_monocyte_batchvisit_tables[["data"]][["fail_vs_cure"]]) %in%
  rownames(cf_monocyte_batchvisit_tables[["data"]][["fail_vs_cure"]])
first <- cf_monocyte_batchvisit_tables[["data"]][["fail_vs_cure"]][shared_ids, ]
second <- cf_tumaco_monocyte_batchvisit_tables[["data"]][["fail_vs_cure"]][rownames(first), ]
cor.test(first[["deseq_logfc"]], second[["deseq_logfc"]])

## 
##  Pearson's product-moment correlation
## 
## data:  first[["deseq_logfc"]] and second[["deseq_logfc"]]
## t = 144, df = 10896, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8023 0.8153
## sample estimates:
##    cor 
## 0.8089

first_sig_names <- rownames(cf_monocyte_batchvisit_sig[["deseq"]][["ups"]][["fail_vs_cure"]])
second_sig_names <- rownames(cf_tumaco_monocyte_batchvisit_sig[["deseq"]][["ups"]][["fail_vs_cure"]])
batch_venn_lst <- list(all_batch = first_sig_names, tumaco_batch = second_sig_names)
batch_venn <- Vennerable::Venn(batch_venn_lst)
Vennerable::plot(batch_venn)

8.1.1.4 Cure/Fail, Neutrophils

cf_neutrophil_de <- all_pairwise(neutrophil_samples, model_batch = "svaseq", filter = TRUE)

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

## Setting 1703 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

cf_neutrophil_tables <- combine_de_tables(
    cf_neutrophil_de,
    keepers = cf_contrasts,
    excel = glue::glue("excel/cf_neutrophil_tables-v{ver}.xlsx"))

## Deleting the file excel/cf_neutrophil_tables-v202205.xlsx before writing the tables.

cf_neutrophil_sig <- extract_significant_genes(
    cf_neutrophil_tables,
    excel = glue::glue("excel/cf_neutrophil_sig-v{ver}.xlsx"))

## Deleting the file excel/cf_neutrophil_sig-v202205.xlsx before writing the tables.

8.1.1.5 Cure/Fail, Eosinophils

cf_eosinophil_de <- all_pairwise(eosinophil_samples, model_batch = "svaseq", filter = TRUE)

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

## Setting 936 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

cf_eosinophil_tables <- combine_de_tables(
    cf_eosinophil_de,
    keepers = cf_contrasts,
    excel = glue::glue("excel/cf_eosinophil_tables-v{ver}.xlsx"))

## Deleting the file excel/cf_eosinophil_tables-v202205.xlsx before writing the tables.

cf_eosinophil_sig <- extract_significant_genes(
    cf_eosinophil_tables,
    excel = glue::glue("excel/cf_eosinophil_sig-v{ver}.xlsx"))

## Deleting the file excel/cf_eosinophil_sig-v202205.xlsx before writing the tables.

8.1.1.6 Cure/Fail clinical (Not biopsies)

cf_nobiopsy_de <- all_pairwise(hs_clinical_nobiop, model_batch = "svaseq", filter = TRUE)

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

## Setting 18546 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

cf_nobiopsy_tables <- combine_de_tables(
    cf_nobiopsy_de,
    keepers = cf_contrasts,
    excel = glue::glue("excel/cf_nobiopsy_tables-v{ver}.xlsx"))

## Deleting the file excel/cf_nobiopsy_tables-v202205.xlsx before writing the tables.

cf_nobiopsy_sig <- extract_significant_genes(
    cf_nobiopsy_tables,
    excel = glue::glue("excel/cf_nobiopsy_sig-v{ver}.xlsx"))

## Deleting the file excel/cf_nobiopsy_sig-v202205.xlsx before writing the tables.

8.1.2 By visit

For these contrasts, we want to see fail_v1 vs. cure_v1, fail_v2 vs. cure_v2 etc. As a result, we will need to juggle the data slightly and add another set of contrasts.

8.1.3 Setting up

visit_cf_expt_factor <- paste0("v", pData(hs_clinical)[["visitnumber"]],
                               pData(hs_clinical)[["condition"]])
visit_cf_expt <- set_expt_conditions(hs_clinical, fact = visit_cf_expt_factor)

visit_cf_eosinophil <- subset_expt(visit_cf_expt, subset="typeofcells=='eosinophils'")

## subset_expt(): There were 192, now there are 42 samples.

visit_cf_tumaco_eosinophil <- subset_expt(visit_cf_eosinophil, subset="clinic=='Tumaco'")

## subset_expt(): There were 42, now there are 27 samples.

visit_cf_monocyte <- subset_expt(visit_cf_expt, subset="typeofcells=='monocytes'")

## subset_expt(): There were 192, now there are 66 samples.

visit_cf_tumaco_monocyte <- subset_expt(visit_cf_monocyte, subset="clinic=='Tumaco'")

## subset_expt(): There were 66, now there are 45 samples.

visit_cf_neutrophil <- subset_expt(visit_cf_expt, subset="typeofcells=='neutrophils'")

## subset_expt(): There were 192, now there are 65 samples.

visit_cf_tumaco_neutrophil <- subset_expt(visit_cf_neutrophil, subset="clinic=='Tumaco'")

## subset_expt(): There were 65, now there are 44 samples.

8.1.3.1 Cure/Fail by visits, all cell types

visit_cf_all_de <- all_pairwise(visit_cf_expt, model_batch = "svaseq", filter = TRUE)

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

## Setting 30450 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

visit_cf_all_tables <- combine_de_tables(
    visit_cf_all_de,
    keepers = visit_cf_contrasts,
    excel = glue::glue("excel/visit_cf_all_tables-v{ver}.xlsx"))

## Deleting the file excel/visit_cf_all_tables-v202205.xlsx before writing the tables.

visit_cf_all_sig <- extract_significant_genes(
    visit_cf_all_tables,
    excel = glue::glue("excel/visit_cf_all_sig-v{ver}.xlsx"))

## Deleting the file excel/visit_cf_all_sig-v202205.xlsx before writing the tables.

8.1.3.2 Cure/Fail by visit, Monocytes

visit_cf_monocyte_de <- all_pairwise(visit_cf_monocyte, model_batch = "svaseq", filter = TRUE)

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

## Setting 1541 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

visit_cf_monocyte_tables <- combine_de_tables(
    visit_cf_monocyte_de,
    keepers = visit_cf_contrasts,
    excel = glue::glue("excel/visit_cf_monocyte_tables-v{ver}.xlsx"))

## Deleting the file excel/visit_cf_monocyte_tables-v202205.xlsx before writing the tables.

visit_cf_monocyte_sig <- extract_significant_genes(
    visit_cf_monocyte_tables,
    excel = glue::glue("excel/visit_cf_monocyte_sig-v{ver}.xlsx"))

## Deleting the file excel/visit_cf_monocyte_sig-v202205.xlsx before writing the tables.

v1fc_deseq_ma <- visit_cf_monocyte_tables[["plots"]][["v1fail_vs_cure"]][["deseq_ma_plots"]][["plot"]]
dev <- pp(file="images/monocyte_cf_de_v1_maplot.png")
v1fc_deseq_ma
closed <- dev.off()
v1fc_deseq_ma

v2fc_deseq_ma <- visit_cf_monocyte_tables[["plots"]][["v2fail_vs_cure"]][["deseq_ma_plots"]][["plot"]]
dev <- pp(file="images/monocyte_cf_de_v2_maplot.png")
v2fc_deseq_ma
closed <- dev.off()
v2fc_deseq_ma

v3fc_deseq_ma <- visit_cf_monocyte_tables[["plots"]][["v3fail_vs_cure"]][["deseq_ma_plots"]][["plot"]]
dev <- pp(file="images/monocyte_cf_de_v3_maplot.png")
v3fc_deseq_ma
closed <- dev.off()
v3fc_deseq_ma

## Repeat for the tumaco subset
visit_cf_tumaco_monocyte_de <- all_pairwise(visit_cf_tumaco_monocyte,
                                          model_batch = "svaseq", filter = TRUE)

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

## Setting 782 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

visit_cf_tumaco_monocyte_tables <- combine_de_tables(
    visit_cf_tumaco_monocyte_de,
    keepers = visit_cf_contrasts,
    excel = glue::glue("excel/visit_cf_tumaco_monocyte_tables-v{ver}.xlsx"))

## Deleting the file excel/visit_cf_tumaco_monocyte_tables-v202205.xlsx before writing the tables.

visit_cf_tumaco_monocyte_sig <- extract_significant_genes(
    visit_cf_tumaco_monocyte_tables,
    excel = glue::glue("excel/visit_cf_tumaco_monocyte_sig-v{ver}.xlsx"))

## Deleting the file excel/visit_cf_tumaco_monocyte_sig-v202205.xlsx before writing the tables.

One query from Alejandro is to look at the genes shared up/down across visits. I am not entirely certain we have enough samples for this to work, but let us find out.

I am thinking this is a good place to use the AUCC curves I learned about thanks to Julie Cridland.

Note that the following is all monocyte samples, this should therefore potentially be moved up and a version of this with only the Tumaco samples put here?

v1fc <- visit_cf_monocyte_tables[["data"]][["v1fail_vs_cure"]]
v2fc <- visit_cf_monocyte_tables[["data"]][["v2fail_vs_cure"]]
v3fc <- visit_cf_monocyte_tables[["data"]][["v3fail_vs_cure"]]

v1_sig <- c(
    rownames(visit_cf_monocyte_sig[["deseq"]][["ups"]][["v1fail_vs_cure"]]),
    rownames(visit_cf_monocyte_sig[["deseq"]][["downs"]][["v1fail_vs_cure"]]))
v2_sig <- c(
    rownames(visit_cf_monocyte_sig[["deseq"]][["ups"]][["v2fail_vs_cure"]]),
    rownames(visit_cf_monocyte_sig[["deseq"]][["downs"]][["v2fail_vs_cure"]]))
v3_sig <- c(
    rownames(visit_cf_monocyte_sig[["deseq"]][["ups"]][["v2fail_vs_cure"]]),
    rownames(visit_cf_monocyte_sig[["deseq"]][["downs"]][["v2fail_vs_cure"]]))

monocyte_visit_aucc_v2v1 <- calculate_aucc(v1fc, tbl2=v2fc, py="deseq_adjp", ly="deseq_logfc")
dev <- pp(file="images/monocyte_visit_v2v1_aucc.png")
monocyte_visit_aucc_v2v1[["plot"]]
closed <- dev.off()
monocyte_visit_aucc_v2v1[["plot"]]

monocyte_visit_aucc_v3v1 <- calculate_aucc(v1fc, tbl2=v3fc, py="deseq_adjp", ly="deseq_logfc")
dev <- pp(file="images/monocyte_visit_v3v1_aucc.png")
monocyte_visit_aucc_v3v1[["plot"]]
closed <- dev.off()
monocyte_visit_aucc_v3v1[["plot"]]

8.1.3.3 Repeat for only the Tumaco samples

v1fc_tumaco <- visit_cf_tumaco_monocyte_tables[["data"]][["v1fail_vs_cure"]]
v2fc_tumaco <- visit_cf_tumaco_monocyte_tables[["data"]][["v2fail_vs_cure"]]
v3fc_tumaco <- visit_cf_tumaco_monocyte_tables[["data"]][["v3fail_vs_cure"]]

v1_tumaco_sig <- c(
    rownames(visit_cf_tumaco_monocyte_sig[["deseq"]][["ups"]][["v1fail_vs_cure"]]),
    rownames(visit_cf_tumaco_monocyte_sig[["deseq"]][["downs"]][["v1fail_vs_cure"]]))
v2_tumaco_sig <- c(
    rownames(visit_cf_tumaco_monocyte_sig[["deseq"]][["ups"]][["v2fail_vs_cure"]]),
    rownames(visit_cf_tumaco_monocyte_sig[["deseq"]][["downs"]][["v2fail_vs_cure"]]))
v3_tumaco_sig <- c(
    rownames(visit_cf_tumaco_monocyte_sig[["deseq"]][["ups"]][["v2fail_vs_cure"]]),
    rownames(visit_cf_tumaco_monocyte_sig[["deseq"]][["downs"]][["v2fail_vs_cure"]]))

monocyte_tumaco_visit_aucc_v2v1 <- calculate_aucc(v1fc_tumaco, tbl2=v2fc_tumaco,
                                                  py="deseq_adjp", ly="deseq_logfc")
dev <- pp(file="images/monocyte_visit_tumaco_v2v1_aucc.png")
monocyte_tumaco_visit_aucc_v2v1[["plot"]]
closed <- dev.off()
monocyte_tumaco_visit_aucc_v2v1[["plot"]]

monocyte_tumaco_visit_aucc_v3v1 <- calculate_aucc(v1fc_tumaco, tbl2=v3fc_tumaco,
                                                  py="deseq_adjp", ly="deseq_logfc")
dev <- pp(file="images/monocyte_visit_tumaco_v3v1_aucc.png")
monocyte_tumaco_visit_aucc_v3v1[["plot"]]
closed <- dev.off()
monocyte_tumaco_visit_aucc_v3v1[["plot"]]

monocyte_tumaco_visit_aucc_v3v2 <- calculate_aucc(v3fc_tumaco, tbl2=v2fc_tumaco,
                                                  py="deseq_adjp", ly="deseq_logfc")
dev <- pp(file="images/monocyte_visit_tumaco_v3v2_aucc.png")
monocyte_tumaco_visit_aucc_v3v2[["plot"]]
closed <- dev.off()
monocyte_tumaco_visit_aucc_v3v2[["plot"]]

8.1.3.4 Cure/Fail by visit, Neutrophils

a

visit_cf_neutrophil_de <- all_pairwise(visit_cf_neutrophil, model_batch = "svaseq", filter = TRUE)

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

## Setting 1612 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

visit_cf_neutrophil_tables <- combine_de_tables(
    visit_cf_neutrophil_de,
    keepers = visit_cf_contrasts,
    excel = glue::glue("excel/visit_cf_neutrophil_tables-v{ver}.xlsx"))

## Deleting the file excel/visit_cf_neutrophil_tables-v202205.xlsx before writing the tables.

visit_cf_neutrophil_sig <- extract_significant_genes(
    visit_cf_neutrophil_tables,
    excel = glue::glue("excel/visit_cf_neutrophil_sig-v{ver}.xlsx"))

## Deleting the file excel/visit_cf_neutrophil_sig-v202205.xlsx before writing the tables.

8.1.3.5 Cure/Fail by visit, Eosinophils

visit_cf_eosinophil_de <- all_pairwise(visit_cf_eosinophil, model_batch = "svaseq", filter = TRUE)

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

## Setting 980 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

visit_cf_eosinophil_tables <- combine_de_tables(
    visit_cf_eosinophil_de,
    keepers = visit_cf_contrasts,
    excel = glue::glue("excel/visit_cf_eosinophil_tables-v{ver}.xlsx"))

## Deleting the file excel/visit_cf_eosinophil_tables-v202205.xlsx before writing the tables.

visit_cf_eosinophil_sig <- extract_significant_genes(
    visit_cf_eosinophil_tables,
    excel = glue::glue("excel/visit_cf_eosinophil_sig-v{ver}.xlsx"))

## Deleting the file excel/visit_cf_eosinophil_sig-v202205.xlsx before writing the tables.

8.2 Persistence in visit 3

Having put some SL read mapping information in the sample sheet, Maria Adelaida added a new column using it with the putative persistence state on a per-sample basis. One question which arised from that: what differences are observable between the persistent yes vs. no samples on a per-cell-type basis among the visit 3 samples.

8.2.1 Setting up

First things first, create the datasets.

persistence_expt <- subset_expt(hs_clinical, subset = "persistence=='Y'|persistence=='N'") %>%
  subset_expt(subset = 'visitnumber==3') %>%
  set_expt_conditions(fact = 'persistence')

## subset_expt(): There were 192, now there are 118 samples.

## subset_expt(): There were 118, now there are 40 samples.

## persistence_biopsy <- subset_expt(persistence_expt, subset = "typeofcells=='biopsy'")
persistence_monocyte <- subset_expt(persistence_expt, subset = "typeofcells=='monocytes'")

## subset_expt(): There were 40, now there are 16 samples.

persistence_neutrophil <- subset_expt(persistence_expt, subset = "typeofcells=='neutrophils'")

## subset_expt(): There were 40, now there are 14 samples.

persistence_eosinophil <- subset_expt(persistence_expt, subset = "typeofcells=='eosinophils'")

## subset_expt(): There were 40, now there are 10 samples.

8.2.2 Take a look

See if there are any patterns which look usable.

## All
persistence_norm <- normalize_expt(persistence_expt, transform = "log2", convert = "cpm",
                                   norm = "quant", filter = TRUE)

## Removing 8378 low-count genes (11545 remaining).

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

plot_pca(persistence_norm)$plot

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

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

## Removing 8378 low-count genes (11545 remaining).

## Setting 2392 low elements to zero.

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

plot_pca(persistence_nb)$plot

## Biopsies
##persistence_biopsy_norm <- normalize_expt(persistence_biopsy, transform = "log2", convert = "cpm",
##                                   norm = "quant", filter = TRUE)
##plot_pca(persistence_biopsy_norm)$plot
## Insufficient data

## Monocytes
persistence_monocyte_norm <- normalize_expt(persistence_monocyte, transform = "log2", convert = "cpm",
                                   norm = "quant", filter = TRUE)

## Removing 9391 low-count genes (10532 remaining).

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

plot_pca(persistence_monocyte_norm)$plot

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

## Removing 9391 low-count genes (10532 remaining).

## Setting 106 low elements to zero.

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

plot_pca(persistence_monocyte_nb)$plot

## Neutrophils
persistence_neutrophil_norm <- normalize_expt(persistence_neutrophil, transform = "log2", convert = "cpm",
                                   norm = "quant", filter = TRUE)

## Removing 11372 low-count genes (8551 remaining).

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

plot_pca(persistence_neutrophil_norm)$plot

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

## Removing 11372 low-count genes (8551 remaining).

## Setting 138 low elements to zero.

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

plot_pca(persistence_neutrophil_nb)$plot

## Eosinophils
persistence_eosinophil_norm <- normalize_expt(persistence_eosinophil, transform = "log2", convert = "cpm",
                                   norm = "quant", filter = TRUE)

## Removing 9757 low-count genes (10166 remaining).

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

plot_pca(persistence_eosinophil_norm)$plot

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

## Removing 9757 low-count genes (10166 remaining).

## Setting 60 low elements to zero.

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

plot_pca(persistence_eosinophil_nb)$plot

8.2.3 persistence DE

persistence_de <- all_pairwise(persistence_expt, filter = TRUE, model_batch = "svaseq")

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

## Setting 2392 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

persistence_table <- combine_de_tables(
    persistence_de,
    excel = glue::glue("excel/persistence_all_de-v{ver}.xlsx"))

## Deleting the file excel/persistence_all_de-v202205.xlsx before writing the tables.

persistence_monocyte_de <- all_pairwise(persistence_monocyte, filter = TRUE, model_batch = "svaseq")

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

## Setting 106 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

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

## Deleting the file excel/persistence_monocyte_de-v202205.xlsx before writing the tables.

persistence_neutrophil_de <- all_pairwise(persistence_neutrophil, filter = TRUE, model_batch = "svaseq")

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

## Setting 138 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

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

## Deleting the file excel/persistence_neutrophil_de-v202205.xlsx before writing the tables.

persistence_eosinophil_de <- all_pairwise(persistence_eosinophil, filter = TRUE, model_batch = "svaseq")

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

## Setting 60 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

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

## Deleting the file excel/persistence_eosinophil_de-v202205.xlsx before writing the tables.

8.3 Comparing visits without regard to cure/fail

8.3.1 All cell types

visit_all_de <- all_pairwise(visit_expt, filter = TRUE, model_batch = "svaseq")

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

## Setting 18454 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

visit_all_table <- combine_de_tables(
    visit_all_de,
    keepers = visit_contrasts,
    excel = glue::glue("excel/visit_all_de-v{ver}.xlsx"))

## Deleting the file excel/visit_all_de-v202205.xlsx before writing the tables.

visit_all_sig <- extract_significant_genes(
    visit_all_table,
    excel = glue::glue("excel/visit_all_sig-v{ver}.xlsx"))

## Deleting the file excel/visit_all_sig-v202205.xlsx before writing the tables.

visit_all_saved <- save(list = "visit_all_table",
                        file = glue::glue("rda/visit_all_table-v{ver}.rda"))

8.3.2 Monocyte samples

visit_monocyte_de <- all_pairwise(visit_monocyte, filter = TRUE, model_batch = "svaseq")

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

## Setting 1324 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

visit_monocyte_table <- combine_de_tables(
    visit_monocyte_de,
    keepers = visit_contrasts,
    excel = glue::glue("excel/visit_monocyte_de-v{ver}.xlsx"))

## Deleting the file excel/visit_monocyte_de-v202205.xlsx before writing the tables.

visit_monocyte_sig <- extract_significant_genes(
    visit_monocyte_table,
    excel = glue::glue("excel/visit_monocyte_sig-v{ver}.xlsx"))

## Deleting the file excel/visit_monocyte_sig-v202205.xlsx before writing the tables.

8.3.3 Neutrophil samples

visit_neutrophil_de <- all_pairwise(visit_neutrophil, filter = TRUE, model_batch = "svaseq")

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

## Setting 1544 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

visit_neutrophil_table <- combine_de_tables(
    visit_neutrophil_de,
    keepers = visit_contrasts,
    excel = glue::glue("excel/visit_neutrophil_de-v{ver}.xlsx"))

## Deleting the file excel/visit_neutrophil_de-v202205.xlsx before writing the tables.

visit_neutrophil_sig <- extract_significant_genes(
    visit_neutrophil_table,
    excel = glue::glue("excel/visit_neutrophil_sig-v{ver}.xlsx"))

## Deleting the file excel/visit_neutrophil_sig-v202205.xlsx before writing the tables.

visit_neutrophil_saved <- save(list = "visit_neutrophil_table",
                               file = glue::glue("rda/visit_neutrophil_table-v{ver}.rda"))

8.3.4 Eosinophil samples

visit_eosinophil_de <- all_pairwise(visit_eosinophil, filter = TRUE, model_batch = "svaseq")

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

## Setting 983 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

visit_eosinophil_table <- combine_de_tables(
    visit_eosinophil_de,
    keepers = visit_contrasts,
    excel = glue::glue("excel/visit_eosinophil_de-v{ver}.xlsx"))

## Deleting the file excel/visit_eosinophil_de-v202205.xlsx before writing the tables.

visit_eosinophil_sig <- extract_significant_genes(
    visit_eosinophil_table,
    excel = glue::glue("excel/visit_eosinophil_sig-v{ver}.xlsx"))

## Deleting the file excel/visit_eosinophil_sig-v202205.xlsx before writing the tables.

visit_eosinophil_saved <- save(list = "visit_eosinophil_table",
                               file = glue::glue("rda/visit_eosinophil_table-v{ver}.rda"))

9 Compare the two clinics directly

monocytes_by_place <- set_expt_conditions(monocyte_samples, fact="clinic")
eosinophils_by_place <- set_expt_conditions(eosinophil_samples, fact="clinic")
neutrophils_by_place <- set_expt_conditions(neutrophil_samples, fact="clinic")
biopsies_by_place <- set_expt_conditions(biopsy_samples, fact="clinic")

monocyte_place_de <- all_pairwise(monocytes_by_place, model_batch=TRUE)

## Using limma's removeBatchEffect to visualize with(out) batch inclusion.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

monocyte_table <- combine_de_tables(
    monocyte_place_de,
    excel="excel/monocytes_by_place-table.xlsx")

## Deleting the file excel/monocytes_by_place-table.xlsx before writing the tables.

monocyte_sig <- extract_significant_genes(
    monocyte_table,
    excel="excel/monocytes_by_place-sig.xlsx")

## Deleting the file excel/monocytes_by_place-sig.xlsx before writing the tables.

eosinophil_place_de <- all_pairwise(eosinophils_by_place, model_batch=TRUE)

## Using limma's removeBatchEffect to visualize with(out) batch inclusion.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

eosinophil_table <- combine_de_tables(
    eosinophil_place_de,
    excel="excel/eosinophils_by_place-table.xlsx")

## Deleting the file excel/eosinophils_by_place-table.xlsx before writing the tables.

eosinophil_sig <- extract_significant_genes(
    eosinophil_table,
    excel="excel/eosinophils_by_place-sig.xlsx")

## Deleting the file excel/eosinophils_by_place-sig.xlsx before writing the tables.

neutrophil_place_de <- all_pairwise(neutrophils_by_place, model_batch=TRUE)

## Using limma's removeBatchEffect to visualize with(out) batch inclusion.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

neutrophil_table <- combine_de_tables(
    neutrophil_place_de,
    excel="excel/neutrophils_by_place-table.xlsx")

## Deleting the file excel/neutrophils_by_place-table.xlsx before writing the tables.

neutrophil_sig <- extract_significant_genes(
    neutrophil_table,
    excel="excel/neutrophils_by_place-sig.xlsx")

## Deleting the file excel/neutrophils_by_place-sig.xlsx before writing the tables.

biopsy_place_de <- all_pairwise(biopsies_by_place, model_batch=TRUE)

## Using limma's removeBatchEffect to visualize with(out) batch inclusion.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

biopsy_table <- combine_de_tables(
    biopsy_place_de,
    excel="biopsies_by_place-table.xlsx")

## Deleting the file biopsies_by_place-table.xlsx before writing the tables.

biopsy_sig <- extract_significant_genes(
    biopsy_table,
    excel="biopsies_by_place-sig.xlsx")

## Deleting the file biopsies_by_place-sig.xlsx before writing the tables.

10 Likelihood Ratio Testing (LRT)

10.1 Shared patterns across visits

clinical_filt <- normalize_expt(hs_clinical, filter = TRUE)

## Removing 5601 low-count genes (14322 remaining).

lrt_visit <- deseq_lrt(clinical_filt, transform = "vst", interaction = FALSE,
                       interactor_column = "visitnumber",
                       interest_column = "clinicaloutcome")

## converting counts to integer mode

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 195 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

## A large number of genes was given-- please, make sure this is not an error. Normally, only DE genes will be useful for this function.

## Working with 5979 genes.

## Working with 5979 genes after filtering: minc > 3

## Joining, by = "merge"
## Joining, by = "merge"

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.

summary(lrt_visit[["favorite_genes"]])

## Length  Class   Mode 
##      0   NULL   NULL

written <- write_xlsx(data=as.data.frame(lrt_visit[["deseq_table"]]),
                      excel=glue::glue("excel/lrt_clinical_visit-v{ver}.xlsx"))

lrt_monocyte_visit <- deseq_lrt(visit_monocyte, transform = "vst",
                                interaction = FALSE,
                                interactor_column = "visitnumber",
                                interest_column = "clinicaloutcome")

## converting counts to integer mode
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 132 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
## Working with 42 genes.
## Working with 41 genes after filtering: minc > 3
## Joining, by = "merge"Joining, by = "merge"

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.

lrt_monocyte_visit$cluster_data$plot

## NULL

lrt_tumaco_monocyte_visit <- deseq_lrt(tumaco_monocyte, transform = "vst",
                                       interaction = FALSE, minc = 1,
                                       interactor_column = "visitnumber",
                                       interest_column = "clinicaloutcome")

## converting counts to integer mode
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 58 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
## Working with 17 genes.
## Working with 17 genes after filtering: minc > 1
## Joining, by = "merge"Joining, by = "merge"

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.

10.2 Shared patterns across cell types

lrt_celltype_clinical_test <- deseq_lrt(hs_clinical, transform = "vst",
                                        interactor_column = "typeofcells",
                                        interest_column = "clinicaloutcome")

## converting counts to integer mode

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 21 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

## Working with 450 genes.

## Working with 446 genes after filtering: minc > 3

## Joining, by = "merge"
## Joining, by = "merge"

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.

deseq_lrt_df <- merge(hs_annot, as.data.frame(lrt_celltype_clinical_test[["deseq_table"]]), all.y=TRUE,
                      by="row.names")
rownames(deseq_lrt_df) <- deseq_lrt_df[["Row.names"]]
deseq_lrt_df[["Row.names"]] <- NULL
written <- write_xlsx(data=deseq_lrt_df,
                      excel=glue::glue("excel/lrt_clinical_celltype-v{ver}.xlsx"))

11 Gene Set Analyses

The gene sets we are most likely to consider are the results of the various preceding differential expression analyses.

11.1 GSEA

11.1.1 Cure/Fail groups

In the context of cure vs. fail, we have mixed and matched the data in quite a few different ways.

11.1.1.1 Cure/Fail, all samples

For the moment, let us assume that the default definition of ‘significant’ is sufficient for these analyses, with the knowledge that is not likely to remain true.

The resulting data structures are organized as:

thing[[“deseq”]][[ups"]][[contrast]]

In addition, we have a few ontology-esque tools available, so let us mix and match and see what pops out.

cf_all_up <- cf_clinical_sig[["deseq"]][["ups"]][["fail_vs_cure"]]
dim(cf_all_up)

## [1] 122  58

cf_all_down <- cf_clinical_sig[["deseq"]][["downs"]][["fail_vs_cure"]]
dim(cf_all_down)

## [1] 50 58

cf_all_up_gp <- simple_gprofiler(cf_all_up)

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

## GO search found 80 hits.

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

## KEGG search found 24 hits.

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

## REAC search found 7 hits.

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

## MI search found 0 hits.

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

## TF search found 8 hits.

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

## CORUM search found 6 hits.

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

## HP search found 0 hits.

cf_all_up_gp[["pvalue_plots"]][["reactome_plot_over"]]

cf_all_up_gp[["pvalue_plots"]][["kegg_plot_over"]]

cf_all_up_gp[["pvalue_plots"]][["mfp_plot_over"]]

cf_all_up_gp[["pvalue_plots"]][["bpp_plot_over"]]

cf_all_down_gp <- simple_gprofiler(cf_all_down)

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

## GO search found 3 hits.

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

## KEGG search found 2 hits.

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

## REAC search found 5 hits.

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

## MI search found 0 hits.

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

## TF search found 0 hits.

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

## CORUM search found 4 hits.

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

## HP search found 0 hits.

cf_all_down_gp[["pvalue_plots"]][["reactome_plot_over"]]

## NULL

cf_all_down_gp[["pvalue_plots"]][["kegg_plot_over"]]

cf_all_down_gp[["pvalue_plots"]][["mfp_plot_over"]]

## NULL

cf_all_down_gp[["pvalue_plots"]][["bpp_plot_over"]]

11.1.1.2 Cure/Fail, Biopsies

There are no genes deemed significant in the biopsy comparisons of cure/fail. Oddly, when I looked manually, I thought I saw one candidate gene. In addition, when I stepped through the function manually I got the same gene…

11.1.1.3 Cure/Fail, Monocytes

Try again with monocytes.

cf_monocyte_up <- cf_monocyte_sva_sig[["deseq"]][["ups"]][["fail_vs_cure"]]
dim(cf_monocyte_up)

## [1] 60 58

cf_monocyte_down <- cf_monocyte_sva_sig[["deseq"]][["downs"]][["fail_vs_cure"]]
dim(cf_monocyte_down)

## [1] 33 58

cf_monocyte_up_gp <- simple_gprofiler(cf_monocyte_up)

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

## GO search found 56 hits.

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

## KEGG search found 14 hits.

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

## REAC search found 9 hits.

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

## MI search found 0 hits.

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

## TF search found 2 hits.

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

## CORUM search found 1 hits.

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

## HP search found 35 hits.

cf_monocyte_up_gp[["pvalue_plots"]][["reactome_plot_over"]]

cf_monocyte_up_gp[["pvalue_plots"]][["kegg_plot_over"]]

cf_monocyte_up_gp[["pvalue_plots"]][["bpp_plot_over"]]

cf_monocyte_down_gp <- simple_gprofiler(cf_monocyte_down)

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

## GO search found 0 hits.

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

## KEGG search found 1 hits.

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

## REAC search found 2 hits.

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

## MI search found 0 hits.

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

## TF search found 0 hits.

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

## CORUM search found 2 hits.

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

## HP search found 0 hits.

cf_monocyte_down_gp[["pvalue_plots"]][["reactome_plot_over"]]

cf_monocyte_down_gp[["pvalue_plots"]][["kegg_plot_over"]]

11.1.1.4 Cure/Fail, Neutrophils

cf_neutrophil_up <- cf_neutrophil_sig[["deseq"]][["ups"]][["fail_vs_cure"]]
dim(cf_neutrophil_up)

## [1] 27 58

cf_neutrophil_down <- cf_neutrophil_sig[["deseq"]][["downs"]][["fail_vs_cure"]]
dim(cf_neutrophil_down)

## [1] 19 58

cf_neutrophil_up_gp <- simple_gprofiler(cf_neutrophil_up)

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

## GO search found 0 hits.

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

## KEGG search found 0 hits.

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

## REAC search found 0 hits.

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

## MI search found 0 hits.

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

## TF search found 11 hits.

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

## CORUM search found 1 hits.

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

## HP search found 0 hits.

cf_neutrophil_down_gp <- simple_gprofiler(cf_neutrophil_down)

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

## GO search found 0 hits.

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

## KEGG search found 1 hits.

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

## REAC search found 0 hits.

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

## MI search found 0 hits.

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

## TF search found 0 hits.

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

## CORUM search found 0 hits.

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

## HP search found 0 hits.

cf_neutrophil_down_gp[["pvalue_plots"]][["kegg_plot_over"]]

11.1.1.5 Cure/Fail, Eosinophils

cf_eosinophil_up <- cf_eosinophil_sig[["deseq"]][["ups"]][["fail_vs_cure"]]
dim(cf_eosinophil_up)

## [1] 207  58

cf_eosinophil_down <- cf_eosinophil_sig[["deseq"]][["downs"]][["fail_vs_cure"]]
dim(cf_eosinophil_down)

## [1] 173  58

cf_eosinophil_up_gp <- simple_gprofiler(cf_eosinophil_up)

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

## GO search found 125 hits.

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

## KEGG search found 8 hits.

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

## REAC search found 10 hits.

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

## MI search found 0 hits.

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

## TF search found 22 hits.

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

## CORUM search found 5 hits.

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

## HP search found 0 hits.

cf_eosinophil_up_gp[["pvalue_plots"]][["kegg_plot_over"]]

cf_eosinophil_up_gp[["pvalue_plots"]][["reactome_plot_over"]]

cf_eosinophil_up_gp[["pvalue_plots"]][["mfp_plot_over"]]

cf_eosinophil_up_gp[["pvalue_plots"]][["bpp_plot_over"]]

cf_eosinophil_down_gp <- simple_gprofiler(cf_eosinophil_down)

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

## GO search found 44 hits.

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

## KEGG search found 1 hits.

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

## REAC search found 19 hits.

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

## MI search found 0 hits.

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

## TF search found 0 hits.

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

## CORUM search found 2 hits.

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

## HP search found 1 hits.

cf_eosinophil_down_gp[["pvalue_plots"]][["reactome_plot_over"]]

cf_eosinophil_down_gp[["pvalue_plots"]][["mfp_plot_over"]]

cf_eosinophil_down_gp[["pvalue_plots"]][["bpp_plot_over"]]

11.1.1.6 Clinical samples

cf_nobiopsy_up <- cf_nobiopsy_sig[["deseq"]][["ups"]][["fail_vs_cure"]]
dim(cf_nobiopsy_up)

## [1] 181  58

cf_nobiopsy_down <- cf_nobiopsy_sig[["deseq"]][["downs"]][["fail_vs_cure"]]
dim(cf_nobiopsy_down)

## [1] 98 58

cf_nobiopsy_up_gp <- simple_gprofiler(cf_nobiopsy_up)

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

## GO search found 131 hits.

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

## KEGG search found 11 hits.

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

## REAC search found 5 hits.

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

## MI search found 0 hits.

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

## TF search found 12 hits.

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

## CORUM search found 1 hits.

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

## HP search found 0 hits.

cf_nobiopsy_up_gp[["pvalue_plots"]][["kegg_plot_over"]]

cf_nobiopsy_up_gp[["pvalue_plots"]][["reactome_plot_over"]]

cf_nobiopsy_up_gp[["pvalue_plots"]][["mfp_plot_over"]]

cf_nobiopsy_up_gp[["pvalue_plots"]][["bpp_plot_over"]]

cf_nobiopsy_down_gp <- simple_gprofiler(cf_nobiopsy_down)

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

## GO search found 5 hits.

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

## KEGG search found 6 hits.

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

## REAC search found 2 hits.

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

## MI search found 0 hits.

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

## TF search found 0 hits.

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

## CORUM search found 0 hits.

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

## HP search found 0 hits.

cf_nobiopsy_down_gp[["pvalue_plots"]][["reactome_plot_over"]]

cf_nobiopsy_down_gp[["pvalue_plots"]][["mfp_plot_over"]]

## NULL

cf_nobiopsy_down_gp[["pvalue_plots"]][["bpp_plot_over"]]

11.1.2 Cell type groups

11.1.3 Visit groups

11.2 GSVA

11.2.1 Clinical samples

hs_celltype_gsva_c2 <- simple_gsva(hs_valid)

## Converting the rownames() of the expressionset to ENTREZID.

## 588 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the expressionset has 19923 entries.

## After conversion, the expressionset has 19495 entries.

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

## Starting limma pairwise comparison.

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Choosing the non-intercept containing model.

## Assuming this data is similar to a micro array and not performign voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Limma step 5/6: Running eBayes with robust = FALSE and trend = FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: failure_vs_cure.  Adjust = BH

## Limma step 6/6: 1/2: Creating table: cure.  Adjust = BH

## Limma step 6/6: 2/2: Creating table: failure.  Adjust = BH

## The factor cure has 132 rows.

## The factor failure has 60 rows.

## plot labels was not set and there are more than 100 samples, disabling it.

## Testing each factor against the others.

## Scoring cure against everything else.

## Scoring failure against everything else.

## Deleting the file excel/individual_celltypes_gsva_c2.xlsx before writing the tables.

hs_celltype_gsva_c2_sig$subset_plot

hs_celltype_gsva_c2_sig$score_plot

broad_c7 <- load_gmt_signatures(signatures="reference/msigdb/c7.all.v7.2.entrez.gmt",
                                signature_category="c7")
hs_celltype_gsva_c7 <- simple_gsva(hs_valid,
                                   signatures="reference/msigdb/c7.all.v7.2.entrez.gmt",
                                   signature_category="c7",
                                   msig_xml="reference/msigdb_v7.2.xml",
                                   cores=10)

## Converting the rownames() of the expressionset to ENTREZID.

## 588 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the expressionset has 19923 entries.

## After conversion, the expressionset has 19495 entries.

## Adding annotations from reference/msigdb_v7.2.xml.

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

## Starting limma pairwise comparison.

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Choosing the non-intercept containing model.

## Assuming this data is similar to a micro array and not performign voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Limma step 5/6: Running eBayes with robust = FALSE and trend = FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: failure_vs_cure.  Adjust = BH

## Limma step 6/6: 1/2: Creating table: cure.  Adjust = BH

## Limma step 6/6: 2/2: Creating table: failure.  Adjust = BH

## The factor cure has 132 rows.

## The factor failure has 60 rows.

## plot labels was not set and there are more than 100 samples, disabling it.

## Testing each factor against the others.

## Scoring cure against everything else.

## Scoring failure against everything else.

## Deleting the file excel/individual_celltypes_gsva_c7.xlsx before writing the tables.

11.2.2 Print some plots of the GSVA outputs

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

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

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

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

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

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

12 Concordance

Let us compare various results to see how well they agreed

monocyte_cor_subset <- cf_monocyte_sva_tables[["data"]][["fail_vs_cure"]][, c("deseq_logfc", "deseq_adjp")]
neutrophil_cor_subset <- cf_neutrophil_tables[["data"]][["fail_vs_cure"]][, c("deseq_logfc", "deseq_adjp")]
eosinophil_cor_subset <- cf_eosinophil_tables[["data"]][["fail_vs_cure"]][, c("deseq_logfc", "deseq_adjp")]
mono_neut_cor <- merge(monocyte_cor_subset, neutrophil_cor_subset, by="row.names")
mono_eo_cor <- merge(monocyte_cor_subset, eosinophil_cor_subset, by="row.names")
neut_eo_cor <- merge(neutrophil_cor_subset, eosinophil_cor_subset, by="row.names")

cor.test(mono_neut_cor[["deseq_logfc.x"]], mono_neut_cor[["deseq_logfc.y"]])

## 
##  Pearson's product-moment correlation
## 
## data:  mono_neut_cor[["deseq_logfc.x"]] and mono_neut_cor[["deseq_logfc.y"]]
## t = 40, df = 8799, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3778 0.4130
## sample estimates:
##    cor 
## 0.3956

monocyte_neutrophil_aucc <- calculate_aucc(cf_monocyte_sva_tables[["data"]][["fail_vs_cure"]],
                                           cf_neutrophil_tables[["data"]][["fail_vs_cure"]],
                                           px="deseq_p", py="deseq_p", ly="deseq_logfc", lx="deseq_logfc")
monocyte_neutrophil_aucc$plot

cor.test(mono_eo_cor[["deseq_logfc.x"]], mono_eo_cor[["deseq_logfc.y"]])

## 
##  Pearson's product-moment correlation
## 
## data:  mono_eo_cor[["deseq_logfc.x"]] and mono_eo_cor[["deseq_logfc.y"]]
## t = 51, df = 10079, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4344 0.4655
## sample estimates:
##    cor 
## 0.4501

monocyte_eosinophil_aucc <- calculate_aucc(cf_monocyte_sva_tables[["data"]][["fail_vs_cure"]],
                                           cf_eosinophil_tables[["data"]][["fail_vs_cure"]],
                                           px="deseq_p", py="deseq_p", ly="deseq_logfc", lx="deseq_logfc")
monocyte_eosinophil_aucc$plot

cor.test(neut_eo_cor[["deseq_logfc.x"]], neut_eo_cor[["deseq_logfc.y"]])

## 
##  Pearson's product-moment correlation
## 
## data:  neut_eo_cor[["deseq_logfc.x"]] and neut_eo_cor[["deseq_logfc.y"]]
## t = 34, df = 8787, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3223 0.3593
## sample estimates:
##   cor 
## 0.341

neutrophil_eosinophil_aucc <- calculate_aucc(cf_neutrophil_tables[["data"]][["fail_vs_cure"]],
                                           cf_eosinophil_tables[["data"]][["fail_vs_cure"]],
                                           px="deseq_p", py="deseq_p", ly="deseq_logfc", lx="deseq_logfc")
neutrophil_eosinophil_aucc$plot

13 Classify me!

I wrote out all the z2.2 and z2.3 specific variants to a couple files, I want to see if I can classify a human sample as infected with 2.2 or 2.3.

z22 <- read.csv("csv/variants_22.csv")
z23 <- read.csv("csv/variants_23.csv")
cure <- read.csv("csv/cure_variants.txt")
fail <- read.csv("csv/fail_variants.txt")
z22_vec <- gsub(pattern="\\-", replacement="_", x=z22[["x"]])
z23_vec <- gsub(pattern="\\-", replacement="_", x=z23[["x"]])
cure_vec <- gsub(pattern="\\-", replacement="_", x=cure)
fail_vec <- gsub(pattern="\\-", replacement="_", x=fail)

classify_zymo <- function(sample) {
  arbitrary_tags <- sm(readr::read_tsv(sample))
  arbitrary_ids <- arbitrary_tags[["position"]]
  message("Length: ", length(arbitrary_ids), ", z22: ",
          sum(arbitrary_ids %in% z22_vec) / (length(z22_vec)), " z23: ",
          sum(arbitrary_ids %in% z23_vec) / (length(z23_vec)))
}

arbitrary_sample <- "preprocessing/TMRC30156/outputs/40freebayes_lpanamensis_v36/all_tags.txt.xz"
classify_zymo(arbitrary_sample)

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))
}

tmp <- loadme(filename=savefile)