RNAseq of L.panamensis: Human Macrophage Sample Estimation.

index.html annotation.html

This document is concerned with analyzing RNAseq data of human and parasite during the infection of human macrophages. A few different strains of L. panamensis were used; the experiment is therefore segregated into groups named ‘self-healing’, ‘chronic’, and ‘uninfected’. Two separate sets of libraries were generated, one earlier set with greater coverage, and a later set including only 1 uninfected sample, and 2-3 chronic samples.

0.1 Extract the macrophage experiment

The following subset operation extract the samples used for the macrophage experiment. The next three lines then change the colors from the defaults.

macrophage_human <- expt_subset(human_expt, subset="condition=='macro_nil'|condition=='macro_sh'|condition=='macro_ch'")
chosen_colors <- c("#009900","#990000", "#000099")
names(chosen_colors) <- c("macro_nil","macro_ch","macro_sh")
macrophage_human <- set_expt_colors(macrophage_human, colors=chosen_colors)

1 Sample metrics

The various metrics used to describe and examine the data come once before, and once after normalization.

1.1 Initial sample metrics

graph_metrics creates a set of plots we can use to look at the data before messing with it.

written <- sm(write_expt(macrophage_human,
                         excel=paste0("excel/macrophage_human_data-v", ver, ".xlsx"),
                         violin=TRUE))

macrophage_metrics <- sm(graph_metrics(macrophage_human))

The two most likely plots one will want to see immediately are the library sizes and data densities.

macrophage_metrics$libsize

## The libraries > 600 in HPGL id have lower coverage as they were sequenced many more to a lane.
macrophage_metrics$density

## This same phenomenon is visible in the density plot
macrophage_metrics$nonzero

## and the non-zero gene / cpm plot.

1.2 Normalization effects on the macrophage data

The plotted library sizes suggest that differences observed in the raw data are due at least partially to less saturation in the new samples. Therefore, low-count filtering is almost certainly to prove an important parameter. Similarly, a batch correction should help alleviate these differences.

nonil_mac <- expt_subset(macrophage_human, subset="condition=='macro_sh'|condition=='macro_ch'")
cpmlow_mac <- sm(normalize_expt(macrophage_human, filter="cbcb", convert="cpm"))
l2cpmlow_mac <- sm(normalize_expt(macrophage_human, filter="cbcb", convert="cpm", transform="log2"))
l2cpmpofa_mac <- sm(normalize_expt(macrophage_human, filter="pofa", convert="cpm", transform="log2"))
cpmrlelow_mac <- sm(normalize_expt(macrophage_human, filter="cbcb", convert="cpm", norm="rle"))
l2cpmrlelow_mac <- sm(normalize_expt(macrophage_human, filter="cbcb", convert="cpm", norm="rle", transform="log2"))

The normalizations in the block above do not attempt any batch correction. Now take a moment to view how they affect pca plots.

macrophage_human_pca <- plot_pca(macrophage_human)
## Start with the unmodified data
## It is clearly a mess.
macrophage_human_pca$plot

knitr::kable(macrophage_human_pca$table)

	sampleid	condition	batch	batch_int	PC1	PC2	colors	labels
hpgl0241	HPGL0241	macro_nil	a	1	-0.5683	-0.3528	#009900	HPGL0241
hpgl0242	HPGL0242	macro_sh	a	1	-0.1753	-0.2831	#000099	HPGL0242
hpgl0243	HPGL0243	macro_sh	a	1	-0.2902	-0.1149	#000099	HPGL0243
hpgl0244	HPGL0244	macro_ch	a	1	-0.0763	0.6999	#990000	HPGL0244
hpgl0245	HPGL0245	macro_ch	a	1	0.0355	0.0086	#990000	HPGL0245
hpgl0246	HPGL0246	macro_ch	a	1	-0.2197	0.3669	#990000	HPGL0246
hpgl0247	HPGL0247	macro_ch	a	1	-0.0294	-0.0494	#990000	HPGL0247
hpgl0248	HPGL0248	macro_ch	a	1	0.1368	0.0229	#990000	HPGL0248
hpgl0637	HPGL0637	macro_nil	bb	2	0.3007	0.1848	#009900	HPGL0637
hpgl0638	HPGL0638	macro_sh	bb	2	0.4894	-0.1976	#000099	HPGL0638
hpgl0639	HPGL0639	macro_sh	bb	2	0.3966	-0.2854	#000099	HPGL0639

cpmlow_mac_pca <- plot_pca(cpmlow_mac)
cpmlow_mac_pca$plot

