RNAseq of L.panamensis: PBMC Infection Differential Expression.

index.html 01_annotation.html 02_estimation_infection.html

1 PBMC Infection Differential Expression, Infection: 20170820

This document turns to the infection of PBMC cells with L.panamensis. This data is particularly strangely affected by the different strains used to infect the cells, and as a result is both useful and troubling.

Given the observations above, we have some ideas of ways to pass the data for differential expression analyses which may or may not be ‘better’. Lets try some and see what happens.

1.1 Create data sets to compare differential expression analyses

Given the above ways to massage the data, lets use a few of them for limma/deseq/edger. The main caveat in this is that those tools really do expect specific distributions of data which we horribly violate if we use log2() data, which is why in the previous blocks I named them l2blahblah, thus we can do the same sets of normalization but without that and forcibly push the resulting data into limma/edger/deseq.

2 The negative control

Everything I did in 02_estimation_infection.html suggests that there are no significant differences visible if one looks just at chronic/self-healing in this data. Further testing has seemingly proven this statement, as a result most of the following analyses will look at chronic/uninfected and self-healing/uninfected followed by attempts to reconcile those results.

2.1 Filter the data

To save some time and annoyance with sva, lets filter the data now. In addition, write down a small function used to extract the sets of significant genes across different contrasts (notably self/uninfected vs. chronic/uninfected).

hs_inf_filt <- sm(normalize_expt(hs_inf, filter=TRUE))
hs_uninf_filt <- sm(normalize_expt(hs_uninf, filter=TRUE))
keepers <- list("sh_nil" = c("sh", "uninf"),
                "ch_nil" = c("chr", "uninf"),
                "ch_sh" = c("chr", "sh"))

subset_significants <- function(hs_sig) {
    sh_nil_up_genes <- rownames(hs_sig[["deseq"]][["ups"]][["sh_vs_uninf"]])
    ch_nil_up_genes <- rownames(hs_sig[["deseq"]][["ups"]][["chr_vs_uninf"]])
    sh_nil_down_genes <- rownames(hs_sig[["deseq"]][["downs"]][["sh_vs_uninf"]])
    ch_nil_down_genes <- rownames(hs_sig[["deseq"]][["downs"]][["chr_vs_uninf"]])

    sh_solo_up_idx <- ! sh_nil_up_genes %in% ch_nil_up_genes
    sh_solo_up <- sh_nil_up_genes[sh_solo_up_idx]
    sh_shared_ch_idx <- sh_nil_up_genes %in% ch_nil_up_genes
    sh_shared_ch_up <- sh_nil_up_genes[sh_shared_ch_idx]
    ch_solo_up_idx <- ! ch_nil_up_genes %in% sh_nil_up_genes
    ch_solo_up <- ch_nil_up_genes[ch_solo_up_idx]

    sh_solo_down_idx <- ! sh_nil_down_genes %in% ch_nil_down_genes
    sh_solo_down <- sh_nil_down_genes[sh_solo_down_idx]
    sh_shared_ch_idx <- sh_nil_down_genes %in% ch_nil_down_genes
    sh_shared_ch_down <- sh_nil_down_genes[sh_shared_ch_idx]
    ch_solo_down_idx <- ! ch_nil_down_genes %in% sh_nil_down_genes
    ch_solo_down <- ch_nil_down_genes[ch_solo_down_idx]

    retlist <- list(
        "sh_solo_up" = hs_sig[["deseq"]][["ups"]][["sh_vs_uninf"]][sh_solo_up, ],
        "ch_solo_up" = hs_sig[["deseq"]][["ups"]][["chr_vs_uninf"]][ch_solo_up, ],
        "sh_shared_ch_up" = hs_sig[["deseq"]][["ups"]][["sh_vs_uninf"]][sh_shared_ch_up, ],
        "sh_solo_down" = hs_sig[["deseq"]][["downs"]][["sh_vs_uninf"]][sh_solo_down, ],
        "ch_solo_down" = hs_sig[["deseq"]][["downs"]][["chr_vs_uninf"]][ch_solo_down, ],
        "sh_shared_ch_down" = hs_sig[["deseq"]][["downs"]][["sh_vs_uninf"]][sh_shared_ch_down, ])
    retlist[["up_weights"]] <- c(0, nrow(retlist[["sh_solo_up"]]),
                                 nrow(retlist[["ch_solo_up"]]), nrow(retlist[["sh_shared_ch_up"]]))
    retlist[["down_weights"]] <- c(0, nrow(retlist[["sh_solo_down"]]),
                                 nrow(retlist[["ch_solo_down"]]), nrow(retlist[["sh_shared_ch_down"]]))
    retlist[["up_venn"]] <- Vennerable::Venn(SetNames = c("sh", "chr"),
                                             Weight = retlist[["up_weights"]])
    retlist[["down_venn"]] <- Vennerable::Venn(SetNames = c("sh", "chr"),
                                               Weight = retlist[["down_weights"]])
    return(retlist)
}

2.2 Do a completely normal limma invocation.

The following probably should not be used.

counts <- exprs(hs_uninf_filt)
design <- pData(hs_uninf_filt)
model <- model.matrix(~ 0 + condition + donor + pathogenstrain, data=design)
voom_weight_result <- limma::voomWithQualityWeights(counts=counts, design=model,
                                                    normalize.method="quantile", plot=TRUE)
voom_result <- limma::voom(counts=counts, design=model, normalize.method="quantile", plot=TRUE)
fitting_weight <- limma::lmFit(object=voom_weight_result, design=model, method="ls")
fitting <- limma::lmFit(object=voom_result, design=model, method="ls")
contrast <- limma::makeContrasts(sh_ch=conditionpbmc_sh-conditionpbmc_ch, levels=model)
contrast_weight <- limma::contrasts.fit(fit=fitting_weight, contrasts=contrast)
contrast <- limma::contrasts.fit(fit=fitting, contrasts=contrast)
ebayes_weighted <- limma::eBayes(contrast_weight)
ebayes <- limma::eBayes(contrast)
toptable_weighted <- limma::topTable(ebayes_weighted)
toptable <- limma::topTable(ebayes)
toptable

limma_test <- limma_pairwise(input=hs_uninf_filt,
                             alt_model="~ 0 + condition + donor + donor:pathogenstrain")
limma_top <- limma::topTable(limma_test$pairwise_comparisons,
                             number=nrow(limma_test$pairwise_comparisons),
                             sort.by="P",
                             coef="pbmc_ch_vs_pbmc_sh")
head(limma_top, n=10)

limma_test2 <- limma_pairwise(input=hs_uninf_filt,
                              alt_model="~ 0 + condition + donor")
limma_top2 <- limma::topTable(limma_test2$pairwise_comparisons,
                              coef="pbmc_ch_vs_pbmc_sh",
                              number=nrow(limma_test2$pairwise_comparisons),
                              sort.by="P")
head(limma_top2)

hs_pairwise_nobatch <- sm(all_pairwise(hs_uninf_filt, model_batch=FALSE))

hs_combined_nobatch <- sm(combine_de_tables(hs_pairwise_nobatch,
                                            excel=paste0("excel/hs_infect_nobatch-v", ver, ".xlsx"),
                                            keepers=keepers))

hs_sig_nobatch <- sm(extract_significant_genes(hs_combined_nobatch,
                                               excel=paste0("excel/hs_infect_nobatch_sig-v", ver, ".xlsx")))
