L. panamensis 202011: A Macrophage Side Experiment

1 Introduction

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

2 Annotation

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

hs_annot <- load_biomart_annotations()

## 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")]

3 Sample Estimation

I used two mapping methods for this data, hisat2 and salmon. Most analyses use hisat2, which is a more traditional map-and-count method. In contrast, salmon uses what may be thought of as a indexed voting method (so that multi-matches are discounted and the votes split among all matches). Salmon also required a pre-existing database of known transcripts (though later versions may actually use mapping from things like hisat), while hisat uses the genome and a database of known transcripts (and optionally can search for splicing junctions to find new transcripts).

3.1 Generate expressionsets

Caveat: This initial section is using salmon quantifications. A majority of analyses used hisat2.

3.1.1 Salmon expressionsets

Currently, I have these disabled.

hs_expt <- sm(create_expt("sample_sheets/tmrc3_samples_20201105.xlsx",
                          file_column="hg3891salmonfile",
                          gene_info=hs_annot, tx_gene_map=tx_gene_map))

libsizes <- plot_libsize(hs_expt)
libsizes$plot
nonzero <- plot_nonzero(hs_expt)
box <- plot_boxplot(hs_expt)

hs_write <- write_expt(hs_expt, excel=glue("excel/hs_written_salmon-v{ver}.xlsx"))

hs_valid <- subset_expt(hs_expt, coverage=100000)
valid_write <- write_expt(hs_valid, excel=glue("excel/hs_valid_salmon-v{ver}.xlsx"))

3.1.2 Hisat2 expressionsets

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

hs_expt <- create_expt("sample_sheets/tmrc3_samples_20201105.xlsx",
                       file_column="hg3891hisatfile", savefile="hs_expt_all.rda",
                       gene_info=hs_annot)

## Reading the sample metadata.

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

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

## Reading count tables.

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

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30001/outputs/hisat2_hg38_91/forward.count.xz contains 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30002/outputs/hisat2_hg38_91/forward.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30003/outputs/hisat2_hg38_91/forward.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30004/outputs/hisat2_hg38_91/forward.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30005/outputs/hisat2_hg38_91/forward.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30006/outputs/hisat2_hg38_91/forward.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30007/outputs/hisat2_hg38_91/forward.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30009/outputs/hisat2_hg38_91/forward.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30010/outputs/hisat2_hg38_91/forward.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30015/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30011/outputs/hisat2_hg38_91/forward.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30012/outputs/hisat2_hg38_91/forward.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30013/outputs/hisat2_hg38_91/forward.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30016/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30017/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30050/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30052/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30071/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30056/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30058/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30018/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30019/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30014/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30021/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30029/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30020/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30038/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30039/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30023/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30025/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30022/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30046/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30047/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30048/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30026/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30030/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30031/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30032/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30024/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30040/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30033/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30049/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30053/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30054/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30037/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30027/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30028/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30034/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30035/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30036/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30044/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30055/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30068/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30070/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30041/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30042/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30043/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30045/outputs/hisat2_hg38_91/concatenated.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30059/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30060/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30061/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30062/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30063/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30051/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30064/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30065/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30066/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30067/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30057/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30069/outputs/hisat2_hg38_91/r1_trimmed.count_hg38_91_sno_gene_gene_id.count.xz contains 58307 rows and merges to 58307 rows.

## Finished reading count data.

## Matched 57768 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 58302 rows and 70 columns.

libsizes <- plot_libsize(hs_expt)

## The scale difference between the smallest and largest
## libraries is > 10. Assuming a log10 scale is better, set scale=FALSE if not.

libsizes$plot

nonzero <- plot_nonzero(hs_expt)
nonzero$plot

box <- plot_boxplot(hs_expt)

## This data will benefit from being displayed on the log scale.

## If this is not desired, set scale='raw'

## Some entries are 0.  We are on log scale, adding 1 to the data.

## Changed 2638993 zero count features.

box

3.2 Minimum coverage sample filtering

I arbitrarily chose 3,000,000 counts as a minimal level of coverage. We may want this to be higher.

hs_valid <- subset_expt(hs_expt, coverage=3000000)