## low-filtering and cpm is insufficient to clarify the data, but it is a start
knitr::kable(cpmlow_mac_pca$table)

	sampleid	condition	batch	batch_int	PC1	PC2	colors	labels
hpgl0241	HPGL0241	macro_nil	a	1	-0.4250	0.0094	#009900	HPGL0241
hpgl0242	HPGL0242	macro_sh	a	1	-0.4189	0.0515	#000099	HPGL0242
hpgl0243	HPGL0243	macro_sh	a	1	-0.2697	-0.0903	#000099	HPGL0243
hpgl0244	HPGL0244	macro_ch	a	1	0.3787	-0.4890	#990000	HPGL0244
hpgl0245	HPGL0245	macro_ch	a	1	-0.0227	-0.1623	#990000	HPGL0245
hpgl0246	HPGL0246	macro_ch	a	1	0.0817	-0.2484	#990000	HPGL0246
hpgl0247	HPGL0247	macro_ch	a	1	-0.1376	-0.1219	#990000	HPGL0247
hpgl0248	HPGL0248	macro_ch	a	1	0.0724	-0.2204	#990000	HPGL0248
hpgl0637	HPGL0637	macro_nil	bb	2	0.5165	0.2654	#009900	HPGL0637
hpgl0638	HPGL0638	macro_sh	bb	2	0.3406	0.3984	#000099	HPGL0638
hpgl0639	HPGL0639	macro_sh	bb	2	-0.1162	0.6077	#000099	HPGL0639

l2cpmlow_mac_pca <- plot_pca(l2cpmlow_mac)
l2cpmlow_mac_pca$plot

## Moving to the log scale makes much clearer the batch effect
## and shows that the data is actually pretty well separated by condition already.
knitr::kable(l2cpmlow_mac_pca$table)

	sampleid	condition	batch	batch_int	PC1	PC2	colors	labels
hpgl0241	HPGL0241	macro_nil	a	1	-0.2739	-0.6083	#009900	HPGL0241
hpgl0242	HPGL0242	macro_sh	a	1	-0.2363	-0.3176	#000099	HPGL0242
hpgl0243	HPGL0243	macro_sh	a	1	-0.1964	-0.1615	#000099	HPGL0243
hpgl0244	HPGL0244	macro_ch	a	1	-0.1604	0.4703	#990000	HPGL0244
hpgl0245	HPGL0245	macro_ch	a	1	-0.1480	0.2028	#990000	HPGL0245
hpgl0246	HPGL0246	macro_ch	a	1	-0.1549	0.2444	#990000	HPGL0246
hpgl0247	HPGL0247	macro_ch	a	1	-0.1324	0.1468	#990000	HPGL0247
hpgl0248	HPGL0248	macro_ch	a	1	-0.1605	0.2111	#990000	HPGL0248
hpgl0637	HPGL0637	macro_nil	bb	2	0.5074	0.1016	#009900	HPGL0637
hpgl0638	HPGL0638	macro_sh	bb	2	0.5063	0.0319	#000099	HPGL0638
hpgl0639	HPGL0639	macro_sh	bb	2	0.4492	-0.3216	#000099	HPGL0639

l2cpmpofa_mac_pca <- plot_pca(l2cpmpofa_mac)
l2cpmpofa_mac_pca$plot

## I have access to a few low-count filters, but they appear to make little or no difference.
knitr::kable(l2cpmpofa_mac_pca$table)

	sampleid	condition	batch	batch_int	PC1	PC2	colors	labels
hpgl0241	HPGL0241	macro_nil	a	1	-0.2784	-0.6197	#009900	HPGL0241
hpgl0242	HPGL0242	macro_sh	a	1	-0.2380	-0.3080	#000099	HPGL0242
hpgl0243	HPGL0243	macro_sh	a	1	-0.1973	-0.1527	#000099	HPGL0243
hpgl0244	HPGL0244	macro_ch	a	1	-0.1575	0.4586	#990000	HPGL0244
hpgl0245	HPGL0245	macro_ch	a	1	-0.1471	0.2057	#990000	HPGL0245
hpgl0246	HPGL0246	macro_ch	a	1	-0.1526	0.2483	#990000	HPGL0246
hpgl0247	HPGL0247	macro_ch	a	1	-0.1320	0.1522	#990000	HPGL0247
hpgl0248	HPGL0248	macro_ch	a	1	-0.1590	0.2123	#990000	HPGL0248
hpgl0637	HPGL0637	macro_nil	bb	2	0.5043	0.0919	#009900	HPGL0637
hpgl0638	HPGL0638	macro_sh	bb	2	0.5071	0.0361	#000099	HPGL0638
hpgl0639	HPGL0639	macro_sh	bb	2	0.4506	-0.3249	#000099	HPGL0639