hs_sig_nobatch$deseq$counts

##              change_counts_up change_counts_down
## sh_vs_uninf              1050                193
## chr_vs_uninf             1052                171
## chr_vs_sh                   0                  0

common_solos_nobatch <- subset_significants(hs_sig_nobatch)
summary(common_solos_nobatch)

##                   Length Class      Mode   
## sh_solo_up        38     data.frame list   
## ch_solo_up        38     data.frame list   
## sh_shared_ch_up   38     data.frame list   
## sh_solo_down      38     data.frame list   
## ch_solo_down      38     data.frame list   
## sh_shared_ch_down 38     data.frame list   
## up_weights         4     -none-     numeric
## down_weights       4     -none-     numeric
## up_venn            1     Venn       S4     
## down_venn          1     Venn       S4

nobatch_up_venn <- Vennerable::plot(common_solos_nobatch$up_venn, doWeights=FALSE)

nobatch_down_venn <- Vennerable::plot(common_solos_nobatch$down_venn, doWeights=FALSE)

3 Add patient to the model

Repeat the previous set of analyses with d107/108/110 in the model.

hs_pairwise_batch <- sm(all_pairwise(hs_uninf_filt, model_batch=TRUE))

hs_combined_batch <- sm(combine_de_tables(hs_pairwise_batch,
                                          excel=paste0("excel/hs_infect_patbatch-v", ver, ".xlsx"),
                                          keepers=keepers))

hs_sig_batch <- sm(extract_significant_genes(hs_combined_batch,
                                             excel=paste0("excel/hs_infect_patbatch_sig-v", ver, ".xlsx")))
hs_sig_batch[["deseq"]][["counts"]]

##              change_counts_up change_counts_down
## sh_vs_uninf              1035                330
## chr_vs_uninf             1060                294
## chr_vs_sh                   0                  0

common_solos_batch <- subset_significants(hs_sig_batch)
summary(common_solos_batch)

##                   Length Class      Mode   
## sh_solo_up        38     data.frame list   
## ch_solo_up        38     data.frame list   
## sh_shared_ch_up   38     data.frame list   
## sh_solo_down      38     data.frame list   
## ch_solo_down      38     data.frame list   
## sh_shared_ch_down 38     data.frame list   
## up_weights         4     -none-     numeric
## down_weights       4     -none-     numeric
## up_venn            1     Venn       S4     
## down_venn          1     Venn       S4

similar <- sm(compare_de_results(hs_combined_nobatch, hs_combined_batch, cor_method="spearman"))
similar[["deseq"]]

## $sh_vs_uninf
## $sh_vs_uninf$logfc
## [1] 0.992
## 
## $sh_vs_uninf$p
## [1] 0.9339
## 
## $sh_vs_uninf$adjp
## [1] 0.9339
## 
## 
## $chr_vs_uninf
## $chr_vs_uninf$logfc
## [1] 0.9926
## 
## $chr_vs_uninf$p
## [1] 0.9361
## 
## $chr_vs_uninf$adjp
## [1] 0.9361
## 
## 
## $chr_vs_sh
## $chr_vs_sh$logfc
## [1] 0.9854
## 
## $chr_vs_sh$p
## [1] 0.9266
## 
## $chr_vs_sh$adjp
## [1] 0.0953

batch_up_venn <- Vennerable::plot(common_solos_batch[["up_venn"]], doWeights=FALSE)

batch_down_venn <- Vennerable::plot(common_solos_batch[["down_venn"]], doWeights=FALSE)

kept_columns <- c("transcriptid", "geneid", "description", "deseq_logfc", "deseq_adjp")
xls_result <- write_xls(data=common_solos_batch[["sh_solo_up"]][, kept_columns], sheet="sh_up_solo")
xls_result <- write_xls(data=common_solos_batch[["ch_solo_up"]][, kept_columns], sheet="chr_up_solo",
                        wb=xls_result[["workbook"]])
xls_result <- write_xls(data=common_solos_batch[["sh_shared_ch_up"]][, kept_columns],
                        sheet="shchr_shared_up", wb=xls_result[["workbook"]])
xls_result <- write_xls(data=common_solos_batch[["sh_solo_down"]][, kept_columns],
                        sheet="sh_down_solo", wb=xls_result[["workbook"]])
xls_result <- write_xls(data=common_solos_batch[["ch_solo_down"]][, kept_columns],
                        sheet="chr_down_solo", wb=xls_result[["workbook"]])
xls_result <- write_xls(data=common_solos_batch[["sh_shared_ch_down"]][, kept_columns],
                        sheet="shchr_shared_down", wb=xls_result[["workbook"]],
                        excel=paste0("excel/figure_5a_stuff-v", ver, ".xlsx"))

## Saving to: excel/figure_5a_stuff-v20170820.xlsx

4 Add sva into the mix

Repeat, this time attmepting to tamp down the variance by person.

hs_pairwise_ssva <- sm(all_pairwise(hs_uninf_filt, model_batch="ssva"))

hs_combined_ssva <- sm(combine_de_tables(hs_pairwise_ssva,
                                         excel=paste0("excel/hs_infect_ssva-v", ver, ".xlsx"),
                                         keepers=keepers))

hs_sig_ssva <- sm(extract_significant_genes(hs_combined_ssva,
                                            excel=paste0("excel/hs_infect_ssva_sig-v", ver, ".xlsx")))
hs_sig_ssva$deseq$counts

##              change_counts_up change_counts_down
## sh_vs_uninf              1151                768
## chr_vs_uninf             1312                822
## chr_vs_sh                   0                  0

common_solos_ssva <- subset_significants(hs_sig_ssva)
similar <- sm(compare_de_results(hs_combined_nobatch, hs_combined_ssva, cor_method="spearman"))
similar$limma

## $sh_vs_uninf
## $sh_vs_uninf$logfc
## [1] 0.3367
## 
## $sh_vs_uninf$p
## [1] 0.04093
## 
## $sh_vs_uninf$adjp
## [1] 0.04093
## 
## 
## $chr_vs_uninf
## $chr_vs_uninf$logfc
## [1] 0.4321
## 
## $chr_vs_uninf$p
## [1] 0.08029
## 
## $chr_vs_uninf$adjp
## [1] 0.08029
## 
## 
## $chr_vs_sh
## $chr_vs_sh$logfc
## [1] 0.9379
## 
## $chr_vs_sh$p
## [1] 0.8481
## 
## $chr_vs_sh$adjp
## [1] 0.07285

similar <- sm(compare_de_results(hs_combined_batch, hs_combined_ssva, cor_method="spearman"))
similar$limma

## $sh_vs_uninf
## $sh_vs_uninf$logfc
## [1] 0.3504
## 
## $sh_vs_uninf$p
## [1] 0.1662
## 
## $sh_vs_uninf$adjp
## [1] 0.1662
## 
## 
## $chr_vs_uninf
## $chr_vs_uninf$logfc
## [1] 0.4144
## 
## $chr_vs_uninf$p
## [1] 0.1923
## 
## $chr_vs_uninf$adjp
## [1] 0.1923
## 
## 
## $chr_vs_sh
## $chr_vs_sh$logfc
## [1] 0.9808
## 
## $chr_vs_sh$p
## [1] 0.9256
## 
## $chr_vs_sh$adjp
## [1] 0.9255

ssva_up_venn <- Vennerable::plot(common_solos_ssva$up_venn, doWeights=FALSE)