## Subsetting given a minimal number of counts/sample.

## There were 70, now there are 64 samples.

plot_libsize(hs_valid)$plot

## The scale difference between the smallest and largest
## libraries is > 10. Assuming a log10 scale is better, set scale=FALSE if not.

3.3 Preparation

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

3.4 Global view

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

4 Macrophages

These samples are rather different from all of the others. The following section is therefore written primarily in response to a separate set of emails from Olga and Maria Adelaida; here is a snippet:

Dear all, about the samples corresponding to infected macrophages with three sensitive (2.2) and three resistant (2.3) clinical strains of L. (V.) panamensis, I send you the results of parasite burden by detection of 7SLRNA. I think these results are interesting, but the sample size is very small. Doctor Najib or Trey could you please send me the quality data and PCA analysis of these samples?

and

Hi Doctor, thank you. These samples corresponding to primary macrophages infected with clinical strains 2.2 (n=3) and 2.3 (n = 3). These information is in the file: TMRC project 3: excel Host TMRC3 v1.1, rows 137 to 150.

Thus I added 3 columns to the end of the sample sheet which seek to include this information. The first is ‘drug’ and encodes both the infection state (no for the two controls and yes for everything else), the second is zymodeme which I took from the tmrc2 sample sheet, the last is drug, which is either no or sb.

macr <- subset_expt(hs_expt, subset="typeofcells=='Macrophages'")

## Using a subset expression.

## There were 70, now there are 12 samples.

macr <- set_expt_conditions(macr, fact="zymodeme")
macr <- set_expt_batches(macr, fact="macrdrug")

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

## This function will replace the expt$expressionset slot with:

## log2(cpm(quant(cbcb(data))))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Not correcting the count-data for batch effects.  If batch is
##  included in EdgerR/limma's model, then this is probably wise; but in extreme
##  batch effects this is a good parameter to play with.

## Step 1: performing count filter with option: cbcb

## Removing 46807 low-count genes (11495 remaining).

## Step 2: normalizing the data with quant.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## Step 5: not doing batch correction.

norm_pca <- plot_pca(macr_norm)
pp(file="macrophage_side_experiment/norm_pca.png", image=norm_pca$plot)

## Writing the image to: macrophage_side_experiment/norm_pca.png and calling dev.off().

plot_3d_pca(norm_pca)

## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

## $plot
## 
## $file
## [1] "3dpca.html"

macr_nb <- normalize_expt(macr, batch="svaseq", filter=TRUE)

## This function will replace the expt$expressionset slot with:

## svaseq(cbcb(data))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Leaving the data in its current base format, keep in mind that
##  some metrics are easier to see when the data is log2 transformed, but
##  EdgeR/DESeq do not accept transformed data.

## Leaving the data unconverted.  It is often advisable to cpm/rpkm
##  the data to normalize for sampling differences, keep in mind though that rpkm
##  has some annoying biases, and voom() by default does a cpm (though hpgl_voom()
##  will try to detect this).

## Leaving the data unnormalized.  This is necessary for DESeq, but
##  EdgeR/limma might benefit from normalization.  Good choices include quantile,
##  size-factor, tmm, etc.

## Step 1: performing count filter with option: cbcb

## Removing 46807 low-count genes (11495 remaining).

## Step 2: not normalizing the data.

## Step 3: not converting the data.

## Step 4: not transforming the data.

## Step 5: doing batch correction with svaseq.

## Using the current state of normalization.

## Passing the data to all_adjusters using the svaseq estimate type.

## batch_counts: Before batch/surrogate estimation, 137726 entries are x>1: 100%.

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

## The be method chose 2 surrogate variables.

## Attempting svaseq estimation with 2 surrogates.

## There are 33 (0%) elements which are < 0 after batch correction.

## Setting low elements to zero.

macr_nb <- normalize_expt(macr_nb, norm="quant", convert="cpm", transform="log2")

## This function will replace the expt$expressionset slot with:

## log2(cpm(quant(data)))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Filter is false, this should likely be set to something, good
##  choices include cbcb, kofa, pofa (anything but FALSE).  If you want this to
##  stay FALSE, keep in mind that if other normalizations are performed, then the
##  resulting libsizes are likely to be strange (potentially negative!)