cpmrlelow_mac_pca <- plot_pca(cpmrlelow_mac)
cpmrlelow_mac_pca$plot

## RLE normalization has an interesting effect on the data
knitr::kable(cpmrlelow_mac_pca$table)

	sampleid	condition	batch	batch_int	PC1	PC2	colors	labels
hpgl0241	HPGL0241	macro_nil	a	1	-0.3438	0.2672	#009900	HPGL0241
hpgl0242	HPGL0242	macro_sh	a	1	-0.4137	-0.0516	#000099	HPGL0242
hpgl0243	HPGL0243	macro_sh	a	1	-0.3644	-0.1527	#000099	HPGL0243
hpgl0244	HPGL0244	macro_ch	a	1	0.4064	0.5885	#990000	HPGL0244
hpgl0245	HPGL0245	macro_ch	a	1	-0.0250	0.1166	#990000	HPGL0245
hpgl0246	HPGL0246	macro_ch	a	1	0.0203	-0.0543	#990000	HPGL0246
hpgl0247	HPGL0247	macro_ch	a	1	-0.1125	0.2384	#990000	HPGL0247
hpgl0248	HPGL0248	macro_ch	a	1	0.0742	0.1733	#990000	HPGL0248
hpgl0637	HPGL0637	macro_nil	bb	2	0.5164	-0.2837	#009900	HPGL0637
hpgl0638	HPGL0638	macro_sh	bb	2	0.3419	-0.3257	#000099	HPGL0638
hpgl0639	HPGL0639	macro_sh	bb	2	-0.0999	-0.5161	#000099	HPGL0639

l2cpmrlelow_mac_pca <- plot_pca(l2cpmrlelow_mac, size_column="pctcategory")
l2cpmrlelow_mac_pca$plot

## Still the batch effect is pretty ridiculous
## This plot resizes the glyphs to indicate the % mapping.
knitr::kable(l2cpmrlelow_mac_pca$table)

	sampleid	condition	batch	batch_int	PC1	PC2	colors	labels	pctcategory
hpgl0241	HPGL0241	macro_nil	a	1	-0.2774	0.6038	#009900	HPGL0241	2
hpgl0242	HPGL0242	macro_sh	a	1	-0.2373	0.3217	#000099	HPGL0242	5
hpgl0243	HPGL0243	macro_sh	a	1	-0.1928	0.1581	#000099	HPGL0243	5
hpgl0244	HPGL0244	macro_ch	a	1	-0.1631	-0.4755	#990000	HPGL0244	3
hpgl0245	HPGL0245	macro_ch	a	1	-0.1475	-0.2064	#990000	HPGL0245	6
hpgl0246	HPGL0246	macro_ch	a	1	-0.1462	-0.2293	#990000	HPGL0246	4
hpgl0247	HPGL0247	macro_ch	a	1	-0.1350	-0.1635	#990000	HPGL0247	4
hpgl0248	HPGL0248	macro_ch	a	1	-0.1627	-0.2027	#990000	HPGL0248	6
hpgl0637	HPGL0637	macro_nil	bb	2	0.5095	-0.0934	#009900	HPGL0637	2
hpgl0638	HPGL0638	macro_sh	bb	2	0.5072	-0.0407	#000099	HPGL0638	7
hpgl0639	HPGL0639	macro_sh	bb	2	0.4453	0.3279	#000099	HPGL0639	5

1.2.1 Look at the principal components

macro_pca_info <- sm(pca_information(macrophage_human,
                                     expt_factors=c("condition","batch","pctmappedparasite",
                                                    "state","pathogenstrain","rnangul")))

##write.csv(macro_pca_info$rsquared_table, file="images/rsquared.csv")
##write.csv(macro_pca_info$anova_p, file="images/anova_p.csv")

macro_pca_info$anova_fstat_heatmap

I could continue the above normalizations and plots, but it quickly becomes boring and redundant. The batch effect is too significant to continue without addressing it. One thing that was brought up in previous discussions: does the mapping percentage act as a surrogate variable? This suggests that if it does, it is pretty weak, and the last PCA plot shows that is cosegregates with the existing conditions: uninfected has the lowest mapping percent to the parasite [0%]; while chronic has the lower percentage of the infections, but not so much lower than self-healing in general.