ssva_down_venn <- Vennerable::plot(common_solos_ssva$down_venn, doWeights=FALSE)

hs_pairwise_fsva <- sm(all_pairwise(hs_uninf_filt, model_batch="fsva"))

hs_combined_fsva <- sm(combine_de_tables(hs_pairwise_fsva,
                                         excel=paste0("excel/hs_infect_fsva-v", ver, ".xlsx"),
                                         keepers=keepers))

hs_sig_fsva <- sm(extract_significant_genes(hs_combined_fsva,
                                            excel=paste0("excel/hs_infect_fsva_sig-v", ver, ".xlsx")))
hs_sig_fsva$deseq$counts

##              change_counts_up change_counts_down
## sh_vs_uninf              1068                328
## chr_vs_uninf             1103                290
## chr_vs_sh                   0                  0

common_solos_fsva <- sm(subset_significants(hs_sig_fsva))
length(common_solos_fsva)

## [1] 10

summary(common_solos_fsva)

##                   Length Class      Mode   
## sh_solo_up        38     data.frame list   
## ch_solo_up        38     data.frame list   
## sh_shared_ch_up   38     data.frame list   
## sh_solo_down      38     data.frame list   
## ch_solo_down      38     data.frame list   
## sh_shared_ch_down 38     data.frame list   
## up_weights         4     -none-     numeric
## down_weights       4     -none-     numeric
## up_venn            1     Venn       S4     
## down_venn          1     Venn       S4

similar <- sm(compare_de_results(hs_combined_nobatch, hs_combined_fsva, cor_method="spearman"))
similar$limma

## $sh_vs_uninf
## $sh_vs_uninf$logfc
## [1] 0.9922
## 
## $sh_vs_uninf$p
## [1] 0.9105
## 
## $sh_vs_uninf$adjp
## [1] 0.9105
## 
## 
## $chr_vs_uninf
## $chr_vs_uninf$logfc
## [1] 0.9865
## 
## $chr_vs_uninf$p
## [1] 0.904
## 
## $chr_vs_uninf$adjp
## [1] 0.904
## 
## 
## $chr_vs_sh
## $chr_vs_sh$logfc
## [1] 0.9401
## 
## $chr_vs_sh$p
## [1] 0.8631
## 
## $chr_vs_sh$adjp
## [1] 0.07755

similar <- sm(compare_de_results(hs_combined_ssva, hs_combined_fsva, cor_method="spearman"))
similar$limma

## $sh_vs_uninf
## $sh_vs_uninf$logfc
## [1] 0.3141
## 
## $sh_vs_uninf$p
## [1] 0.1535
## 
## $sh_vs_uninf$adjp
## [1] 0.1535
## 
## 
## $chr_vs_uninf
## $chr_vs_uninf$logfc
## [1] 0.3568
## 
## $chr_vs_uninf$p
## [1] 0.1718
## 
## $chr_vs_uninf$adjp
## [1] 0.1718
## 
## 
## $chr_vs_sh
## $chr_vs_sh$logfc
## [1] 0.9893
## 
## $chr_vs_sh$p
## [1] 0.9666
## 
## $chr_vs_sh$adjp
## [1] 0.9664

fsva_up_venn <- Vennerable::plot(common_solos_fsva$up_venn, doWeights=FALSE)

fsva_down_venn <- Vennerable::plot(common_solos_fsva$down_venn, doWeights=FALSE)

5 Try with the combat modified data

old_condition <- hs_uninf$design$condition
names(old_condition) <- hs_uninf$design$sampleid
new_condition <- paste0(hs_uninf$design$state, '_', hs_uninf$design$donor)
combat_input <- set_expt_factors(hs_uninf_filt, batch="pathogenstrain", condition=new_condition)
combat_input <- sm(normalize_expt(combat_input, batch="combat_scale"))
combat_input <- set_expt_condition(combat_input, fact=old_condition)

hs_pairwise_combatpath <- sm(all_pairwise(combat_input, model_batch=FALSE, force=TRUE))

hs_combined_combatpath <- sm(combine_de_tables(hs_pairwise_combatpath,
                                               excel=paste0("excel/hs_infect_combatpath-v", ver, ".xlsx"),
                                               keepers=keepers))

hs_sig_combatpath <- sm(extract_significant_genes(hs_combined_combatpath,
                                                  excel=paste0("excel/hs_infect_combatpath_sig-v", ver, ".xlsx")))
hs_sig_combatpath$deseq$counts

##              change_counts_up change_counts_down
## sh_vs_uninf              2055                  4
## chr_vs_uninf             3017               2389
## chr_vs_sh                 672               1129

common_solos_combatpath <- subset_significants(hs_sig_combatpath)
summary(common_solos_combatpath)

##                   Length Class      Mode   
## sh_solo_up        38     data.frame list   
## ch_solo_up        38     data.frame list   
## sh_shared_ch_up   38     data.frame list   
## sh_solo_down      38     data.frame list   
## ch_solo_down      38     data.frame list   
## sh_shared_ch_down 38     data.frame list   
## up_weights         4     -none-     numeric
## down_weights       4     -none-     numeric
## up_venn            1     Venn       S4     
## down_venn          1     Venn       S4

similar <- sm(compare_de_results(hs_combined_nobatch, hs_combined_combatpath, cor_method="spearman"))
similar$limma

## $sh_vs_uninf
## $sh_vs_uninf$logfc
## [1] -0.03092
## 
## $sh_vs_uninf$p
## [1] -0.09456
## 
## $sh_vs_uninf$adjp
## [1] -0.09457
## 
## 
## $chr_vs_uninf
## $chr_vs_uninf$logfc
## [1] -0.1478
## 
## $chr_vs_uninf$p
## [1] -0.04167
## 
## $chr_vs_uninf$adjp
## [1] -0.04167
## 
## 
## $chr_vs_sh
## $chr_vs_sh$logfc
## [1] -0.1485
## 
## $chr_vs_sh$p
## [1] -0.177
## 
## $chr_vs_sh$adjp
## [1] -0.01484

similar <- sm(compare_de_results(hs_combined_fsva, hs_combined_combatpath, cor_method="spearman"))
similar$limma

## $sh_vs_uninf
## $sh_vs_uninf$logfc
## [1] -0.07044
## 
## $sh_vs_uninf$p
## [1] 0.03195
## 
## $sh_vs_uninf$adjp
## [1] 0.03195
## 
## 
## $chr_vs_uninf
## $chr_vs_uninf$logfc
## [1] -0.1593
## 
## $chr_vs_uninf$p
## [1] 0.1013
## 
## $chr_vs_uninf$adjp
## [1] 0.1013
## 
## 
## $chr_vs_sh
## $chr_vs_sh$logfc
## [1] -0.07519
## 
## $chr_vs_sh$p
## [1] -0.02281
## 
## $chr_vs_sh$adjp
## [1] -0.02287

## OUCH!

combat_up_venn <- Vennerable::plot(common_solos_combatpath$up_venn, doWeights=FALSE)

combat_down_venn <- Vennerable::plot(common_solos_combatpath$down_venn, doWeights=FALSE)

## inf_hsstr is the data set that provided the relatively 'pretty' PCA plots in infection_estimation
hs_inf_strbatch$notes