## Not correcting the count-data for batch effects.  If batch is
##  included in EdgerR/limma's model, then this is probably wise; but in extreme
##  batch effects this is a good parameter to play with.

## Step 1: not doing count filtering.

## Step 2: normalizing the data with quant.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## Step 5: not doing batch correction.

nb_pca <- plot_pca(macr_nb)
pp(file="macrophage_side_experiment/normbatch_pca.png", image=nb_pca$plot)

## Writing the image to: macrophage_side_experiment/normbatch_pca.png and calling dev.off().

macr_written <- write_expt(macr, excel="macrophage_side_experiment/macrophage_expt.xlsx")

## Writing the first sheet, containing a legend and some summary data.

## Writing the raw reads.

## Graphing the raw reads.

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete

## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete

## Attempting mixed linear model with: ~  condition + batch

## Fitting the expressionset to the model, this is slow.

## Loading required package: Matrix

## 
## Total:160 s

## Placing factor: condition at the beginning of the model.

## Writing the normalized reads.

## Graphing the normalized reads.

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete

## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete

## Attempting mixed linear model with: ~  condition + batch

## Fitting the expressionset to the model, this is slow.

## 
## Total:67 s

## Placing factor: condition at the beginning of the model.

## Writing the median reads by factor.

##tmp <- normalize_expt(macr, filter=TRUE)
##zy_de <- deseq_pairwise(tmp, model_batch="svaseq")
##zy_ed <- edger_pairwise(tmp, model_batch="svaseq")
zymo_de <- all_pairwise(macr, model_batch="svaseq", filter=TRUE, parallel=FALSE)

## batch_counts: Before batch/surrogate estimation, 137726 entries are x>1: 100%.

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

## The be method chose 2 surrogate variables.

## Attempting svaseq estimation with 2 surrogates.

## Plotting a PCA before surrogate/batch inclusion.

## Not putting labels on the PC plot.

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

## Performing a test normalization with: raw

## This function will replace the expt$expressionset slot with:

## log2(svaseq(cpm(cbcb(data))))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Leaving the data unnormalized.  This is necessary for DESeq, but
##  EdgeR/limma might benefit from normalization.  Good choices include quantile,
##  size-factor, tmm, etc.

## Step 1: performing count filter with option: cbcb

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

## Step 2: not normalizing the data.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

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

## Step 5: doing batch correction with svaseq.

## Using the current state of normalization.

## Passing the data to all_adjusters using the svaseq estimate type.

## batch_counts: Before batch/surrogate estimation, 133046 entries are x>1: 96%.

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

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

## The be method chose 2 surrogate variables.

## Attempting svaseq estimation with 2 surrogates.

## There are 238 (0%) elements which are < 0 after batch correction.

## Setting low elements to zero.

## Not putting labels on the PC plot.

## Starting basic_pairwise().

## Starting basic pairwise comparison.

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## Basic step 0/3: Transforming data.

## Basic step 1/3: Creating mean and variance tables.

## Basic step 2/3: Performing 6 comparisons.

## Basic step 3/3: Creating faux DE Tables.

## Basic: Returning tables.

## Starting deseq_pairwise().

## Starting DESeq2 pairwise comparisons.

## The data should be suitable for EdgeR/DESeq/EBSeq. If they freak out, check the state of the count table and ensure that it is in integer counts.

## Including batch estimates from sva/ruv/pca in the model.

## Choosing the non-intercept containing model.

## DESeq2 step 1/5: Including a matrix of batch estimates in the deseq model.

## converting counts to integer mode

## DESeq2 step 2/5: Estimate size factors.

## DESeq2 step 3/5: Estimate dispersions.

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## Using a parametric fitting seems to have worked.

## DESeq2 step 4/5: nbinomWaldTest.

## Starting ebseq_pairwise().

## The data should be suitable for EdgeR/DESeq/EBSeq. If they freak out, check the state of the count table and ensure that it is in integer counts.

## Starting EBSeq pairwise subset.

## Choosing the non-intercept containing model.

## Starting EBTest of undef vs. z22.