1.2.2 Normalization with batch effect minimization

We have good information that suggests a strong batch effect in this data. Let us put that information to use in the following block.

How do these modifications affect the separation of the data?

cpmlowcbcb_mac <- sm(normalize_expt(macrophage_human, filter="cbcb", convert="cpm", batch=TRUE))
l2cpmlowcom_mac <- sm(normalize_expt(macrophage_human, filter="cbcb", convert="cpm",
                                     batch="combat_scale", transform="log2"))

l2cpmrlelowcom_mac <- sm(normalize_expt(macrophage_human, filter="cbcb", convert="cpm",
                                        batch="combat_scale", transform="log2", norm="rle"))

l2cpmquantlowcmod_mac <- sm(normalize_expt(macrophage_human, filter="cbcb", convert="cpm",
                                           batch="combatmod", transform="log2", norm="quant"))
l2cpmquantlowcom_mac <- sm(normalize_expt(macrophage_human, filter="cbcb", convert="cpm",
                                          batch="combat_scale", transform="log2", norm="quant"))

l2cpmquantlowsva_mac <- sm(normalize_expt(macrophage_human, filter="cbcb", convert="cpm",
                                          batch="svaseq", norm="quant"))
l2cpmquantlowfsva_mac <- sm(normalize_expt(macrophage_human, filter="cbcb", convert="cpm",
                                          batch="fsva", norm="quant"))
l2cpmquantlowruv_mac <- sm(normalize_expt(macrophage_human, filter="cbcb", convert="cpm",
                                          batch="combat_scale", transform="log2", norm="quant"))

tt <- sm(graph_metrics(l2cpmquantlowruv_mac))

cpmlowcbcb_mac_pca <- plot_pca(cpmlowcbcb_mac)
cpmlowcbcb_mac_pca$plot

## cbcb's batch correction is insufficient, at least without log2
cpmlowcbcb_mac_pca$table

##          sampleid condition batch batch_int      PC1      PC2  colors   labels
## hpgl0241 HPGL0241 macro_nil     a         1 -0.09337  0.66888 #009900 HPGL0241
## hpgl0242 HPGL0242  macro_sh     a         1 -0.23685 -0.15485 #000099 HPGL0242
## hpgl0243 HPGL0243  macro_sh     a         1 -0.06401 -0.01482 #000099 HPGL0243
## hpgl0244 HPGL0244  macro_ch     a         1  0.60514 -0.04000 #990000 HPGL0244
## hpgl0245 HPGL0245  macro_ch     a         1  0.09059 -0.20559 #990000 HPGL0245
## hpgl0246 HPGL0246  macro_ch     a         1  0.21245 -0.28906 #990000 HPGL0246
## hpgl0247 HPGL0247  macro_ch     a         1 -0.05050 -0.20619 #990000 HPGL0247
## hpgl0248 HPGL0248  macro_ch     a         1  0.19782 -0.22250 #990000 HPGL0248
## hpgl0637 HPGL0637 macro_nil    bb         2  0.13221  0.51916 #009900 HPGL0637
## hpgl0638 HPGL0638  macro_sh    bb         2 -0.13449  0.11016 #000099 HPGL0638
## hpgl0639 HPGL0639  macro_sh    bb         2 -0.65901 -0.16518 #000099 HPGL0639

l2cpmlowcom_mac_pca <- plot_pca(l2cpmlowcom_mac)
l2cpmlowcom_mac_pca$plot

## oooo pretty
knitr::kable(l2cpmlowcom_mac_pca$table)

	sampleid	condition	batch	batch_int	PC1	PC2	colors	labels
hpgl0241	HPGL0241	macro_nil	a	1	-0.5127	0.0450	#009900	HPGL0241
hpgl0242	HPGL0242	macro_sh	a	1	-0.2641	-0.1765	#000099	HPGL0242
hpgl0243	HPGL0243	macro_sh	a	1	-0.1075	0.2984	#000099	HPGL0243
hpgl0244	HPGL0244	macro_ch	a	1	0.4668	-0.1982	#990000	HPGL0244
hpgl0245	HPGL0245	macro_ch	a	1	0.2520	0.2666	#990000	HPGL0245
hpgl0246	HPGL0246	macro_ch	a	1	0.2854	-0.4778	#990000	HPGL0246
hpgl0247	HPGL0247	macro_ch	a	1	0.2180	0.4763	#990000	HPGL0247
hpgl0248	HPGL0248	macro_ch	a	1	0.2529	-0.2399	#990000	HPGL0248
hpgl0637	HPGL0637	macro_nil	bb	2	-0.0447	0.2786	#009900	HPGL0637
hpgl0638	HPGL0638	macro_sh	bb	2	-0.1499	0.1284	#000099	HPGL0638
hpgl0639	HPGL0639	macro_sh	bb	2	-0.3962	-0.4009	#000099	HPGL0639