## [1] "New experimental design factors by snp added 2016-09-20Subsetted with experimentname=='infection' on Thu Aug 31 10:46:18 2017.\nSubsetted with condition!='uninf' on Thu Aug 31 10:46:19 2017.\n"

5.1 Compare DE results

For each of the following, perform a simple DE and see what happens: 1. no uninfected strain as batch, try to compare each of the 3 patients chronic/self 2. no uninfected strain as batch, try to compare chronic/self for all 3. and 4. Repeat with uninfected

5.1.1 DE: include uninfected, use strain as batch

The data used in the following is the quantile(cpm(filter())) where the condition was set to a concatenation of patient and healing state, combat was also performed, so we no longer want batch in the experimental model and also we need to pass ‘force=TRUE’ because deseq/edger will need to be coerced into accepting these modified data.

hs_inf$condition

## chr_5430_d108 chr_5397_d108 chr_2504_d108  sh_2272_d108  sh_1022_d108  sh_2189_d108 
##         "chr"         "chr"         "chr"          "sh"          "sh"          "sh" 
## chr_5430_d110 chr_5397_d110 chr_2504_d110  sh_2272_d110  sh_1022_d110  sh_2189_d110 
##         "chr"         "chr"         "chr"          "sh"          "sh"          "sh" 
## chr_5430_d107 chr_5397_d107 chr_2504_d107  sh_2272_d107  sh_1022_d107  sh_2189_d107 
##         "chr"         "chr"         "chr"          "sh"          "sh"          "sh"

## Start by leaving the data relatively alone, especially noting that we do not have a usable batch
## by default.
hs_uninf_filtv2 <- hs_uninf_filt
donor_state <- paste0(hs_uninf_filtv2$design$state, "_", hs_uninf_filtv2$design$donor)
hs_uninf_filtv2 <- set_expt_factors(hs_uninf_filtv2, condition=donor_state)

uninf_patient_keepers <- list(
    "d107_chun" = c("chronic_d107", "uninfected_d107"),
    "d107_shun" = c("self_heal_d107", "uninfected_d107"),
    "d107_chsh" = c("chronic_d107", "self_heal_d107"),
    "d108_chun" = c("chronic_d108", "uninfected_d108"),
    "d108_shun" = c("self_heal_d108", "uninfected_d108"),
    "d108_chsh" = c("chronic_d108", "self_heal_d108"),
    "d110_chun" = c("chronic_d110", "uninfected_d110"),
    "d110_shun" = c("self_heal_d110", "uninfected_d110"),
    "d110_chsh" = c("chronic_d110", "self_heal_d110"))

hs_uninf_filtv2_pairwise <- sm(all_pairwise(hs_uninf_filtv2, model_batch=FALSE))

hs_uninf_filtv2_combined <- sm(combine_de_tables(hs_uninf_filtv2_pairwise,
                                                 keepers=uninf_patient_keepers,
                                                 excel=paste0("excel/hs_infect_patient_nobatch-v", ver, ".xlsx")))

hs_uninf_filtv2_sig <- sm(extract_significant_genes(hs_uninf_filtv2_combined,
                                                    excel=paste0("excel/hs_infect_patient_nobatch_sig-v", ver, ".xlsx")))

##hs_uninf_filtv3_pairwise <- all_pairwise(hs_uninf_filtv2, model_batch="svaseq", surrogates=1)
##hs_uninf_filtv3_combined <- sm(combine_de_tables(hs_uninf_filtv3_pairwise,
##                                                 keepers=uninf_patient_keepers,
##                                                 excel=paste0("excel/hs_infect_patient_fsva-v", ver, ".xlsx")))
##hs_uninf_filtv3_sig <- sm(extract_significant_genes(hs_uninf_filtv3_combined,
##                                                    excel=paste0("excel/hs_infect_patient_fsva_sig-v", ver, ".xlsx")))

5.2 Make some Venns

They want venns, I’ll given them venns…

comp_people <- function(f, s, t) {
    ddf <- f
    ddf[["rn"]] <- rownames(ddf)
    ddf <- ddf[, c("rn", "deseq_logfc")]
    colnames(ddf) <- c("rn", "first")
    ddf <- merge(ddf, s, by="row.names", all=TRUE)
    rownames(ddf) <- ddf[["Row.names"]]
    ddf <- ddf[, -1]
    ddf <- ddf[, c("first", "deseq_logfc")]
    colnames(ddf) <- c("first", "second")
    ddf <- merge(ddf, t, by="row.names", all=TRUE)
    rownames(ddf) <- ddf[["Row.names"]]
    ddf <- ddf[, -1]
    ddf <- ddf[, c("first", "second", "deseq_logfc")]
    colnames(ddf) <- c("first", "second", "third")
    return(ddf)
}

up_sig <- hs_uninf_filtv2_sig$deseq$ups
sh_un_sig <- comp_people(f=up_sig[["self_heal_d107_vs_uninfected_d107"]],
                         s=up_sig[["self_heal_d108_vs_uninfected_d108"]],
                         t=up_sig[["self_heal_d110_vs_uninfected_d110"]])
upa <- sum(!is.na(sh_un_sig[["first"]]) & is.na(sh_un_sig[["second"]]) & is.na(sh_un_sig[["third"]]))
upb <- sum(is.na(sh_un_sig[["first"]]) & !is.na(sh_un_sig[["second"]]) & is.na(sh_un_sig[["third"]]))
upc <- sum(is.na(sh_un_sig[["first"]]) & is.na(sh_un_sig[["second"]]) & !is.na(sh_un_sig[["third"]]))
upab <- sum(!is.na(sh_un_sig[["first"]]) & !is.na(sh_un_sig[["second"]]) & is.na(sh_un_sig[["third"]]))
upbc <- sum(is.na(sh_un_sig[["first"]]) & !is.na(sh_un_sig[["second"]]) & !is.na(sh_un_sig[["third"]]))
upac <- sum(!is.na(sh_un_sig[["first"]]) & is.na(sh_un_sig[["second"]]) & !is.na(sh_un_sig[["third"]]))
upabc <- sum(!is.na(sh_un_sig[["first"]]) & !is.na(sh_un_sig[["second"]]) & !is.na(sh_un_sig[["third"]]))
up_ones <- c("a" = upa, "b" = upb, "c" = upc)
up_twos <- c("a&b" = upab, "b&c" = upbc, "a&c" = upac)
up_threes <- c("a&b&c" <- upabc)
sh_up_venn <- Vennerable::Venn(SetNames = c("a", "b", "c"),
                               Weight = c(0, upa, upb, upc,
                                          upab, upbc, upac,
                                          upabc))
up_res <- Vennerable::plot(sh_up_venn, doWeights=FALSE)

shared_up_sh <- complete.cases(sh_un_sig)
shared_up_sh <- rownames(sh_un_sig[shared_up_sh, ])

de_table_shared_up_sh_first <- hs_uninf_filtv2_combined[["data"]][["self_heal_d107_vs_uninfected_d107"]][shared_up_sh, ]
de_table_shared_up_sh_second <- hs_uninf_filtv2_combined[["data"]][["self_heal_d108_vs_uninfected_d108"]][shared_up_sh, ]
de_table_shared_up_sh_third <- hs_uninf_filtv2_combined[["data"]][["self_heal_d110_vs_uninfected_d110"]][shared_up_sh, ]
de_table_shared_up_sh_all <- merge(de_table_shared_up_sh_first[, c("description", "deseq_logfc", "deseq_adjp")],
                                   de_table_shared_up_sh_second[, c("deseq_logfc", "deseq_adjp")],
                                   by="row.names")