## Copying ppee values as ajusted p-values until I figure out how to deal with them.

## Starting EBTest of undef vs. z23.

## Copying ppee values as ajusted p-values until I figure out how to deal with them.

## Starting EBTest of z22 vs. z23.

## Copying ppee values as ajusted p-values until I figure out how to deal with them.

## Starting edger_pairwise().

## Starting edgeR pairwise comparisons.

## The data should be suitable for EdgeR/DESeq/EBSeq. If they freak out, check the state of the count table and ensure that it is in integer counts.

## Including batch estimates from sva/ruv/pca in the model.

## Choosing the non-intercept containing model.

## EdgeR step 1/9: Importing and normalizing data.

## EdgeR step 2/9: Estimating the common dispersion.

## EdgeR step 3/9: Estimating dispersion across genes.

## EdgeR step 4/9: Estimating GLM Common dispersion.

## EdgeR step 5/9: Estimating GLM Trended dispersion.

## EdgeR step 6/9: Estimating GLM Tagged dispersion.

## EdgeR step 7/9: Running glmFit, switch to glmQLFit by changing the argument 'edger_test'.

## EdgeR step 8/9: Making pairwise contrasts.

## Starting limma_pairwise().

## Starting limma pairwise comparison.

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$best_libsize.

## Limma step 1/6: choosing model.

## Including batch estimates from sva/ruv/pca in the model.

## Choosing the non-intercept containing model.

## Limma step 2/6: running limma::voom(), switch with the argument 'which_voom'.

## Using normalize.method=quantile for 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/3: Creating table: z22_vs_undef.  Adjust=BH

## Limma step 6/6: 2/3: Creating table: z23_vs_undef.  Adjust=BH

## Limma step 6/6: 3/3: Creating table: z23_vs_z22.  Adjust=BH

## Limma step 6/6: 1/3: Creating table: undef.  Adjust=BH

## Limma step 6/6: 2/3: Creating table: z22.  Adjust=BH

## Limma step 6/6: 3/3: Creating table: z23.  Adjust=BH

## Comparing analyses.

zymo_table <- combine_de_tables(zymo_de, excel="macrophage_side_experiment/macrophage_de.xlsx", parallel=FALSE)

## Deleting the file macrophage_side_experiment/macrophage_de.xlsx before writing the tables.

## Writing a legend of columns.

## Printing a pca plot before/after surrogates/batch estimation.

## Working on table 1/3: z22_vs_undef

## Working on table 2/3: z23_vs_undef

## Working on table 3/3: z23_vs_z22

## Adding venn plots for z22_vs_undef.

## Limma expression coefficients for z22_vs_undef; R^2: 0.963; equation: y = 0.973x + 0.146

## Deseq expression coefficients for z22_vs_undef; R^2: 0.962; equation: y = 1.02x - 0.178

## Edger expression coefficients for z22_vs_undef; R^2: 0.962; equation: y = 1.02x - 0.159

## Adding venn plots for z23_vs_undef.

## Limma expression coefficients for z23_vs_undef; R^2: 0.956; equation: y = 0.976x + 0.152

## Deseq expression coefficients for z23_vs_undef; R^2: 0.953; equation: y = 0.971x + 0.307

## Edger expression coefficients for z23_vs_undef; R^2: 0.953; equation: y = 0.971x + 0.298

## Adding venn plots for z23_vs_z22.

## Limma expression coefficients for z23_vs_z22; R^2: 0.969; equation: y = 0.99x + 0.0703

## Deseq expression coefficients for z23_vs_z22; R^2: 0.966; equation: y = 0.946x + 0.587

## Edger expression coefficients for z23_vs_z22; R^2: 0.966; equation: y = 0.946x + 0.515

## Writing summary information, compare_plot is: TRUE.

## Performing save of macrophage_side_experiment/macrophage_de.xlsx.

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

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

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

## > git reset d6114393775c44346befc5916efed02675b9bd7f

## This is hpgltools commit: Fri Nov 6 15:01:37 2020 -0500: d6114393775c44346befc5916efed02675b9bd7f

## Saving to 01_annotation_v202011.rda.xz

tmp <- loadme(filename=savefile)