l2cpmrlelowcom_mac_pca <- plot_pca(l2cpmrlelowcom_mac)
l2cpmrlelowcom_mac_pca$plot

## also looking good
knitr::kable(l2cpmrlelowcom_mac_pca$table)

	sampleid	condition	batch	batch_int	PC1	PC2	colors	labels
hpgl0241	HPGL0241	macro_nil	a	1	-0.4928	0.2216	#009900	HPGL0241
hpgl0242	HPGL0242	macro_sh	a	1	-0.2838	-0.1072	#000099	HPGL0242
hpgl0243	HPGL0243	macro_sh	a	1	-0.1208	0.1463	#000099	HPGL0243
hpgl0244	HPGL0244	macro_ch	a	1	0.4726	-0.0200	#990000	HPGL0244
hpgl0245	HPGL0245	macro_ch	a	1	0.2551	0.0956	#990000	HPGL0245
hpgl0246	HPGL0246	macro_ch	a	1	0.2789	-0.3640	#990000	HPGL0246
hpgl0247	HPGL0247	macro_ch	a	1	0.2226	0.4386	#990000	HPGL0247
hpgl0248	HPGL0248	macro_ch	a	1	0.2500	-0.3830	#990000	HPGL0248
hpgl0637	HPGL0637	macro_nil	bb	2	-0.0367	0.3580	#009900	HPGL0637
hpgl0638	HPGL0638	macro_sh	bb	2	-0.1424	0.1492	#000099	HPGL0638
hpgl0639	HPGL0639	macro_sh	bb	2	-0.4027	-0.5352	#000099	HPGL0639

l2cpmquantlowcmod_mac_pca <- plot_pca(l2cpmquantlowcmod_mac)
l2cpmquantlowcmod_mac_pca$plot

## not as good
knitr::kable(l2cpmquantlowcmod_mac_pca$table)

	sampleid	condition	batch	batch_int	PC1	PC2	colors	labels
hpgl0241	HPGL0241	macro_nil	a	1	-0.6808	0.1257	#009900	HPGL0241
hpgl0242	HPGL0242	macro_sh	a	1	-0.3746	-0.1335	#000099	HPGL0242
hpgl0243	HPGL0243	macro_sh	a	1	-0.2021	0.2806	#000099	HPGL0243
hpgl0244	HPGL0244	macro_ch	a	1	0.3934	-0.1607	#990000	HPGL0244
hpgl0245	HPGL0245	macro_ch	a	1	0.1805	0.2191	#990000	HPGL0245
hpgl0246	HPGL0246	macro_ch	a	1	0.1726	-0.4126	#990000	HPGL0246
hpgl0247	HPGL0247	macro_ch	a	1	0.1340	0.3985	#990000	HPGL0247
hpgl0248	HPGL0248	macro_ch	a	1	0.1736	-0.2283	#990000	HPGL0248
hpgl0637	HPGL0637	macro_nil	bb	2	0.2009	0.3443	#009900	HPGL0637
hpgl0638	HPGL0638	macro_sh	bb	2	0.1591	0.1152	#000099	HPGL0638
hpgl0639	HPGL0639	macro_sh	bb	2	-0.1566	-0.5484	#000099	HPGL0639

l2cpmquantlowfsva_mac_pca <- plot_pca(l2cpmquantlowfsva_mac)
l2cpmquantlowfsva_mac_pca$plot

pretty_pca <- plot_pca(l2cpmquantlowcom_mac)
pretty_pca$plot

## perhaps the best?
knitr::kable(pretty_pca$table)

	sampleid	condition	batch	batch_int	PC1	PC2	colors	labels