de_table_shared_up_sh_all <- merge(de_table_shared_up_sh_all,
                                   de_table_shared_up_sh_third[, c("deseq_logfc", "deseq_adjp")],
                                   by.x="Row.names", by.y="row.names")
rownames(de_table_shared_up_sh_all) <- de_table_shared_up_sh_all[["Row.names"]]
de_table_shared_up_sh_all <- de_table_shared_up_sh_all[, -1]
colnames(de_table_shared_up_sh_all) <- c("description", "logfc_107", "adjp_107", "logfc_108", "adjp_108", "logfc_110", "adjp_110")
write.csv(de_table_shared_up_sh_all, file="images/de_table_shared_up_sh_all.csv")

up_sig <- hs_uninf_filtv2_sig$deseq$ups
ch_un_sig <- comp_people(f=up_sig[["chronic_d107_vs_uninfected_d107"]],
                         s=up_sig[["chronic_d108_vs_uninfected_d108"]],
                         t=up_sig[["chronic_d110_vs_uninfected_d110"]])
upa <- sum(!is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
upb <- sum(is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
upc <- sum(is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
upab <- sum(!is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
upbc <- sum(is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
upac <- sum(!is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
upabc <- sum(!is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
up_ones <- c("a" = upa, "b" = upb, "c" = upc)
up_twos <- c("a&b" = upab, "b&c" = upbc, "a&c" = upac)
up_threes <- c("a&b&c" <- upabc)
chr_up_venn <- Vennerable::Venn(SetNames = c("a", "b", "c"),
                                Weight = c(0, upa, upb, upc,
                                           upab, upbc, upac,
                                           upabc))
up_res <- Vennerable::plot(chr_up_venn, doWeights=FALSE)

shared_up_ch <- complete.cases(ch_un_sig)
shared_up_ch <- rownames(ch_un_sig[shared_up_ch, ])

de_table_shared_up_ch_first <- hs_uninf_filtv2_combined[["data"]][["chronic_d107_vs_uninfected_d107"]][shared_up_ch, ]
de_table_shared_up_ch_second <- hs_uninf_filtv2_combined[["data"]][["chronic_d108_vs_uninfected_d108"]][shared_up_ch, ]
de_table_shared_up_ch_third <- hs_uninf_filtv2_combined[["data"]][["chronic_d110_vs_uninfected_d110"]][shared_up_ch, ]
de_table_shared_up_ch_all <- merge(de_table_shared_up_ch_first[, c("description", "deseq_logfc", "deseq_adjp")],
                                   de_table_shared_up_ch_second[, c("deseq_logfc", "deseq_adjp")],
                                   by="row.names")
de_table_shared_up_ch_all <- merge(de_table_shared_up_ch_all,
                                   de_table_shared_up_ch_third[, c("deseq_logfc", "deseq_adjp")],
                                   by.x="Row.names", by.y="row.names")
rownames(de_table_shared_up_ch_all) <- de_table_shared_up_ch_all[["Row.names"]]
de_table_shared_up_ch_all <- de_table_shared_up_ch_all[, -1]
colnames(de_table_shared_up_ch_all) <- c("logfc_107", "adjp_107", "logfc_108", "adjp_108", "logfc_110", "adjp_110")
write.csv(de_table_shared_up_ch_all, file="images/de_table_shared_up_ch_all.csv")

down_sig <- hs_uninf_filtv2_sig$deseq$downs
sh_un_sig <- comp_people(f=down_sig[["self_heal_d107_vs_uninfected_d107"]],
                         s=down_sig[["self_heal_d108_vs_uninfected_d108"]],
                         t=down_sig[["self_heal_d110_vs_uninfected_d110"]])
downa <- sum(!is.na(sh_un_sig[["first"]]) & is.na(sh_un_sig[["second"]]) & is.na(sh_un_sig[["third"]]))
downb <- sum(is.na(sh_un_sig[["first"]]) & !is.na(sh_un_sig[["second"]]) & is.na(sh_un_sig[["third"]]))
downc <- sum(is.na(sh_un_sig[["first"]]) & is.na(sh_un_sig[["second"]]) & !is.na(sh_un_sig[["third"]]))
downab <- sum(!is.na(sh_un_sig[["first"]]) & !is.na(sh_un_sig[["second"]]) & is.na(sh_un_sig[["third"]]))
downbc <- sum(is.na(sh_un_sig[["first"]]) & !is.na(sh_un_sig[["second"]]) & !is.na(sh_un_sig[["third"]]))
downac <- sum(!is.na(sh_un_sig[["first"]]) & is.na(sh_un_sig[["second"]]) & !is.na(sh_un_sig[["third"]]))
downabc <- sum(!is.na(sh_un_sig[["first"]]) & !is.na(sh_un_sig[["second"]]) & !is.na(sh_un_sig[["third"]]))
down_ones <- c("a" = downa, "b" = downb, "c" = downc)
down_twos <- c("a&b" = downab, "b&c" = downbc, "a&c" = downac)
down_threes <- c("a&b&c" <- downabc)
sh_down_venn <- Vennerable::Venn(SetNames = c("a", "b", "c"),
                                Weight = c(0, downa, downb, downc,
                                           downab, downbc, downac,
                                           downabc))
sh_down_res <- Vennerable::plot(sh_down_venn, doWeights=FALSE)

shared_down_sh <- complete.cases(sh_un_sig)
shared_down_sh <- rownames(sh_un_sig[shared_down_sh, ])

de_table_shared_down_sh_first <- hs_uninf_filtv2_combined[["data"]][["self_heal_d107_vs_uninfected_d107"]][shared_down_sh, ]
de_table_shared_down_sh_second <- hs_uninf_filtv2_combined[["data"]][["self_heal_d108_vs_uninfected_d108"]][shared_down_sh, ]
de_table_shared_down_sh_third <- hs_uninf_filtv2_combined[["data"]][["self_heal_d110_vs_uninfected_d110"]][shared_down_sh, ]
de_table_shared_down_sh_all <- merge(de_table_shared_down_sh_first[, c("description", "deseq_logfc", "deseq_adjp")],
                                   de_table_shared_down_sh_second[, c("deseq_logfc", "deseq_adjp")],
                                   by="row.names")
de_table_shared_down_sh_all <- merge(de_table_shared_down_sh_all,
                                   de_table_shared_down_sh_third[, c("deseq_logfc", "deseq_adjp")],
                                   by.x="Row.names", by.y="row.names")
rownames(de_table_shared_down_sh_all) <- de_table_shared_down_sh_all[["Row.names"]]
de_table_shared_down_sh_all <- de_table_shared_down_sh_all[, -1]
colnames(de_table_shared_down_sh_all) <- c("description", "logfc_107", "adjp_107", "logfc_108", "adjp_108", "logfc_110", "adjp_110")
write.csv(de_table_shared_down_sh_all, file="images/de_table_shared_down_sh_all.csv")

down_sig <- hs_uninf_filtv2_sig$deseq$downs
ch_un_sig <- comp_people(f=down_sig[["chronic_d107_vs_uninfected_d107"]],
                         s=down_sig[["chronic_d108_vs_uninfected_d108"]],
                         t=down_sig[["chronic_d110_vs_uninfected_d110"]])
downa <- sum(!is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
downb <- sum(is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
downc <- sum(is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
downab <- sum(!is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
downbc <- sum(is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
downac <- sum(!is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
downabc <- sum(!is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
down_ones <- c("a" = downa, "b" = downb, "c" = downc)
down_twos <- c("a&b" = downab, "b&c" = downbc, "a&c" = downac)
down_threes <- c("a&b&c" <- downabc)
chr_down_venn <- Vennerable::Venn(SetNames = c("a", "b", "c"),
                                Weight = c(0, downa, downb, downc,
                                           downab, downbc, downac,
                                           downabc))
chr_down_res <- Vennerable::plot(chr_down_venn, doWeights=FALSE)

shared_down_ch <- complete.cases(ch_un_sig)
shared_down_ch <- rownames(ch_un_sig[shared_down_ch, ])

de_table_shared_down_ch_first <- hs_uninf_filtv2_combined[["data"]][["chronic_d107_vs_uninfected_d107"]][shared_down_ch, ]
de_table_shared_down_ch_second <- hs_uninf_filtv2_combined[["data"]][["chronic_d108_vs_uninfected_d108"]][shared_down_ch, ]
de_table_shared_down_ch_third <- hs_uninf_filtv2_combined[["data"]][["chronic_d110_vs_uninfected_d110"]][shared_down_ch, ]
de_table_shared_down_ch_all <- merge(de_table_shared_down_ch_first[, c("description", "deseq_logfc", "deseq_adjp")],
                                   de_table_shared_down_ch_second[, c("deseq_logfc", "deseq_adjp")],
                                   by="row.names")
de_table_shared_down_ch_all <- merge(de_table_shared_down_ch_all,
                                   de_table_shared_down_ch_third[, c("deseq_logfc", "deseq_adjp")],
                                   by.x="Row.names", by.y="row.names")
rownames(de_table_shared_down_ch_all) <- de_table_shared_down_ch_all[["Row.names"]]
de_table_shared_down_ch_all <- de_table_shared_down_ch_all[, -1]
colnames(de_table_shared_down_ch_all) <- c("logfc_107", "adjp_107", "logfc_108", "adjp_108", "logfc_110", "adjp_110")
write.csv(de_table_shared_down_ch_all, file="images/de_table_shared_down_ch_all.csv")

upupgenes <- intersect(shared_up_sh, shared_up_ch)
upup <- length(upupgenes)
upsh_notch <- !shared_up_sh %in% shared_up_ch
upsh_notch <- shared_up_sh[upsh_notch]
upch_notsh <- !shared_up_ch %in% shared_up_sh
upch_notsh <- shared_up_ch[upch_notsh]
## upnot and length(upsh_notch) should be equivalent.
upnot <- sum(! shared_up_sh %in% shared_up_ch)
notup <- sum(! shared_up_ch %in% shared_up_sh)
shared_up <- Vennerable::Venn(SetNames=c("up_sh", "up_ch"),
                              Weight=c(0, upnot, notup, upup))
Vennerable::plot(shared_up, doWeights=FALSE)

downdowngenes <- intersect(shared_down_sh, shared_down_ch)
downdown <- length(downdowngenes)
downsh_notch <- !shared_down_sh %in% shared_down_ch
downsh_notch <- shared_down_sh[downsh_notch]
downch_notsh <- !shared_down_ch %in% shared_down_sh
downch_notsh <- shared_down_ch[downch_notsh]
## downnot and length(downsh_notch) should be equivalent.
downnot <- sum(! shared_down_sh %in% shared_down_ch)
notdown <- sum(! shared_down_ch %in% shared_down_sh)
shared_down <- Vennerable::Venn(SetNames=c("down_sh", "down_ch"),
                              Weight=c(0, downnot, notdown, downdown))
Vennerable::plot(shared_down, doWeights=FALSE)

kept_columns <- c("transcriptid", "geneid", "description", "deseq_logfc", "deseq_adjp")
xls_result <- write_xls(data=up_sig[["self_heal_d107_vs_uninfected_d107"]][upupgenes, kept_columns],
                        sheet="upsh_upch")
xls_result <- write_xls(data=up_sig[["self_heal_d107_vs_uninfected_d107"]][upsh_notch, kept_columns],
                        sheet="upsh_noch",
                        wb=xls_result[["workbook"]])
xls_result <- write_xls(data=up_sig[["chronic_d107_vs_uninfected_d107"]][upch_notsh, kept_columns],
                        sheet="upch_nosh", wb=xls_result[["workbook"]])
xls_result <- write_xls(data=down_sig[["self_heal_d107_vs_uninfected_d107"]][downdowngenes, kept_columns],
                        sheet="downsh_downch", wb=xls_result[["workbook"]])
xls_result <- write_xls(data=down_sig[["self_heal_d107_vs_uninfected_d107"]][downsh_notch, kept_columns],
                        sheet="downsh_noch", wb=xls_result[["workbook"]])
xls_result <- write_xls(data=down_sig[["chronic_d107_vs_uninfected_d107"]][downch_notsh, kept_columns],
                        sheet="downch_nosh", wb=xls_result[["workbook"]],
                        excel=paste0("excel/figure_5c_stuff-v", ver, ".xlsx"))

## Saving to: excel/figure_5c_stuff-v20170820.xlsx

up_sig <- hs_uninf_filtv2_sig$deseq$ups
ch_un_sig <- comp_people(f=up_sig[["self_heal_d107_vs_uninfected_d107"]],
                         s=up_sig[["self_heal_d108_vs_uninfected_d108"]],
                         t=up_sig[["self_heal_d110_vs_uninfected_d110"]])
upa <- sum(!is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
upb <- sum(is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
upc <- sum(is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
upab <- sum(!is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
upbc <- sum(is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
upac <- sum(!is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
upabc <- sum(!is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
up_ones <- c("a" = upa, "b" = upb, "c" = upc)
up_twos <- c("a&b" = upab, "b&c" = upbc, "a&c" = upac)
up_threes <- c("a&b&c" <- upabc)
up_venn <- Vennerable::Venn(SetNames = c("a", "b", "c"),
                                Weight = c(0, upa, upb, upc,
                                           upab, upbc, upac,
                                           upabc))
up_res <- Vennerable::plot(up_venn, doWeights=FALSE)

down_sig <- hs_uninf_filtv2_sig$deseq$downs
ch_un_sig <- comp_people(f=down_sig[["self_heal_d107_vs_uninfected_d107"]],
                         s=down_sig[["self_heal_d108_vs_uninfected_d108"]],
                         t=down_sig[["self_heal_d110_vs_uninfected_d110"]])
downa <- sum(!is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
downb <- sum(is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
downc <- sum(is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
downab <- sum(!is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
downbc <- sum(is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
downac <- sum(!is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
downabc <- sum(!is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
down_ones <- c("a" = downa, "b" = downb, "c" = downc)
down_twos <- c("a&b" = downab, "b&c" = downbc, "a&c" = downac)
down_threes <- c("a&b&c" <- downabc)
down_venn <- Vennerable::Venn(SetNames = c("a", "b", "c"),
                                Weight = c(0, downa, downb, downc,
                                           downab, downbc, downac,
                                           downabc))
down_res <- Vennerable::plot(down_venn, doWeights=FALSE)

up_sig <- hs_uninf_filtv2_sig$deseq$ups
ch_un_sig <- comp_people(f=up_sig[["chronic_d107_vs_self_heal_d107"]],
                         s=up_sig[["chronic_d108_vs_self_heal_d108"]],
                         t=up_sig[["chronic_d110_vs_self_heal_d110"]])
upa <- sum(!is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
upb <- sum(is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
upc <- sum(is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
upab <- sum(!is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
upbc <- sum(is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
upac <- sum(!is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
upabc <- sum(!is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
up_ones <- c("a" = upa, "b" = upb, "c" = upc)
up_twos <- c("a&b" = upab, "b&c" = upbc, "a&c" = upac)
up_threes <- c("a&b&c" <- upabc)
up_venn <- Vennerable::Venn(SetNames = c("a", "b", "c"),
                                Weight = c(0, upa, upb, upc,
                                           upab, upbc, upac,
                                           upabc))
up_res <- Vennerable::plot(up_venn, doWeights=FALSE)

down_sig <- hs_uninf_filtv2_sig$deseq$downs
ch_un_sig <- comp_people(f=down_sig[["chronic_d107_vs_self_heal_d107"]],
                         s=down_sig[["chronic_d108_vs_self_heal_d108"]],
                         t=down_sig[["chronic_d110_vs_self_heal_d110"]])
downa <- sum(!is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
downb <- sum(is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
downc <- sum(is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
downab <- sum(!is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
downbc <- sum(is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
downac <- sum(!is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
downabc <- sum(!is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
down_ones <- c("a" = downa, "b" = downb, "c" = downc)
down_twos <- c("a&b" = downab, "b&c" = downbc, "a&c" = downac)
down_threes <- c("a&b&c" <- downabc)
down_venn <- Vennerable::Venn(SetNames = c("a", "b", "c"),
                                Weight = c(0, downa, downb, downc,
                                           downab, downbc, downac,
                                           downabc))
down_res <- Vennerable::plot(down_venn, doWeights=FALSE)

up_sig <- hs_uninf_filtv3_sig$deseq$ups
ch_un_sig <- comp_people(f=up_sig[["chronic_d107_vs_uninfected_d107"]],
                         s=up_sig[["chronic_d108_vs_uninfected_d108"]],
                         t=up_sig[["chronic_d110_vs_uninfected_d110"]])
upa <- sum(!is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
upb <- sum(is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
upc <- sum(is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
upab <- sum(!is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
upbc <- sum(is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
upac <- sum(!is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
upabc <- sum(!is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
up_ones <- c("a" = upa, "b" = upb, "c" = upc)
up_twos <- c("a&b" = upab, "b&c" = upbc, "a&c" = upac)
up_threes <- c("a&b&c" <- upabc)
up_venn <- Vennerable::Venn(SetNames = c("a", "b", "c"),
                                Weight = c(0, upa, upb, upc,
                                           upab, upbc, upac,
                                           upabc))
up_res <- Vennerable::plot(up_venn, doWeights=FALSE)

down_sig <- hs_uninf_filtv3_sig$deseq$downs
ch_un_sig <- comp_people(f=down_sig[["chronic_d107_vs_uninfected_d107"]],
                         s=down_sig[["chronic_d108_vs_uninfected_d108"]],
                         t=down_sig[["chronic_d110_vs_uninfected_d110"]])
downa <- sum(!is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
downb <- sum(is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
downc <- sum(is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
downab <- sum(!is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & is.na(ch_un_sig[["third"]]))
downbc <- sum(is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
downac <- sum(!is.na(ch_un_sig[["first"]]) & is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
downabc <- sum(!is.na(ch_un_sig[["first"]]) & !is.na(ch_un_sig[["second"]]) & !is.na(ch_un_sig[["third"]]))
down_ones <- c("a" = downa, "b" = downb, "c" = downc)
down_twos <- c("a&b" = downab, "b&c" = downbc, "a&c" = downac)
down_threes <- c("a&b&c" <- downabc)
down_venn <- Vennerable::Venn(SetNames = c("a", "b", "c"),
                                Weight = c(0, downa, downb, downc,
                                           downab, downbc, downac,
                                           downabc))
down_res <- Vennerable::plot(down_venn, doWeights=FALSE)

5.2.1 DE: include uninfected, repeat with condition simplified to chronic/self

Do I think the exact same thing as in the previous comparison, but now simplify the ‘condition’ factor to just self-healing vs. chronic and see what happens.

uninf_strainbatch_qcf <- set_expt_batch(expt=uninf_cqf, fact="pathogenstrain")
uninf_strainbatch_qcf <- set_expt_condition(expt=uninf_strainbatch_qcf, fact="state")
withuninf_strainbatch_pairs_chsh <- all_pairwise(uninf_strainbatch_qcf,
                                                 model_batch=FALSE, force=TRUE)
chsh_keepers <- list(
    "ch_sh" = c("chronic", "self_heal"),
    "ch_nil" = c("chronic", "uninfected"),
    "sh_nil" = c("self_heal", "uninfected"))
withuninf_strainbatch_tables_chsh <- sm(combine_de_tables(withuninf_strainbatch_pairs_chsh,
                                                          keepers=chsh_keepers,
                                                          excel=paste0("excel/withuninf_strainbatch_chsh_pairwise-v", ver, ".xlsx")))
withuninf_strainbatch_sig_chsh <- extract_significant_genes(withuninf_strainbatch_tables_chsh,
                                                            excel=paste0("excel/withuninf_strainbatch_chsh_sig-v", ver, ".xlsx"))
##withuninf_strainbatch_tables_chsh$limma_plots
##withuninf_strainbatch_tables_chsh$edger_plots
##withuninf_strainbatch_tables_chsh$deseq_plots

6 Figure 5

Generate DE lists of each donor for all contrasts for PBMCs.

1. Venn sh/uninf vs chr/uninf 2 venn diagram up. (donor in model)
2. Venn sh/uninf vs chr/uninf 2 venn diagram down.
3. Venn Sh/uninf up genes 3 venn diagram.
4. Venn Sh/uninf down genes 3 venn.
5. Venn Chr/uninf up genes 3 venn.
6. Venn Chr/uninf down genes 3 venn.
7. 2 way venn of (common up in 3 venn sh/uninf) vs. (common up in 3 venn chr/uninf)
8. 2 way venn of (common down in 3 venn sh/uninf) vs. (common down in 3 venn chr/uninf)

pp(file="images/fig_05a-sh_uninf_vs_chr_uninf_up.pdf")

## NULL

Vennerable::plot(common_solos_batch$up_venn, doWeights=FALSE)
dev.off()

## png 
##   2

pp(file="images/fig_05b-sh_uninf_vs_chr_uninf_down.pdf")

## NULL

Vennerable::plot(common_solos_batch$down_venn, doWeights=FALSE)
dev.off()

## png 
##   2

pp(file="images/fig_05c-sh_uninf_donors_up.pdf")

## NULL

Vennerable::plot(sh_up_venn, doWeights=FALSE)
dev.off()

## png 
##   2

pp(file="images/fig_05d-sh_uninf_donors_down.pdf")

## NULL

Vennerable::plot(sh_down_venn, doWeights=FALSE)
dev.off()

## png 
##   2

pp(file="images/fig_05e-chr_uninf_donors_up.pdf")

## NULL

Vennerable::plot(chr_up_venn, doWeights=FALSE)
dev.off()

## png 
##   2

pp(file="images/fig_05f-chr_uninf_donors_down.pdf")

## NULL

Vennerable::plot(chr_down_venn, doWeights=FALSE)
dev.off()

## png 
##   2

pp(file="images/fig_05g-up_common.pdf")

## NULL

Vennerable::plot(shared_up, doWeights=FALSE)
dev.off()

## png 
##   2

pp(file="images/fig_05h-down_common.pdf")

## NULL

Vennerable::plot(shared_down, doWeights=FALSE)
dev.off()

## png 
##   2

7 Table Szzz

Tables of the stuff in figure 5

funkytown <- hpgltools:::compare_logfc_plots(withuninf_strainbatch_tables_chsh)
funkytown$ch_nil$le
funkytown$ch_nil$ld
funkytown$ch_nil$de

8 Try again on the parasite data

8.1 Remember our data set

lp_inf_filt <- sm(normalize_expt(lp_inf, filter=TRUE))

lp_pairwise_nobatch <- sm(all_pairwise(lp_inf_filt, model_batch=FALSE))
lp_combined_nobatch <- sm(combine_de_tables(lp_pairwise_nobatch,
                                            excel=paste0("excel/lp_infect_nobatch-v", ver, ".xlsx")))

lp_sig_nobatch <- sm(extract_significant_genes(lp_combined_nobatch,
                                               excel=paste0("excel/lp_infect_nobatch_sig-v", ver, ".xlsx")))

lp_pairwise_batch <- sm(all_pairwise(lp_inf_filt, model_batch=TRUE))
lp_combined_batch <- sm(combine_de_tables(lp_pairwise_batch,
                                          excel=paste0("excel/lp_infect_batch-v", ver, ".xlsx")))

lp_sig_batch <- sm(extract_significant_genes(lp_combined_batch,
                                             excel=paste0("excel/lp_infect_batch_sig-v", ver, ".xlsx")))

lp_pairwise_ssva <- sm(all_pairwise(lp_inf_filt, model_batch="ssva"))
lp_combined_ssva <- sm(combine_de_tables(lp_pairwise_ssva,
                                         excel=paste0("excel/lp_infect_ssva-v", ver, ".xlsx")))

lp_sig_ssva <- sm(extract_significant_genes(lp_combined_ssva,
                                             excel=paste0("excel/lp_infect_ssva_sig-v", ver, ".xlsx")))

lp_pairwise_fsva <- sm(all_pairwise(lp_inf_filt, model_batch="fsva"))
lp_combined_fsva <- sm(combine_de_tables(lp_pairwise_fsva,
                                         excel=paste0("excel/lp_infect_fsva-v", ver, ".xlsx")))

lp_sig_fsva <- sm(extract_significant_genes(lp_combined_fsva,
                                             excel=paste0("excel/lp_infect_fsva_sig-v", ver, ".xlsx")))

index.html 02_infection_estimation.html

pander::pander(sessionInfo())

R version 3.4.1 (2017-06-30)

**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: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: ruv(v.0.9.6), Vennerable(v.3.1.0.9000) and hpgltools(v.2017.01)

loaded via a namespace (and not attached): minqa(v.1.2.4), colorspace(v.1.3-2), colorRamps(v.2.3), rprojroot(v.1.2), htmlTable(v.1.9), corpcor(v.1.6.9), XVector(v.0.16.0), GenomicRanges(v.1.28.4), base64enc(v.0.1-3), roxygen2(v.6.0.1), ggrepel(v.0.6.5), bit64(v.0.9-7), AnnotationDbi(v.1.38.2), xml2(v.1.1.1), codetools(v.0.2-15), splines(v.3.4.1), doParallel(v.1.0.10), geneplotter(v.1.54.0), knitr(v.1.17), Formula(v.1.2-2), nloptr(v.1.0.4), Rsamtools(v.1.28.0), pbkrtest(v.0.4-7), annotate(v.1.54.0), cluster(v.2.0.6), graph(v.1.54.0), compiler(v.3.4.1), backports(v.1.1.0), Matrix(v.1.2-11), lazyeval(v.0.2.0), limma(v.3.32.5), acepack(v.1.4.1), htmltools(v.0.3.6), tools(v.3.4.1), gtable(v.0.2.0), GenomeInfoDbData(v.0.99.0), reshape2(v.1.4.2), Rcpp(v.0.12.12), Biobase(v.2.36.2), Biostrings(v.2.44.2), gdata(v.2.18.0), preprocessCore(v.1.38.1), nlme(v.3.1-131), rtracklayer(v.1.36.4), iterators(v.1.0.8), stringr(v.1.2.0), openxlsx(v.4.0.17), testthat(v.1.0.2), lme4(v.1.1-13), gtools(v.3.5.0), devtools(v.1.13.3), XML(v.3.98-1.9), edgeR(v.3.18.1), directlabels(v.2017.03.31), zlibbioc(v.1.22.0), MASS(v.7.3-47), scales(v.0.5.0), BiocInstaller(v.1.26.0), parallel(v.3.4.1), SummarizedExperiment(v.1.6.3), RBGL(v.1.52.0), RColorBrewer(v.1.1-2), yaml(v.2.1.14), memoise(v.1.1.0), gridExtra(v.2.2.1), pander(v.0.6.1), ggplot2(v.2.2.1), biomaRt(v.2.32.1), rpart(v.4.1-11), latticeExtra(v.0.6-28), stringi(v.1.1.5), RSQLite(v.2.0), genefilter(v.1.58.1), S4Vectors(v.0.14.3), foreach(v.1.4.3), checkmate(v.1.8.3), GenomicFeatures(v.1.28.4), caTools(v.1.17.1), BiocGenerics(v.0.22.0), BiocParallel(v.1.10.1), GenomeInfoDb(v.1.12.2), rlang(v.0.1.2), commonmark(v.1.2), matrixStats(v.0.52.2), bitops(v.1.0-6), evaluate(v.0.10.1), lattice(v.0.20-35), GenomicAlignments(v.1.12.2), htmlwidgets(v.0.9), labeling(v.0.3), bit(v.1.1-12), plyr(v.1.8.4), magrittr(v.1.5), variancePartition(v.1.6.0), DESeq2(v.1.16.1), R6(v.2.2.2), IRanges(v.2.10.2), gplots(v.3.0.1), Hmisc(v.4.0-3), DelayedArray(v.0.2.7), DBI(v.0.7), mgcv(v.1.8-19), foreign(v.0.8-69), withr(v.2.0.0), survival(v.2.41-3), RCurl(v.1.95-4.8), nnet(v.7.3-12), tibble(v.1.3.4), crayon(v.1.3.2), KernSmooth(v.2.23-15), OrganismDbi(v.1.18.0), rmarkdown(v.1.6), locfit(v.1.5-9.1), grid(v.3.4.1), sva(v.3.24.4), data.table(v.1.10.4), blob(v.1.1.0), digest(v.0.6.12), xtable(v.1.8-2), stats4(v.3.4.1), munsell(v.0.4.3) and quadprog(v.1.5-5)

message(paste0("This is hpgltools commit: ", get_git_commit()))

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

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

## > git reset 739e4a10345a89efb17b37e7de07c5811491ccea

## R> packrat::restore()

## This is hpgltools commit: Thu Aug 31 11:24:06 2017 -0400: 739e4a10345a89efb17b37e7de07c5811491ccea

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

## Saving to 03_expression_infection-v20170820.rda.xz

tmp <- sm(saveme(filename=this_save))