L. panamensis 20191001: Sample Estimation

1 Sample Estimation

1.1 Generate expressionsets

Caveat: This is using hisat2 quantifications.

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

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

## I think samples 7,10 should be removed at minimum, probably also 9,11
nonzero <- plot_nonzero(hs_expt)
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 23422 zero count features.

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

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

## Writing the raw reads.

## Graphing the raw reads.

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

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

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

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

## Writing the normalized reads.

## Error in density.default(x, adjust = adj) : 'x' contains missing values
## Error in density.default(x, adjust = adj) : 'x' contains missing values

## Graphing the normalized reads.

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

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

## Writing the median reads by factor.

hs_valid <- subset_expt(hs_expt, coverage=100000)

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

## There were 13, now there are 11 samples.

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

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

## Writing the raw reads.

## Graphing the raw reads.

## Warning: Removed 3 rows containing missing values (geom_text).

## Warning: Removed 3 rows containing missing values (geom_text).

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

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

## Writing the normalized reads.

## Error in density.default(x, adjust = adj) : 'x' contains missing values
## Error in density.default(x, adjust = adj) : 'x' contains missing values

## Graphing the normalized reads.

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

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

## Writing the median reads by factor.

hs_expt <- create_expt("sample_sheets/tmrc3_samples_20191001.xlsx",
                       file_column="hg3891hisatfile",
                       gene_info=hs_annot)

## Reading the sample metadata.

## The sample definitions comprises: 17 rows(samples) and 62 columns(metadata fields).

## Reading count tables.

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

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

## /mnt/sshfs/cbcbsub/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/cbcbsub/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/cbcbsub/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/cbcbsub/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/cbcbsub/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/cbcbsub/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/cbcbsub/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc30008/outputs/hisat2_hg38_91/forward.count.xz contains 58307 rows and merges to 58307 rows.

## /mnt/sshfs/cbcbsub/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/cbcbsub/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/cbcbsub/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/cbcbsub/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/cbcbsub/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.

## Finished reading count tables.

## Matched 58243 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 13 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

## I think samples 7,10 should be removed at minimum, probably also 9,11
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 494937 zero count features.

box

## This is causing segmentation faults in R
##hs_write <- write_expt(hs_expt, excel=glue("excel/hs_hisat2_written-v{ver}.xlsx"))

hs_valid <- subset_expt(hs_expt, coverage=3000000)

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

## There were 13, now there are 11 samples.

plot_libsize(hs_valid)$plot

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

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

## Writing the raw reads.

## Graphing the raw reads.

## Warning: Removed 3 rows containing missing values (geom_text).

## Warning: Removed 3 rows containing missing values (geom_text).

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

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

## Writing the normalized reads.

## Graphing the normalized reads.

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

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

## Writing the median reads by factor.

2 Questions from Maria Adelaida

The following comes from an email 20190830 from Maria Adelaida.

1. Samples WT1010 and WT1011 PBMCs from two healthy donors processed 2h, 7h and 12h after sample procurement. This is an analysis to explore the time-effect on gene expression and define steps for data analysis for patient samples considering time-dependent effects.

   1. An initial PCA on the raw data would be very useful to see if there is clustering based on time or (as usual), mostly a donor-specific effect. Then I think a hierarchical clustering of genes based on time-dependent modifications to see what is mostly affected (if any) - like what you guys did for T.cruzi.

2. Samples from SU1017, SU1034 Samples from TMRC CL patients. m= monocyte, n= neutrophil. Samples labeled “1” are taken before treatment and those “2” mid way through treatment. This is exiting, because these will be our first neutrophil transcriptomes.

In an attempt to poke at these questions, I mapped the reads to hg38_91 using salmon and hisat2. It is very noteworthy that the salmon mappings are exhibiting some serious problems and should be looked into further. The hisat2 mappings are significantly more ‘normal’. Having said that, two samples remain basically unusable: tmrc30009 (1034n1) and (to a smaller degree) tmrc30007 (1017n1) have too few reads as shown above.

2.1 Preparation

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

2.2 Global view

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