hpgl0241	HPGL0241	macro_nil	a	1	-0.5168	0.0403	#009900	HPGL0241
hpgl0242	HPGL0242	macro_sh	a	1	-0.2662	-0.1860	#000099	HPGL0242
hpgl0243	HPGL0243	macro_sh	a	1	-0.1145	0.3031	#000099	HPGL0243
hpgl0244	HPGL0244	macro_ch	a	1	0.4741	-0.1917	#990000	HPGL0244
hpgl0245	HPGL0245	macro_ch	a	1	0.2622	0.2748	#990000	HPGL0245
hpgl0246	HPGL0246	macro_ch	a	1	0.2692	-0.4869	#990000	HPGL0246
hpgl0247	HPGL0247	macro_ch	a	1	0.2171	0.4809	#990000	HPGL0247
hpgl0248	HPGL0248	macro_ch	a	1	0.2625	-0.2296	#990000	HPGL0248
hpgl0637	HPGL0637	macro_nil	bb	2	-0.0613	0.2625	#009900	HPGL0637
hpgl0638	HPGL0638	macro_sh	bb	2	-0.1519	0.1251	#000099	HPGL0638
hpgl0639	HPGL0639	macro_sh	bb	2	-0.3742	-0.3925	#000099	HPGL0639

##write.csv(pretty_pca$table, file="excel/pretty_pca.csv")

1.2.3 Surrogate estimations

## Now try and combine this with
## https://github.com/jtleek/svaseq/blob/master/recount.Rmd
## To see if things become clearer
human_filt <- sm(normalize_expt(expt=macrophage_human, filter=TRUE))
surrogate_estimations <- compare_surrogate_estimates(human_filt)

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

## Attempting pca surrogate estimation.

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

## Attempting sva supervised surrogate estimation.

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

## Attempting sva unsupervised surrogate estimation.

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

## Attempting ruvseq supervised surrogate estimation.

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

## Attempting ruvseq residual surrogate estimation.

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

## Attempting ruvseq empirical surrogate estimation.

## 1/8: Performing lmFit(data) etc. with null in the model.

## 2/8: Performing lmFit(data) etc. with + batch_adjustments$batch in the model.

## 3/8: Performing lmFit(data) etc. with + batch_adjustments$pca in the model.

## 4/8: Performing lmFit(data) etc. with + batch_adjustments$sva_sup in the model.

## 5/8: Performing lmFit(data) etc. with + batch_adjustments$sva_unsup in the model.

## 6/8: Performing lmFit(data) etc. with + batch_adjustments$ruv_sup in the model.

## 7/8: Performing lmFit(data) etc. with + batch_adjustments$ruv_resid in the model.

## 8/8: Performing lmFit(data) etc. with + batch_adjustments$ruv_emp in the model.

## I think my implementation of ruv_empirical may be borken.
surrogate_estimations$pca_plots$ruvsup

surrogate_estimations$ruv_supervised_adjust$svs_sample

## NULL

surrogate_estimations$pca_plots$svasup

surrogate_estimations$sva_supervised_adjust$svs_sample

## NULL

## The new uninfected sample still seems oddly placed
surrogate_estimations$pca_plots$pca

surrogate_estimations$pca_adjust$svs_sample

## NULL

## Ditto here
surrogate_estimations$pca_plots$svaunsup

## This seems pretty reasonable to me
surrogate_estimations$pca_plots$ruvresid

## This looks a lot like pca and ruvsup, but the corrplot suggested to me that it would look weirder.

1.2.4 Effect on other metrics by batch correction

Lets see how these modification affect the visualizations

batchnorm_metrics <- sm(graph_metrics(l2cpmquantlowcom_mac))

print(batchnorm_metrics$corheat)

## The combat adjusted counts seem to split  relatively clearly between types.
print(batchnorm_metrics$smc)

## This one sample is still problematic, though.
print(batchnorm_metrics$pcaplot)

## But quantile normalization really does smash them into the same distribution.

Some other analyses performed suggest a possible switch between two samples. We can artifically see what the data would look like if that is true by switching those two samples in a separate experimental design.

2 TODO 201611 #4

PBMC: variance partition human, “Are there any genes with variance we can ascribe to condition?” That set of genes is the most interesting. – use the variance parition table and pull the top 200. Report back the % variance by condition for all of these genes.

index.html

3 Variance Partition

In the following, I will attempt to find the variance associated with different experimental factors in the macrophage data with mappings against the human transcriptome.

3.1 Setting up

As in macrophage_estimation_human.html, I need to pull out the relevant samples:

human_expt$notes

## [1] "New experimental design factors by snp added 2016-09-20"

colnames(human_expt$design)

##  [1] "sampleid"                        "experimentname"                 
##  [3] "tubelabel"                       "alias"                          
##  [5] "condition"                       "batch"                          
##  [7] "anotherbatch"                    "snpclade"                       
##  [9] "snpcladev2"                      "snpcladev3"                     
## [11] "pathogenstrain"                  "donor"                          
## [13] "time"                            "pctmappedparasite"              
## [15] "pctcategory"                     "state"                          
## [17] "sourcelab"                       "expperson"                      
## [19] "pathogen"                        "host"                           
## [21] "hostcelltype"                    "noofhostcells"                  
## [23] "infectionperiodhpitimeofharvest" "moiexposure"                    
## [25] "parasitespercell"                "pctinf"                         
## [27] "rnangul"                         "rnaqcpassed"                    
## [29] "libraryconst"                    "libqcpassed"                    
## [31] "index"                           "descriptonandremarks"           
## [33] "observation"                     "lowercaseid"                    
## [35] "humanfile"                       "parasitefile"                   
## [37] "file"

## make sure that I am getting the material from ~ 2016-09-20
macrophage_human <- expt_subset(human_expt, subset="experimentname=='macrophage'")
macrophage_human <- set_expt_colors(macrophage_human, colors=c("#005500","#880000","#0000FF"))
macronorm_human <- sm(normalize_expt(macrophage_human, filter=TRUE))

4 Perform variancePartition analyses

4.1 Condition + batch

Start out with just condition and batch

## The default arguments are to query ~ condition + (1|batch)
## Which, because condition is categorical, will end badly.
vp_cb <- sm(varpart(macronorm_human,
                    predictor=NULL,
                    factors=c("condition","batch")))
vp_cb$percent_plot

vp_cb$partition_plot

## Batch is dominant, but there is still quite a lot left in residuals.

4.2 Condition + batch + snpclades

Try using the snp information and see if it helps. Unfortunately, I have not yet processed one of the strains, that should be addressed asap.

The following all fail with “Downdated VtV is not positive definite” wtf ever that means. My google searches have so far proven useless for this.

vp_cbs3 <- sm(varpart(macronorm_human,
                      predictor=NULL,
                      factors=c("condition","batch","snpcladev3")))
vp_cbs3$percent_plot

vp_cbs3$partition_plot

nob <- expt_subset(expt=macronorm_human, subset="batch!='b'")
vp_nob <- sm(varpart(expt=nob, predictor=NULL, factors=c("condition", "snpcladev3")))
knitr::kable(varpart_human_nob$sorted_df, n=100)

## Error in inherits(x, "list"): object 'varpart_human_nob' not found

png(file="images/withb_violin.png")
vp_cbs3$partition_plot
dev.off()

## X11cairo 
##        2

png(file="images/withb_percent.png")
vp_cbs3$percent_plot
dev.off()

## X11cairo 
##        2

png(file="images/nob_violin.png")
vp_nob$partition_plot
dev.off()

## X11cairo 
##        2

png(file="images/nob_percent.png")
vp_nob$percent_plot
dev.off()

## X11cairo 
##        2

index.html annotation.html infection_estimation.html

this_save <- paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
message(paste0("Saving to: ", this_save))

## Saving to: macrophage_estimation_human-v20170612.rda.xz

tmp <- sm(saveme(filename=this_save))
library("pander")
pander(sessionInfo())

R version 3.3.3 (2017-03-06)

**Platform:** x86_64-pc-linux-gnu (64-bit)

locale: LC_CTYPE=en_US.utf8, LC_NUMERIC=C, LC_TIME=en_US.utf8, LC_COLLATE=en_US.utf8, LC_MONETARY=en_US.utf8, LC_MESSAGES=en_US.utf8, LC_PAPER=en_US.utf8, LC_NAME=C, LC_ADDRESS=C, LC_TELEPHONE=C, LC_MEASUREMENT=en_US.utf8 and LC_IDENTIFICATION=C

attached base packages: parallel, stats4, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: hpgltools(v.2017.01), ruv(v.0.9.6), ggplot2(v.2.2.1), directlabels(v.2017.03.31), AnnotationHub(v.2.6.5), Homo.sapiens(v.1.3.1), TxDb.Hsapiens.UCSC.hg19.knownGene(v.3.2.2), org.Hs.eg.db(v.3.4.0), Leishmania.braziliensis.MHOMBR75M2904(v.28.0), TxDb.lbraziliensis.MHOMBR75M2904.TriTrypDB28(v.28.0), org.Lbraziliensis.MHOMBR75M2904.eg.db(v.28.0), Leishmania.panamensis.MHOMCOL81L13(v.27.0), GO.db(v.3.4.0), OrganismDbi(v.1.16.0), TxDb.lpanamensis.MHOMCOL81L13.TriTrypDB27(v.27.0), GenomicFeatures(v.1.26.4), GenomicRanges(v.1.26.4), GenomeInfoDb(v.1.10.3), org.Lpanamensis.MHOMCOL81L13.eg.db(v.27.0), AnnotationDbi(v.1.36.2), IRanges(v.2.8.2), S4Vectors(v.0.12.2), Biobase(v.2.34.0), BiocGenerics(v.0.20.0) and pander(v.0.6.0)