hs_valid <- set_expt_batches(hs_valid, fact="donor")
all_norm <- normalize_expt(hs_valid, norm="quant", transform="log2", convert="cpm", batch=FALSE,
                           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 43798 low-count genes (14504 remaining).

## Step 2: normalizing the data with quant.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

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

## Step 5: not doing batch correction.

all_pca <- plot_pca(all_norm)
knitr::kable(all_pca$table)

	sampleid	condition	batch	batch_int	colors	labels	PC1	PC2	pc_1	pc_2	pc_3	pc_4	pc_5	pc_6	pc_7	pc_8	pc_9	pc_10
tmrc30001	TMRC30001	healthy_2hr	d1010	1	#1B9E77	TMRC30001	-0.2057	0.0785	-0.2057	0.0785	0.0752	0.6722	-0.0517	0.0514	-0.3254	0.2793	-0.0402	0.4607
tmrc30002	TMRC30002	healthy_7hr	d1010	1	#D95F02	TMRC30002	-0.2087	0.1218	-0.2087	0.1218	-0.4252	0.2430	-0.0752	-0.0399	0.7361	-0.2157	-0.1169	0.0400
tmrc30003	TMRC30003	healthy_12hr	d1010	1	#7570B3	TMRC30003	-0.2305	0.2114	-0.2305	0.2114	-0.5807	-0.0727	0.0416	0.0496	-0.4850	-0.0226	0.1595	-0.4510
tmrc30004	TMRC30004	healthy_2hr	d1011	2	#1B9E77	TMRC30004	-0.2060	0.0790	-0.2060	0.0790	0.5984	0.1779	-0.0955	-0.0400	0.0991	-0.2649	0.4738	-0.3943
tmrc30005	TMRC30005	healthy_7hr	d1011	2	#D95F02	TMRC30005	-0.2494	0.2307	-0.2494	0.2307	0.3347	-0.2517	0.0114	0.0054	-0.0267	0.1245	-0.7510	-0.1945
tmrc30006	TMRC30006	healthy_12hr	d1011	2	#7570B3	TMRC30006	-0.2290	0.2256	-0.2290	0.2256	0.0337	-0.6042	0.0076	-0.0105	0.0326	0.0510	0.3456	0.5624
tmrc30008	TMRC30008	monocyte_early	d1017	3	#66A61E	TMRC30008	0.0039	-0.4277	0.0039	-0.4277	0.0147	0.0058	0.8215	-0.2164	0.0431	0.0449	0.0131	-0.0129
tmrc30009	TMRC30009	neutrophil_early	d1034	4	#E7298A	TMRC30009	0.5796	0.2219	0.5796	0.2219	-0.0215	0.0132	-0.1785	-0.6973	-0.0659	0.0212	-0.0203	0.0045
tmrc30011	TMRC30011	neutrophil_late	d1034	4	#E6AB02	TMRC30011	0.5975	0.2746	0.5975	0.2746	0.0419	0.0260	0.1989	0.6550	0.0663	-0.0305	0.0188	-0.0060
tmrc30012	TMRC30012	monocyte_late	d1034	4	#A6761D	TMRC30012	0.0739	-0.5233	0.0739	-0.5233	-0.0284	-0.0839	-0.2995	0.1131	-0.2517	-0.6194	-0.1945	0.1863
tmrc30013	TMRC30013	monocyte_late	d1034	4	#A6761D	TMRC30013	0.0743	-0.4925	0.0743	-0.4925	-0.0430	-0.1256	-0.3804	0.1296	0.1775	0.6322	0.1123	-0.1953

plot_corheat(all_norm)$plot

plot_topn(hs_valid)$plot

2.3 Question 1

I interpreted question 1 as: pull out tmrc3000[1-6] and look at them.

I am not entirely certain what is meant by the heirarchical clustering request. I can see a couple of possibilities for what this means. The most similar thing I recall in the cruzi context was a post-DE visualization of some fairly specific genes.

hs_q1 <- subset_expt(hs_valid, subset="donor=='d1010'|donor=='d1011'")

## Using a subset expression.

## There were 11, now there are 6 samples.

q1_norm <- normalize_expt(hs_q1, norm="quant", transform="log2", convert="cpm", batch=FALSE,
                           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 44614 low-count genes (13688 remaining).

## Step 2: normalizing the data with quant.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

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

## Step 5: not doing batch correction.

q1_pca <- plot_pca(q1_norm)
q1_pca$plot

knitr::kable(q1_pca$table)

	sampleid	condition	batch	batch_int	colors	labels	PC1	PC2	pc_1	pc_2	pc_3	pc_4	pc_5
tmrc30001	TMRC30001	healthy_2hr	d1010	1	#1B9E77	TMRC30001	-0.5759	0.2704	-0.5759	0.2704	0.4004	-0.0566	0.5148
tmrc30002	TMRC30002	healthy_7hr	d1010	1	#D95F02	TMRC30002	0.0048	0.4828	0.0048	0.4828	-0.6648	0.3977	-0.0112
tmrc30003	TMRC30003	healthy_12hr	d1010	1	#7570B3	TMRC30003	0.4107	0.4514	0.4107	0.4514	0.4033	-0.3107	-0.4492
tmrc30004	TMRC30004	healthy_2hr	d1011	2	#1B9E77	TMRC30004	-0.4764	-0.4166	-0.4764	-0.4166	-0.2714	-0.3745	-0.4679
tmrc30005	TMRC30005	healthy_7hr	d1011	2	#D95F02	TMRC30005	0.1313	-0.4486	0.1313	-0.4486	0.3425	0.6930	-0.1315
tmrc30006	TMRC30006	healthy_12hr	d1011	2	#7570B3	TMRC30006	0.5054	-0.3394	0.5054	-0.3394	-0.2100	-0.3488	0.5449

## Looks like PC1 == Time and PC2 is donor, and they have pretty much the same amount of variance

hs_time <- set_expt_conditions(hs_q1, fact="time")
time_norm <- normalize_expt(hs_time, transform="log2",
                            batch="svaseq", filter=TRUE)

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

## log2(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 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 44614 low-count genes (13688 remaining).

## Step 2: not normalizing the data.

## Step 3: not converting the data.

## Step 4: transforming the data with log2.

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

## Step 5: doing batch correction with svaseq.

## Note to self:  If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.

## Passing off to all_adjusters.

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

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

## The be method chose 1 surrogate variable(s).

## Attempting svaseq estimation with 1 surrogates.

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

## Setting low elements to zero.

## Get a set of genes with high PC loads for PC1(time), and PC2(donor):
high_scores <- pca_highscores(time_norm, n=40)

time_stuff <- high_scores$highest[, c(1, 2)]
time_stuff

##       Comp.1                  Comp.2                 
##  [1,] "11.95:ENSG00000073756" "11.76:ENSG00000242041"
##  [2,] "11.2:ENSG00000276107"  "11.28:ENSG00000100298"
##  [3,] "10.4:ENSG00000169245"  "10.98:ENSG00000223609"
##  [4,] "10.11:ENSG00000239948" "9.871:ENSG00000211803"
##  [5,] "9.636:ENSG00000125538" "9.823:ENSG00000278030"
##  [6,] "9.048:ENSG00000235576" "8.979:ENSG00000178429"
##  [7,] "8.984:ENSG00000093134" "8.701:ENSG00000106341"
##  [8,] "8.747:ENSG00000213876" "8.684:ENSG00000224376"
##  [9,] "8.638:ENSG00000112299" "8.345:ENSG00000185499"
## [10,] "7.741:ENSG00000267731" "8.165:ENSG00000011426"
## [11,] "7.646:ENSG00000174946" "7.765:ENSG00000026559"
## [12,] "7.637:ENSG00000123700" "7.61:ENSG00000167874" 
## [13,] "7.596:ENSG00000012817" "7.437:ENSG00000012817"
## [14,] "7.575:ENSG00000274536" "7.426:ENSG00000129824"
## [15,] "7.56:ENSG00000168329"  "7.42:ENSG00000227081" 
## [16,] "7.333:ENSG00000171049" "7.042:ENSG00000051180"
## [17,] "7.225:ENSG00000245954" "6.946:ENSG00000128641"
## [18,] "7.202:ENSG00000140092" "6.703:ENSG00000077152"
## [19,] "7.155:ENSG00000187474" "6.637:ENSG00000196169"
## [20,] "7.065:ENSG00000126243" "6.609:ENSG00000125538"
## [21,] "7.043:ENSG00000086570" "6.539:ENSG00000217130"
## [22,] "6.995:ENSG00000206503" "6.526:ENSG00000206077"
## [23,] "6.818:ENSG00000277496" "6.487:ENSG00000266402"
## [24,] "6.806:ENSG00000173114" "6.376:ENSG00000243716"
## [25,] "6.789:ENSG00000170915" "6.317:ENSG00000144152"
## [26,] "6.786:ENSG00000125703" "6.214:ENSG00000108556"
## [27,] "6.708:ENSG00000211817" "6.155:ENSG00000167191"
## [28,] "6.449:ENSG00000261799" "6.121:ENSG00000206503"
## [29,] "6.25:ENSG00000263002"  "6.111:ENSG00000186354"
## [30,] "6.239:ENSG00000091106" "6.055:ENSG00000235618"
## [31,] "6.217:ENSG00000162739" "6.039:ENSG00000169245"
## [32,] "6.132:ENSG00000143365" "6.013:ENSG00000104951"
## [33,] "6.096:ENSG00000204469" "5.895:ENSG00000154736"
## [34,] "6.031:ENSG00000211793" "5.827:ENSG00000164877"
## [35,] "6.007:ENSG00000120526" "5.785:ENSG00000067048"
## [36,] "6:ENSG00000164484"     "5.777:ENSG00000114374"
## [37,] "5.997:ENSG00000268001" "5.763:ENSG00000114315"
## [38,] "5.978:ENSG00000168310" "5.69:ENSG00000131002" 
## [39,] "5.957:ENSG00000120262" "5.601:ENSG00000162383"
## [40,] "5.954:ENSG00000272632" "5.528:ENSG00000228140"

## Column 1 should be winners vs. time and column 2 by donor.
time_genes <- gsub(x=time_stuff[, "Comp.1"], pattern="^.*:(.*)$", replacement="\\1")
head(time_genes)

## [1] "ENSG00000073756" "ENSG00000276107" "ENSG00000169245" "ENSG00000239948"
## [5] "ENSG00000125538" "ENSG00000235576"

time_subset <- exprs(time_norm)[time_genes, ]
plot_sample_heatmap(time_norm, row_label=rownames(time_norm))

hs_donor <- set_expt_conditions(hs_q1, fact="donor")
donor_norm <- normalize_expt(hs_donor, convert="cpm", transform="log2",
                            batch="svaseq", filter=TRUE)

## 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 44614 low-count genes (13688 remaining).

## Step 2: not normalizing the data.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

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

## Step 5: doing batch correction with svaseq.

## Note to self:  If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.

## Passing off to all_adjusters.

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

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

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

## The be method chose 1 surrogate variable(s).

## Attempting svaseq estimation with 1 surrogates.

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

## Setting low elements to zero.

## Get a set of genes with high PC loads for PC1(donor), and PC2(donor):
high_scores <- pca_highscores(donor_norm, n=40)

donor_stuff <- high_scores$highest[, c(1, 2)]
donor_stuff

##       Comp.1                  Comp.2                 
##  [1,] "37.81:ENSG00000129824" "7.091:ENSG00000230291"
##  [2,] "36.01:ENSG00000067048" "6.637:ENSG00000227081"
##  [3,] "34.87:ENSG00000012817" "5.616:ENSG00000217130"
##  [4,] "30.89:ENSG00000179344" "5.47:ENSG00000235174" 
##  [5,] "30.13:ENSG00000114374" "5.365:ENSG00000230897"
##  [6,] "29.9:ENSG00000206503"  "5.023:ENSG00000210049"
##  [7,] "29.47:ENSG00000131002" "4.946:ENSG00000266402"
##  [8,] "28.5:ENSG00000196735"  "4.319:ENSG00000256937"
##  [9,] "27.6:ENSG00000239830"  "4.023:ENSG00000129824"
## [10,] "26.2:ENSG00000183878"  "3.919:ENSG00000206077"
## [11,] "26.09:ENSG00000234745" "3.769:ENSG00000245748"
## [12,] "25.51:ENSG00000099725" "3.7:ENSG00000225616"  
## [13,] "23.88:ENSG00000213058" "3.678:ENSG00000242041"
## [14,] "23.44:ENSG00000067646" "3.663:ENSG00000144152"
## [15,] "23.36:ENSG00000233864" "3.639:ENSG00000108556"
## [16,] "19.63:ENSG00000198692" "3.541:ENSG00000240616"
## [17,] "19.28:ENSG00000259045" "3.472:ENSG00000243716"
## [18,] "17.79:ENSG00000204642" "3.465:ENSG00000112212"
## [19,] "16.05:ENSG00000215580" "3.421:ENSG00000037280"
## [20,] "14.9:ENSG00000112139"  "3.418:ENSG00000172602"
## [21,] "13.58:ENSG00000271207" "3.417:ENSG00000105856"
## [22,] "12.56:ENSG00000049247" "3.414:ENSG00000145850"
## [23,] "12.37:ENSG00000223534" "3.364:ENSG00000224376"
## [24,] "12.18:ENSG00000115602" "3.351:ENSG00000159842"
## [25,] "11.57:ENSG00000131042" "3.234:ENSG00000235618"
## [26,] "11.57:ENSG00000129422" "3.224:ENSG00000278993"
## [27,] "11.54:ENSG00000228314" "3.208:ENSG00000228887"
## [28,] "11.49:ENSG00000224114" "3.195:ENSG00000179344"
## [29,] "11.17:ENSG00000127022" "3.187:ENSG00000143786"
## [30,] "10.75:ENSG00000254681" "3.182:ENSG00000213638"
## [31,] "10.44:ENSG00000144115" "3.144:ENSG00000211797"
## [32,] "10.34:ENSG00000196436" "3.144:ENSG00000061656"
## [33,] "10.13:ENSG00000230327" "3.137:ENSG00000204620"
## [34,] "10.06:ENSG00000196126" "3.102:ENSG00000095397"
## [35,] "9.997:ENSG00000154620" "3.1:ENSG00000239532"  
## [36,] "9.851:ENSG00000225630" "3.056:ENSG00000099725"
## [37,] "9.56:ENSG00000241945"  "3.021:ENSG00000214756"
## [38,] "9.31:ENSG00000109321"  "3.013:ENSG00000258315"
## [39,] "8.724:ENSG00000182534" "3.012:ENSG00000230202"
## [40,] "8.413:ENSG00000204001" "3.001:ENSG00000158966"

## Column 1 should be winners vs. donor and column 2 by donor.
donor_genes <- gsub(x=donor_stuff[, "Comp.1"], pattern="^.*:(.*)$", replacement="\\1")
head(donor_genes)

## [1] "ENSG00000129824" "ENSG00000067048" "ENSG00000012817" "ENSG00000179344"
## [5] "ENSG00000114374" "ENSG00000206503"

donor_subset <- exprs(donor_norm)[donor_genes, ]
plot_sample_heatmap(donor_norm, row_label=rownames(donor_norm))

time_keepers <- list(
  "2hr_vs_7hr" = c("t2hr", "t7hr"),
  "2hr_vs_12hr" = c("t2hr", "t12hr"),
  "7hr_vs_12hr" = c("t7hr", "t12hr"))

q1_time <- set_expt_conditions(hs_q1, fact="time")
time_de <- all_pairwise(q1_time, model_batch=FALSE, parallel=FALSE)

## 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 44614 low-count genes (13688 remaining).

## Step 2: normalizing the data with quant.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

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

## Step 5: not doing batch correction.

## Plotting a PCA before surrogates/batch inclusion.

## Not putting labels on the plot.

## Assuming no batch in model for testing pca.

## Not putting labels on the plot.

## Starting basic_pairwise().

## Starting basic pairwise comparison.

## Basic step 0/3: Filtering data.

## 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.

## Choosing the non-intercept containing model.

## DESeq2 step 1/5: Including only condition 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 t12hr vs. t2hr.

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

## Starting EBTest of t12hr vs. t7hr.

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

## Starting EBTest of t2hr vs. t7hr.

## 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.

## 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$libsize.

## Limma step 1/6: choosing 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: t2hr_vs_t12hr.  Adjust=BH

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

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

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

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

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

## Comparing analyses.

time_all_tables <- combine_de_tables(time_de, excel=glue::glue("excel/time_de_tables-v{ver}.xlsx"))

## Writing a legend of columns.

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

## Working on table 1/3: t2hr_vs_t12hr

## Working on table 2/3: t7hr_vs_t12hr

## Working on table 3/3: t7hr_vs_t2hr

## Adding venn plots for t2hr_vs_t12hr.

## Limma expression coefficients for t2hr_vs_t12hr; R^2: 0.967; equation: y = 0.982x - 0.0165

## Deseq expression coefficients for t2hr_vs_t12hr; R^2: 0.925; equation: y = 0.954x + 0.206

## Edger expression coefficients for t2hr_vs_t12hr; R^2: 0.94; equation: y = 0.972x + 0.0868

## Adding venn plots for t7hr_vs_t12hr.

## Limma expression coefficients for t7hr_vs_t12hr; R^2: 0.967; equation: y = 0.982x - 0.0165

## Deseq expression coefficients for t7hr_vs_t12hr; R^2: 0.925; equation: y = 0.954x + 0.206

## Edger expression coefficients for t7hr_vs_t12hr; R^2: 0.94; equation: y = 0.972x + 0.0868

## Adding venn plots for t7hr_vs_t2hr.

## Limma expression coefficients for t7hr_vs_t2hr; R^2: 0.967; equation: y = 0.982x - 0.0165

## Deseq expression coefficients for t7hr_vs_t2hr; R^2: 0.925; equation: y = 0.954x + 0.206

## Edger expression coefficients for t7hr_vs_t2hr; R^2: 0.94; equation: y = 0.972x + 0.0868

## Writing summary information, compare_plot is: TRUE.

## Performing save of excel/time_de_tables-v20191001.xlsx.

time_all_tables_all <- combine_de_tables(time_de,
                                         keepers="all",
                                         excel=glue::glue("excel/time_de_all_tables-v{ver}.xlsx"))

## Deleting the file excel/time_de_all_tables-v20191001.xlsx before writing the tables.

## Writing a legend of columns.

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

## Working on table 1/3: t2hr_vs_t12hr

## Working on table 2/3: t7hr_vs_t12hr

## Working on table 3/3: t7hr_vs_t2hr

## Adding venn plots for t2hr_vs_t12hr.

## Limma expression coefficients for t2hr_vs_t12hr; R^2: 0.967; equation: y = 0.982x - 0.0165

## Deseq expression coefficients for t2hr_vs_t12hr; R^2: 0.925; equation: y = 0.954x + 0.206

## Edger expression coefficients for t2hr_vs_t12hr; R^2: 0.94; equation: y = 0.972x + 0.0868

## Adding venn plots for t7hr_vs_t12hr.

## Limma expression coefficients for t7hr_vs_t12hr; R^2: 0.967; equation: y = 0.982x - 0.0165

## Deseq expression coefficients for t7hr_vs_t12hr; R^2: 0.925; equation: y = 0.954x + 0.206

## Edger expression coefficients for t7hr_vs_t12hr; R^2: 0.94; equation: y = 0.972x + 0.0868

## Adding venn plots for t7hr_vs_t2hr.

## Limma expression coefficients for t7hr_vs_t2hr; R^2: 0.967; equation: y = 0.982x - 0.0165

## Deseq expression coefficients for t7hr_vs_t2hr; R^2: 0.925; equation: y = 0.954x + 0.206

## Edger expression coefficients for t7hr_vs_t2hr; R^2: 0.94; equation: y = 0.972x + 0.0868

## Writing summary information, compare_plot is: TRUE.

## Performing save of excel/time_de_all_tables-v20191001.xlsx.

time_all_tables <- combine_de_tables(time_de,
                                     keepers=time_keepers,
                                     excel=glue::glue("excel/{rundate}-time_de_tables-v{ver}.xlsx"))

## Writing a legend of columns.

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

## Working on 1/3: 2hr_vs_7hr which is: t2hr/t7hr.

## Found inverse table with t7hr_vs_t2hr

## Working on 2/3: 2hr_vs_12hr which is: t2hr/t12hr.

## Found table with t2hr_vs_t12hr

## Working on 3/3: 7hr_vs_12hr which is: t7hr/t12hr.

## Found table with t7hr_vs_t12hr

## Adding venn plots for 2hr_vs_7hr.

## Limma expression coefficients for 2hr_vs_7hr; R^2: 0.967; equation: y = 0.982x - 0.0165

## Deseq expression coefficients for 2hr_vs_7hr; R^2: 0.925; equation: y = 0.954x + 0.206

## Edger expression coefficients for 2hr_vs_7hr; R^2: 0.94; equation: y = 0.972x + 0.0868

## Adding venn plots for 2hr_vs_12hr.

## Limma expression coefficients for 2hr_vs_12hr; R^2: 0.967; equation: y = 0.982x - 0.0165

## Deseq expression coefficients for 2hr_vs_12hr; R^2: 0.925; equation: y = 0.954x + 0.206

## Edger expression coefficients for 2hr_vs_12hr; R^2: 0.94; equation: y = 0.972x + 0.0868

## Adding venn plots for 7hr_vs_12hr.

## Limma expression coefficients for 7hr_vs_12hr; R^2: 0.967; equation: y = 0.982x - 0.0165

## Deseq expression coefficients for 7hr_vs_12hr; R^2: 0.925; equation: y = 0.954x + 0.206

## Edger expression coefficients for 7hr_vs_12hr; R^2: 0.94; equation: y = 0.972x + 0.0868

## Writing summary information, compare_plot is: TRUE.

## Performing save of excel/20191211-time_de_tables-v20191001.xlsx.

q1_donor <- set_expt_conditions(hs_q1, fact="donor")
donor_de <- all_pairwise(q1_donor, model_batch=FALSE)

## 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 44614 low-count genes (13688 remaining).

## Step 2: normalizing the data with quant.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

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

## Step 5: not doing batch correction.

## Plotting a PCA before surrogates/batch inclusion.

## Not putting labels on the plot.

## Assuming no batch in model for testing pca.

## Not putting labels on the plot.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

donor_tables <- combine_de_tables(donor_de, excel=glue::glue("excel/donor_de_tables-v{ver}.xlsx"))

## Writing a legend of columns.

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

## Working on table 1/1: d1011_vs_d1010

## Adding venn plots for d1011_vs_d1010.

## Limma expression coefficients for d1011_vs_d1010; R^2: 0.975; equation: y = 0.986x + 0.00062

## Deseq expression coefficients for d1011_vs_d1010; R^2: 0.954; equation: y = 0.989x + 0.0249

## Edger expression coefficients for d1011_vs_d1010; R^2: 1; equation: y = 1x - 1.72e-09

## Writing summary information, compare_plot is: TRUE.

## Performing save of excel/donor_de_tables-v20191001.xlsx.

hs_q2 <- subset_expt(hs_valid, subset="donor!='d1010'&donor!='d1011'")

## Using a subset expression.

## There were 11, now there are 5 samples.

q2_norm <- normalize_expt(hs_q2, transform="log2", convert="cpm", norm="quant", 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 46075 low-count genes (12227 remaining).

## Step 2: normalizing the data with quant.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

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

## Step 5: not doing batch correction.

q2_pca <- plot_pca(q2_norm)
knitr::kable(q2_pca$table)

	sampleid	condition	batch	batch_int	colors	labels	PC1	PC2	pc_1	pc_2	pc_3	pc_4
tmrc30008	TMRC30008	monocyte_early	d1017	1	#66A61E	TMRC30008	-0.3841	-0.8033	-0.3841	-0.8033	0.0827	-0.0204
tmrc30009	TMRC30009	neutrophil_early	d1034	2	#E7298A	TMRC30009	0.5272	0.0713	0.5272	0.0713	0.7188	-0.0148
tmrc30011	TMRC30011	neutrophil_late	d1034	2	#E6AB02	TMRC30011	0.5674	-0.0929	0.5674	-0.0929	-0.6852	0.0017
tmrc30012	TMRC30012	monocyte_late	d1034	2	#A6761D	TMRC30012	-0.3630	0.4329	-0.3630	0.4329	-0.0691	-0.6899
tmrc30013	TMRC30013	monocyte_late	d1034	2	#A6761D	TMRC30013	-0.3475	0.3919	-0.3475	0.3919	-0.0472	0.7234

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 ca7dfe3c77056f495da05ef0c95b6208370a23f1commit!

## This is hpgltools commit: Sat Nov 23 19:41:15 2019 -0500: ca7dfe3c77056f495da05ef0c95b6208370a23f1This is hpgltools commit: Sat Nov 23 19:41:15 2019 -0500: commit!

## Saving to 01_annotation_v20191001.rda.xz

tmp <- loadme(filename=savefile)