loaded via a namespace (and not attached): backports(v.1.1.0), corrplot(v.0.77), aroma.light(v.3.4.0), plyr(v.1.8.4), lazyeval(v.0.2.0), survJamda.data(v.1.0.2), splines(v.3.3.3), BiocParallel(v.1.8.2), sva(v.3.22.0), digest(v.0.6.12), SuppDists(v.1.1-9.4), foreach(v.1.4.3), BiocInstaller(v.1.24.0), htmltools(v.0.3.6), gdata(v.2.18.0), magrittr(v.1.5), memoise(v.1.1.0), survJamda(v.1.1.4), doParallel(v.1.0.10), openxlsx(v.4.0.17), limma(v.3.30.13), Biostrings(v.2.42.1), annotate(v.1.52.1), matrixStats(v.0.52.2), R.utils(v.2.5.0), colorspace(v.1.3-2), blob(v.1.1.0), ggrepel(v.0.6.5), crayon(v.1.3.2), RCurl(v.1.95-4.8), graph(v.1.52.0), roxygen2(v.6.0.1), genefilter(v.1.56.0), lme4(v.1.1-13), survival(v.2.41-3), iterators(v.1.0.8), gtable(v.0.2.0), zlibbioc(v.1.20.0), XVector(v.0.14.1), DESeq(v.1.26.0), scales(v.0.4.1), DBI(v.0.7), edgeR(v.3.16.5), Rcpp(v.0.12.11), xtable(v.1.8-2), bit(v.1.1-12), preprocessCore(v.1.36.0), lava(v.1.5), prodlim(v.1.6.1), ecodist(v.1.2.9), httr(v.1.2.1), gplots(v.3.0.1), RColorBrewer(v.1.1-2), R.methodsS3(v.1.7.1), pkgconfig(v.2.0.1), XML(v.3.98-1.9), locfit(v.1.5-9.1), labeling(v.0.3), rlang(v.0.1.1), reshape2(v.1.4.2), munsell(v.0.4.3), tools(v.3.3.3), RSQLite(v.2.0), devtools(v.1.13.2), evaluate(v.0.10), stringr(v.1.2.0), yaml(v.2.1.14), bootstrap(v.2017.2), knitr(v.1.16), bit64(v.0.9-7), caTools(v.1.17.1), EDASeq(v.2.8.0), RBGL(v.1.50.0), nlme(v.3.1-131), mime(v.0.5), R.oo(v.1.21.0), xml2(v.1.1.1), biomaRt(v.2.30.0), compiler(v.3.3.3), pbkrtest(v.0.4-7), curl(v.2.6), interactiveDisplayBase(v.1.12.0), variancePartition(v.1.4.2), testthat(v.1.0.2), statmod(v.1.4.30), geneplotter(v.1.52.0), tibble(v.1.3.3), stringi(v.1.1.5), highr(v.0.6), lattice(v.0.20-35), Matrix(v.1.2-10), commonmark(v.1.2), survcomp(v.1.24.0), nloptr(v.1.0.4), survivalROC(v.1.0.3), data.table(v.1.10.4), bitops(v.1.0-6), corpcor(v.1.6.9), httpuv(v.1.3.3), rtracklayer(v.1.34.2), colorRamps(v.2.3), latticeExtra(v.0.6-28), hwriter(v.1.3.2), R6(v.2.2.2), ShortRead(v.1.32.1), KernSmooth(v.2.23-15), codetools(v.0.2-15), MASS(v.7.3-47), gtools(v.3.5.0), SummarizedExperiment(v.1.4.0), rprojroot(v.1.2), RUVSeq(v.1.8.0), withr(v.1.0.2), GenomicAlignments(v.1.10.1), Rsamtools(v.1.26.2), mgcv(v.1.8-17), quadprog(v.1.5-5), grid(v.3.3.3), minqa(v.1.2.4), rmarkdown(v.1.6), shiny(v.1.0.3), base64enc(v.0.1-3) and rmeta(v.2.16)