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.
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.
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.
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).
## There were 18, now there are 6 samples.
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
## Error in combine_de_tables(d107_pairwise): object 'd107_pairwise' not found
## Error in summary(d107_table$data): object 'd107_table' not found
## Error in extract_significant_genes(d107_table, according_to = "deseq", : object 'd107_table' not found
## Error in head(d107_sig$deseq$ups[[1]]): object 'd107_sig' not found
## Error in extract_de_plots(d107_pairwise, p = 0.1, p_type = "raw"): object 'd107_pairwise' not found
## Error in eval(expr, envir, enclos): object 'd107_ma' not found
## Error in plot_histogram(d107_table[["data"]][["sh_vs_chr"]][, c("deseq_p")]): object 'd107_table' not found
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
## There were 18, now there are 6 samples.
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
## Error in combine_de_tables(d108_pairwise): object 'd108_pairwise' not found
## Error in summary(d108_table$data): object 'd108_table' not found
## Error in extract_significant_genes(d108_table, according_to = "deseq", : object 'd108_table' not found
## Error in head(d108_sig$deseq$ups[[1]]): object 'd108_sig' not found
## Error in extract_de_plots(d108_pairwise, p = 0.1, p_type = "raw"): object 'd108_pairwise' not found
## Error in eval(expr, envir, enclos): object 'd108_ma' not found
## Error in plot_histogram(d108_table[["data"]][["sh_vs_chr"]][, c("deseq_p")]): object 'd108_table' not found
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
## There were 18, now there are 6 samples.
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
## Error in combine_de_tables(d110_pairwise): object 'd110_pairwise' not found
## Error in summary(d110_table$data): object 'd110_table' not found
## Error in extract_significant_genes(d110_table, according_to = "deseq", : object 'd110_table' not found
## Error in head(d110_sig$deseq$ups[[1]]): object 'd110_sig' not found
## Error in extract_de_plots(d110_pairwise, p = 0.1, p_type = "raw"): object 'd110_pairwise' not found
## Error in eval(expr, envir, enclos): object 'd110_ma' not found
## Error in plot_histogram(d110_table[["data"]][["sh_vs_chr"]][, c("deseq_p")]): object 'd110_table' not found
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
hs_tmp <- set_expt_batches(hs_inf_filt, fact="pathogenstrain")
hs_tmp2 <- normalize_expt(hs_tmp, transform="log2", convert="cpm", norm="quant", filter=TRUE)
## This function will replace the expt$expressionset slot with:
## log2(cpm(quant(hpgl(data))))
## It backs up the current data into a slot named:
## expt$backup_expressionset. It will also save copies of each step along the way
## in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
## when invoking limma. The appropriate libsize is the 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: hpgl
## Removing 0 low-count genes (11944 remaining).
## Step 2: normalizing the data with quant.
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
## Error in plot_pca(hs_tmp2): object 'hs_tmp2' not found
## There were 18, now there are 6 samples.
## Error in set_expt_conditions(d107, fact = "sh", ids = ): The provided factor is not in the design matrix.
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
## Error in combine_de_tables(d107_pairwise): object 'd107_pairwise' not found
## Error in summary(d107_table$data): object 'd107_table' not found
## Error in extract_significant_genes(d107_table, according_to = "deseq", : object 'd107_table' not found
## Error in head(d107_sig$deseq$ups[[1]]): object 'd107_sig' not found
## Error in extract_de_plots(d107_pairwise, p = 0.1, p_type = "raw"): object 'd107_pairwise' not found
## Error in eval(expr, envir, enclos): object 'd107_ma' not found
## Error in plot_histogram(d107_table[["data"]][["sh_vs_chr"]][, c("deseq_p")]): object 'd107_table' not found
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
## There were 18, now there are 6 samples.
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
## Error in combine_de_tables(d108_pairwise): object 'd108_pairwise' not found
## Error in summary(d108_table$data): object 'd108_table' not found
## Error in extract_significant_genes(d108_table, according_to = "deseq", : object 'd108_table' not found
## Error in head(d108_sig$deseq$ups[[1]]): object 'd108_sig' not found
## Error in extract_de_plots(d108_pairwise, p = 0.1, p_type = "raw"): object 'd108_pairwise' not found
## Error in eval(expr, envir, enclos): object 'd108_ma' not found
## Error in plot_histogram(d108_table[["data"]][["sh_vs_chr"]][, c("deseq_p")]): object 'd108_table' not found
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
## There were 18, now there are 6 samples.
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
## Error in combine_de_tables(d110_pairwise): object 'd110_pairwise' not found
## Error in summary(d110_table$data): object 'd110_table' not found
## Error in extract_significant_genes(d110_table, according_to = "deseq", : object 'd110_table' not found
## Error in head(d110_sig$deseq$ups[[1]]): object 'd110_sig' not found
## Error in extract_de_plots(d110_pairwise, p = 0.1, p_type = "raw"): object 'd110_pairwise' not found
## Error in eval(expr, envir, enclos): object 'd110_ma' not found
## Error in plot_histogram(d110_table[["data"]][["sh_vs_chr"]][, c("deseq_p")]): object 'd110_table' not found
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
## There were 18, now there are 15 samples.
## This function will replace the expt$expressionset slot with:
## hpgl(data)
## It backs up the current data into a slot named:
## expt$backup_expressionset. It will also save copies of each step along the way
## in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
## when invoking limma. The appropriate libsize is the non-log(cpm(normalized)).
## This is most likely kept at:
## 'new_expt$normalized$intermediate_counts$normalization$libsizes'
## A copy of this may also be found at:
## new_expt$best_libsize
## Leaving the data in its current base format, keep in mind that
## some metrics are easier to see when the data is log2 transformed, but
## EdgeR/DESeq do not accept transformed data.
## Leaving the data unconverted. It is often advisable to cpm/rpkm
## the data to normalize for sampling differences, keep in mind though that rpkm
## has some annoying biases, and voom() by default does a cpm (though hpgl_voom()
## will try to detect this).
## Leaving the data unnormalized. This is necessary for DESeq, but
## EdgeR/limma might benefit from normalization. Good choices include quantile,
## size-factor, tmm, etc.
## 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: hpgl
## Removing 38 low-count genes (11906 remaining).
## Step 2: not normalizing the data.
## Step 3: not converting the data.
## Step 4: not transforming the data.
## Step 5: not doing batch correction.
## There were 18, now there are 15 samples.
## This function will replace the expt$expressionset slot with:
## hpgl(data)
## It backs up the current data into a slot named:
## expt$backup_expressionset. It will also save copies of each step along the way
## in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
## when invoking limma. The appropriate libsize is the non-log(cpm(normalized)).
## This is most likely kept at:
## 'new_expt$normalized$intermediate_counts$normalization$libsizes'
## A copy of this may also be found at:
## new_expt$best_libsize
## Leaving the data in its current base format, keep in mind that
## some metrics are easier to see when the data is log2 transformed, but
## EdgeR/DESeq do not accept transformed data.
## Leaving the data unconverted. It is often advisable to cpm/rpkm
## the data to normalize for sampling differences, keep in mind though that rpkm
## has some annoying biases, and voom() by default does a cpm (though hpgl_voom()
## will try to detect this).
## Leaving the data unnormalized. This is necessary for DESeq, but
## EdgeR/limma might benefit from normalization. Good choices include quantile,
## size-factor, tmm, etc.
## 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: hpgl
## Removing 51 low-count genes (11893 remaining).
## Step 2: not normalizing the data.
## Step 3: not converting the data.
## Step 4: not transforming the data.
## Step 5: not doing batch correction.
Once using batch in model, once with svaseq.
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
excel_file <- glue::glue("excel/{rundate}_hs_infect_removeone_batchmodel_contr-v{ver}.xlsx")
remove_one_table <- combine_de_tables(remove_one_de, excel=excel_file,
keepers=keepers_inf)
## Error in combine_de_tables(remove_one_de, excel = excel_file, keepers = keepers_inf): object 'remove_one_de' not found
excel_file <- glue::glue("excel/{rundate}_hs_infect_removeone_batchmodel_sig-v{ver}.xlsx")
remove_one_sig <- extract_significant_genes(remove_one_table, excel=excel_file,
according_to="deseq", p=0.1, lfc=0.6)
## Error in extract_significant_genes(remove_one_table, excel = excel_file, : object 'remove_one_table' not found
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
excel_file <- glue::glue("excel/{rundate}_hs_infect_removeone_svaseq_contr-v{ver}.xlsx")
remove_one_table_sva <- combine_de_tables(remove_one_de_sva, excel=excel_file,
keepers=keepers_inf)
## Error in combine_de_tables(remove_one_de_sva, excel = excel_file, keepers = keepers_inf): object 'remove_one_de_sva' not found
excel_file <- glue::glue("excel/{rundate}_hs_infect_removeone_svaseq_sig-v{ver}.xlsx")
remove_one_sig_sva <- extract_significant_genes(remove_one_table_sva, excel=excel_file,
according_to="deseq", p=0.1, lfc=0.6)
## Error in extract_significant_genes(remove_one_table_sva, excel = excel_file, : object 'remove_one_table_sva' not found
Repeat as above, removing both weirdo samples.
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
excel_file <- glue::glue("excel/{rundate}_hs_infect_removeone_batchmodel_contr-v{ver}.xlsx")
remove_two_table <- combine_de_tables(remove_two_de, excel=excel_file,
keepers=keepers_inf)
## Error in combine_de_tables(remove_two_de, excel = excel_file, keepers = keepers_inf): object 'remove_two_de' not found
excel_file <- glue::glue("excel/{rundate}_hs_infect_removeone_batchmodel_sig-v{ver}.xlsx")
remove_two_sig <- extract_significant_genes(remove_two_table, excel=excel_file,
according_to="deseq", p=0.1, lfc=0.6)
## Error in extract_significant_genes(remove_two_table, excel = excel_file, : object 'remove_two_table' not found
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
excel_file <- glue::glue("excel/{rundate}_hs_infect_removeone_svaseq_contr-v{ver}.xlsx")
remove_two_table_sva <- combine_de_tables(remove_two_de_sva, excel=excel_file,
keepers=keepers_inf)
## Error in combine_de_tables(remove_two_de_sva, excel = excel_file, keepers = keepers_inf): object 'remove_two_de_sva' not found
excel_file <- glue::glue("excel/{rundate}_hs_infect_removeone_svaseq_sig-v{ver}.xlsx")
remove_two_sig_sva <- extract_significant_genes(remove_two_table_sva, excel=excel_file,
according_to="deseq", p=0.1, lfc=0.6)
## Error in extract_significant_genes(remove_two_table_sva, excel = excel_file, : object 'remove_two_table_sva' not found
Try building a donor*state interaction model.
I think I realized my confusion in this: the only interaction models I have used before were a experimental factor*experimental floating-point observation. If it is the case, that the following models are all x vs. factor_reference, then I guess it makes sense; but if not, then I have no clue what is going on in these.
First make sure that self-healing is the reference.
Then make some designs of potential interest.
data <- edgeR::DGEList(exprs(hs_inf_filt))
data <- edgeR::calcNormFactors(data)
design <- pData(hs_inf_filt)
design[["condition"]] <- relevel(design[["condition"]], ref="sh")
design[["donor"]] <- as.factor(design[["donor"]])
design[["donor"]] <- relevel(design[["donor"]], ref="d107")
chsh <- model.matrix(object=~0+condition, data=design)
head(chsh) ## Looks good to me.
## conditionsh conditionchr
## chr_5430_d108 0 1
## chr_5397_d108 0 1
## chr_2504_d108 0 1
## sh_2272_d108 1 0
## sh_1022_d108 1 0
## sh_2189_d108 1 0
## donord107 donord108 donord110
## chr_5430_d108 0 1 0
## chr_5397_d108 0 1 0
## chr_2504_d108 0 1 0
## sh_2272_d108 0 1 0
## sh_1022_d108 0 1 0
## sh_2189_d108 0 1 0
cond_donor_inter <- model.matrix(object=~0+condition+donor+condition:donor, data=design)
head(cond_donor_inter)
## conditionsh conditionchr donord108 donord110
## chr_5430_d108 0 1 1 0
## chr_5397_d108 0 1 1 0
## chr_2504_d108 0 1 1 0
## sh_2272_d108 1 0 1 0
## sh_1022_d108 1 0 1 0
## sh_2189_d108 1 0 1 0
## conditionchr:donord108 conditionchr:donord110
## chr_5430_d108 1 0
## chr_5397_d108 1 0
## chr_2504_d108 1 0
## sh_2272_d108 0 0
## sh_1022_d108 0 0
## sh_2189_d108 0 0
## Why are there no conditionchr:donor107 conditionsh:donor107 columns?
## I know that if I relevel the donor factor to set the reference to another donor, that changes
## but this still does not make sense to me.
## I think the above model should be the same as:
cond_donor_interv2 <- model.matrix(object=~0+condition*donor, data=design)
head(cond_donor_interv2)
## conditionsh conditionchr donord108 donord110
## chr_5430_d108 0 1 1 0
## chr_5397_d108 0 1 1 0
## chr_2504_d108 0 1 1 0
## sh_2272_d108 1 0 1 0
## sh_1022_d108 1 0 1 0
## sh_2189_d108 1 0 1 0
## conditionchr:donord108 conditionchr:donord110
## chr_5430_d108 1 0
## chr_5397_d108 1 0
## chr_2504_d108 1 0
## sh_2272_d108 0 0
## sh_1022_d108 0 0
## sh_2189_d108 0 0
testthat::expect_equal(cond_donor_inter, cond_donor_interv2)
## So as I understand it, if I do a test of coef='conditionchr', this is looking at
## chr vs. sh
## If I do coef='conditionchr:donor110' this is the donor effect for chr/sh for donor110?
## If this is true, then how do I get this for donor107?
## Perhaps my understanding is backwards?
donor_cond_inter <- model.matrix(object=~donor+condition+donor:condition, data=design)
donor_cond_interv2 <- model.matrix(object=~donor*condition, data=design)
testthat::expect_equal(donor_cond_inter, donor_cond_interv2)
head(donor_cond_inter)
## (Intercept) donord108 donord110 conditionchr
## chr_5430_d108 1 1 0 1
## chr_5397_d108 1 1 0 1
## chr_2504_d108 1 1 0 1
## sh_2272_d108 1 1 0 0
## sh_1022_d108 1 1 0 0
## sh_2189_d108 1 1 0 0
## donord108:conditionchr donord110:conditionchr
## chr_5430_d108 1 0
## chr_5397_d108 1 0
## chr_2504_d108 1 0
## sh_2272_d108 0 0
## sh_1022_d108 0 0
## sh_2189_d108 0 0
disp <- edgeR::estimateDisp(data, design=cond_donor_inter)
fit <- edgeR::glmQLFit(disp, cond_donor_inter)
## I think this is my global chr/sh
chr_sh_qlf <- edgeR::glmQLFTest(fit, coef="conditionchr")
chr_sh_result <- edgeR::topTags(chr_sh_qlf)
summary(chr_sh_result$table)
## logFC logCPM F PValue
## Min. :-14.9 Min. :5.04 Min. :21787 Min. :2.58e-28
## 1st Qu.:-14.3 1st Qu.:5.70 1st Qu.:22428 1st Qu.:7.45e-28
## Median :-13.6 Median :6.35 Median :23358 Median :1.01e-27
## Mean :-13.8 Mean :6.12 Mean :23693 Mean :1.05e-27
## 3rd Qu.:-13.4 3rd Qu.:6.51 3rd Qu.:24216 3rd Qu.:1.42e-27
## Max. :-13.1 Max. :6.87 Max. :27493 Max. :1.80e-27
## FDR
## Min. :8.93e-25
## 1st Qu.:8.93e-25
## Median :8.93e-25
## Mean :8.93e-25
## 3rd Qu.:8.93e-25
## Max. :8.93e-25
## [1] "conditionchr:donord108"
d110_ch_qlf <- edgeR::glmQLFTest(fit, coef="conditionchr:donord110")
d110_ch_result <- edgeR::topTags(d110_ch_qlf)
summary(d110_ch_result$table)
## logFC logCPM F PValue
## Min. :-1.305 Min. :0.762 Min. :10.6 Min. :0.000729
## 1st Qu.:-0.659 1st Qu.:2.330 1st Qu.:10.9 1st Qu.:0.001731
## Median :-0.484 Median :3.884 Median :12.9 Median :0.002407
## Mean :-0.393 Mean :3.428 Mean :13.0 Mean :0.002769
## 3rd Qu.:-0.395 3rd Qu.:4.252 3rd Qu.:13.9 3rd Qu.:0.004282
## Max. : 1.133 Max. :5.837 Max. :17.0 Max. :0.004813
## FDR
## Min. :1
## 1st Qu.:1
## Median :1
## Mean :1
## 3rd Qu.:1
## Max. :1
summarized <- DESeq2::DESeqDataSetFromMatrix(countData=exprs(hs_inf_filt),
colData=design,
design=~condition*donor)
## converting counts to integer mode
dataset <- DESeq2::DESeqDataSet(se=summarized, design=~condition*donor)
run <- DESeq2::DESeq(dataset)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## [1] "Intercept" "condition_chr_vs_sh"
## [3] "donor_d108_vs_d107" "donor_d110_vs_d107"
## [5] "conditionchr.donord108" "conditionchr.donord110"
## Length Class Mode
## 6 DESeqResults S4
## log2 fold change (MLE): condition chr vs sh
## Wald test p-value: condition chr vs sh
## DataFrame with 6 rows and 6 columns
## baseMean log2FoldChange lfcSE
## <numeric> <numeric> <numeric>
## ENSG00000000419 214.886619915588 -0.0934303786250254 0.12300393077645
## ENSG00000000457 318.115512380986 -0.068815027486115 0.120178565849457
## ENSG00000000460 102.031345099905 0.197578912145072 0.203833151567596
## ENSG00000000938 5676.81912456484 0.0269256146826752 0.197204584725487
## ENSG00000000971 53.5624271330181 0.344010769796456 0.293375111918845
## ENSG00000001036 460.076692835547 0.0156360632261453 0.148896600722685
## stat pvalue padj
## <numeric> <numeric> <numeric>
## ENSG00000000419 -0.759572300131027 0.447510281800118 0.997804083570582
## ENSG00000000457 -0.572606496006257 0.566911160357569 0.997804083570582
## ENSG00000000460 0.969316868358142 0.332387115533936 0.997804083570582
## ENSG00000000938 0.136536453856568 0.891397208363037 0.997804083570582
## ENSG00000000971 1.17259697847723 0.240957461669005 0.997804083570582
## ENSG00000001036 0.105012895863667 0.916365575881328 0.997804083570582
## Maybe my thinking is just a little wrong: is this the interaction of chr/sh with d110/d107?
## That would make more sense I think.
res_donor110_chrsh <- DESeq2::results(run, name="conditionchr.donord110")
summary(res_donor110_chrsh)
## Length Class Mode
## 6 DESeqResults S4
## log2 fold change (MLE): conditionchr.donord110
## Wald test p-value: conditionchr.donord110
## DataFrame with 6 rows and 6 columns
## baseMean log2FoldChange lfcSE
## <numeric> <numeric> <numeric>
## ENSG00000000419 214.886619915588 0.295162897632013 0.172516639204641
## ENSG00000000457 318.115512380986 0.108194027119009 0.172637547612825
## ENSG00000000460 102.031345099905 0.376710409795487 0.293516835681214
## ENSG00000000938 5676.81912456484 0.0530822519196353 0.278697424287239
## ENSG00000000971 53.5624271330181 -0.388785003193931 0.436191950074646
## ENSG00000001036 460.076692835547 0.10657602357023 0.206872518783626
## stat pvalue padj
## <numeric> <numeric> <numeric>
## ENSG00000000419 1.71092422732562 0.0870951015388544 0.999936246663048
## ENSG00000000457 0.626712025368067 0.530848019781775 0.999936246663048
## ENSG00000000460 1.2834371456792 0.199338967059447 0.999936246663048
## ENSG00000000938 0.19046552746367 0.848944353796083 0.999936246663048
## ENSG00000000971 -0.891316318715644 0.372759496477825 0.999936246663048
## ENSG00000001036 0.515177289844383 0.606429137146659 0.999936246663048
res_donor108_chrsh <- DESeq2::results(run, name="conditionchr.donord108")
summary(res_donor110_chrsh)
## Length Class Mode
## 6 DESeqResults S4
## log2 fold change (MLE): conditionchr.donord110
## Wald test p-value: conditionchr.donord110
## DataFrame with 6 rows and 6 columns
## baseMean log2FoldChange lfcSE
## <numeric> <numeric> <numeric>
## ENSG00000000419 214.886619915588 0.295162897632013 0.172516639204641
## ENSG00000000457 318.115512380986 0.108194027119009 0.172637547612825
## ENSG00000000460 102.031345099905 0.376710409795487 0.293516835681214
## ENSG00000000938 5676.81912456484 0.0530822519196353 0.278697424287239
## ENSG00000000971 53.5624271330181 -0.388785003193931 0.436191950074646
## ENSG00000001036 460.076692835547 0.10657602357023 0.206872518783626
## stat pvalue padj
## <numeric> <numeric> <numeric>
## ENSG00000000419 1.71092422732562 0.0870951015388544 0.999936246663048
## ENSG00000000457 0.626712025368067 0.530848019781775 0.999936246663048
## ENSG00000000460 1.2834371456792 0.199338967059447 0.999936246663048
## ENSG00000000938 0.19046552746367 0.848944353796083 0.999936246663048
## ENSG00000000971 -0.891316318715644 0.372759496477825 0.999936246663048
## ENSG00000001036 0.515177289844383 0.606429137146659 0.999936246663048
## Length Class Mode
## 6 DESeqResults S4
## log2 fold change (MLE): donor d110 vs d107
## Wald test p-value: donor d110 vs d107
## DataFrame with 6 rows and 6 columns
## baseMean log2FoldChange lfcSE
## <numeric> <numeric> <numeric>
## ENSG00000000419 214.886619915588 -0.113766010123853 0.120483955620923
## ENSG00000000457 318.115512380986 -0.456556848180769 0.121122399221616
## ENSG00000000460 102.031345099905 -0.5010401182924 0.210579406027488
## ENSG00000000938 5676.81912456484 0.505694298463597 0.197026331645792
## ENSG00000000971 53.5624271330181 -0.842424805273729 0.308759375189539
## ENSG00000001036 460.076692835547 0.568390991546202 0.145482651751354
## stat pvalue
## <numeric> <numeric>
## ENSG00000000419 -0.944241990873823 0.345046001807488
## ENSG00000000457 -3.76938411982257 0.000163650866259048
## ENSG00000000460 -2.3793405430491 0.0173436449623682
## ENSG00000000938 2.56663307000365 0.0102691214578726
## ENSG00000000971 -2.72841854520721 0.00636388048892608
## ENSG00000001036 3.90693312710334 9.3475010176665e-05
## padj
## <numeric>
## ENSG00000000419 0.476062375176819
## ENSG00000000457 0.000739090391073861
## ENSG00000000460 0.0422654009174772
## ENSG00000000938 0.0271711897470546
## ENSG00000000971 0.0181060512075926
## ENSG00000001036 0.000448806683671478
## Comp.1 Comp.2
## [1,] "31.7:ENSG00000129824" "22.16:ENSG00000100079"
## [2,] "30.31:ENSG00000012817" "16.86:ENSG00000128283"
## [3,] "30:ENSG00000185736" "15.92:ENSG00000170439"
## [4,] "29.65:ENSG00000067048" "13.44:ENSG00000156113"
## [5,] "27.68:ENSG00000114374" "13.39:ENSG00000007350"
## [6,] "25.19:ENSG00000183878" "13.25:ENSG00000123219"
## [7,] "21.54:ENSG00000198692" "13.11:ENSG00000205809"
## [8,] "20.86:ENSG00000067646" "13.01:ENSG00000165810"
## [9,] "18.84:ENSG00000074410" "12.6:ENSG00000150687"
## [10,] "16.59:ENSG00000179344" "12.6:ENSG00000167880"
## [11,] "15.01:ENSG00000134184" "12.36:ENSG00000150510"
## [12,] "14.9:ENSG00000196735" "12.1:ENSG00000167094"
## [13,] "14.46:ENSG00000122641" "12.09:ENSG00000146918"
## [14,] "14.37:ENSG00000121380" "11.77:ENSG00000173372"
## [15,] "14.37:ENSG00000161055" "11.67:ENSG00000214944"
## [16,] "13.88:ENSG00000154620" "11.34:ENSG00000174945"
## [17,] "13.76:ENSG00000166923" "11.01:ENSG00000240247"
## [18,] "13.5:ENSG00000123453" "10.99:ENSG00000205810"
## [19,] "13.06:ENSG00000123454" "10.95:ENSG00000149635"
## [20,] "13.02:ENSG00000108688" "10.72:ENSG00000102837"
## [21,] "12.49:ENSG00000131435" "10.69:ENSG00000109956"
## [22,] "12.47:ENSG00000136449" "10.68:ENSG00000061337"
## [23,] "12.22:ENSG00000110092" "10.6:ENSG00000206047"
## [24,] "12.05:ENSG00000120645" "10.51:ENSG00000205846"
## [25,] "11.83:ENSG00000205057" "10.5:ENSG00000183542"
## [26,] "11.54:ENSG00000123689" "10.04:ENSG00000113721"
## [27,] "11.34:ENSG00000170160" "9.932:ENSG00000101188"
## [28,] "11.32:ENSG00000113302" "9.836:ENSG00000113657"
## [29,] "11.2:ENSG00000102970" "9.783:ENSG00000196421"
## [30,] "11.08:ENSG00000198203" "9.572:ENSG00000145708"
## [31,] "11.08:ENSG00000184371" "9.241:ENSG00000133962"
## [32,] "11.07:ENSG00000149635" "9.213:ENSG00000148488"
## [33,] "10.94:ENSG00000204644" "8.953:ENSG00000125780"
## [34,] "10.83:ENSG00000050165" "8.802:ENSG00000243772"
## [35,] "10.66:ENSG00000150510" "8.796:ENSG00000170160"
## [36,] "10.47:ENSG00000182871" "8.76:ENSG00000162383"
## [37,] "10.42:ENSG00000134668" "8.719:ENSG00000108370"
## [38,] "10.23:ENSG00000164530" "8.705:ENSG00000054793"
## [39,] "10.22:ENSG00000178860" "8.564:ENSG00000168329"
## [40,] "10.19:ENSG00000109471" "8.503:ENSG00000184349"
## [41,] "10.17:ENSG00000197121" "8.377:ENSG00000100450"
## [42,] "10.15:ENSG00000164400" "8.224:ENSG00000173239"
## [43,] "9.967:ENSG00000180616" "8.088:ENSG00000086548"
## [44,] "9.874:ENSG00000069482" "8.001:ENSG00000111199"
## [45,] "9.8:ENSG00000133687" "7.934:ENSG00000239839"
## [46,] "9.786:ENSG00000134716" "7.928:ENSG00000139567"
## [47,] "9.778:ENSG00000115009" "7.858:ENSG00000188820"
## [48,] "9.737:ENSG00000197822" "7.744:ENSG00000101425"
## [49,] "9.705:ENSG00000172724" "7.624:ENSG00000176083"
## [50,] "9.686:ENSG00000181634" "7.616:ENSG00000169908"
annot <- fData(hs_inf_filt)
comp1_genes <- unlist(strsplit(x=as.character(scores$highest[, 1]), split=":"))
idx <- grep(pattern="^ENSG", x=comp1_genes)
comp1_genes <- comp1_genes[idx]
comp1_high_scores <- exclude_genes_expt(hs_inf_filt, ids=comp1_genes, method="keep")
## Before removal, there were 11944 entries.
## Now there are 50 entries.
## Percent kept: 0.140, 0.119, 0.075, 0.131, 0.083, 0.119, 0.017, 0.016, 0.009, 0.019, 0.010, 0.012, 0.210, 0.213, 0.183, 0.252, 0.162, 0.221
## Percent removed: 99.860, 99.881, 99.925, 99.869, 99.917, 99.881, 99.983, 99.984, 99.991, 99.981, 99.990, 99.988, 99.790, 99.787, 99.817, 99.748, 99.838, 99.779
plot_sample_heatmap(
data=sm(normalize_expt(comp1_high_scores, transform="log2", convert="cpm", norm="quant")),
row_label=NULL)
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
ensembl_gene_id | description | |
---|---|---|
ENSG00000129824 | ENSG00000129824 | ribosomal protein S4, Y-linked 1 [Source:HGNC Symbol;Acc:HGNC:10425] |
ENSG00000012817 | ENSG00000012817 | lysine demethylase 5D [Source:HGNC Symbol;Acc:HGNC:11115] |
ENSG00000185736 | ENSG00000185736 | adenosine deaminase, RNA specific B2 (inactive) [Source:HGNC Symbol;Acc:HGNC:227] |
ENSG00000067048 | ENSG00000067048 | DEAD-box helicase 3, Y-linked [Source:HGNC Symbol;Acc:HGNC:2699] |
ENSG00000114374 | ENSG00000114374 | ubiquitin specific peptidase 9, Y-linked [Source:HGNC Symbol;Acc:HGNC:12633] |
ENSG00000183878 | ENSG00000183878 | ubiquitously transcribed tetratricopeptide repeat containing, Y-linked [Source:HGNC Symbol;Acc:HGNC:12638] |
ENSG00000198692 | ENSG00000198692 | eukaryotic translation initiation factor 1A, Y-linked [Source:HGNC Symbol;Acc:HGNC:3252] |
ENSG00000067646 | ENSG00000067646 | zinc finger protein, Y-linked [Source:HGNC Symbol;Acc:HGNC:12870] |
ENSG00000074410 | ENSG00000074410 | carbonic anhydrase 12 [Source:HGNC Symbol;Acc:HGNC:1371] |
ENSG00000179344 | ENSG00000179344 | major histocompatibility complex, class II, DQ beta 1 [Source:HGNC Symbol;Acc:HGNC:4944] |
ENSG00000134184 | ENSG00000134184 | glutathione S-transferase mu 1 [Source:HGNC Symbol;Acc:HGNC:4632] |
ENSG00000196735 | ENSG00000196735 | major histocompatibility complex, class II, DQ alpha 1 [Source:HGNC Symbol;Acc:HGNC:4942] |
ENSG00000122641 | ENSG00000122641 | inhibin beta A subunit [Source:HGNC Symbol;Acc:HGNC:6066] |
ENSG00000121380 | ENSG00000121380 | BCL2 like 14 [Source:HGNC Symbol;Acc:HGNC:16657] |
ENSG00000161055 | ENSG00000161055 | secretoglobin family 3A member 1 [Source:HGNC Symbol;Acc:HGNC:18384] |
ENSG00000154620 | ENSG00000154620 | thymosin beta 4, Y-linked [Source:HGNC Symbol;Acc:HGNC:11882] |
ENSG00000166923 | ENSG00000166923 | gremlin 1, DAN family BMP antagonist [Source:HGNC Symbol;Acc:HGNC:2001] |
ENSG00000123453 | ENSG00000123453 | sarcosine dehydrogenase [Source:HGNC Symbol;Acc:HGNC:10536] |
ENSG00000123454 | ENSG00000123454 | dopamine beta-hydroxylase [Source:HGNC Symbol;Acc:HGNC:2689] |
ENSG00000108688 | ENSG00000108688 | C-C motif chemokine ligand 7 [Source:HGNC Symbol;Acc:HGNC:10634] |
ENSG00000131435 | ENSG00000131435 | PDZ and LIM domain 4 [Source:HGNC Symbol;Acc:HGNC:16501] |
ENSG00000136449 | ENSG00000136449 | MYCBP associated protein [Source:HGNC Symbol;Acc:HGNC:19677] |
ENSG00000110092 | ENSG00000110092 | cyclin D1 [Source:HGNC Symbol;Acc:HGNC:1582] |
ENSG00000120645 | ENSG00000120645 | IQ motif and Sec7 domain 3 [Source:HGNC Symbol;Acc:HGNC:29193] |
ENSG00000205057 | ENSG00000205057 | chronic lymphocytic leukemia up-regulated 1 opposite strand [Source:HGNC Symbol;Acc:HGNC:24070] |
ENSG00000123689 | ENSG00000123689 | G0/G1 switch 2 [Source:HGNC Symbol;Acc:HGNC:30229] |
ENSG00000170160 | ENSG00000170160 | coiled-coil domain containing 144A [Source:HGNC Symbol;Acc:HGNC:29072] |
ENSG00000113302 | ENSG00000113302 | interleukin 12B [Source:HGNC Symbol;Acc:HGNC:5970] |
ENSG00000102970 | ENSG00000102970 | C-C motif chemokine ligand 17 [Source:HGNC Symbol;Acc:HGNC:10615] |
ENSG00000198203 | ENSG00000198203 | sulfotransferase family 1C member 2 [Source:HGNC Symbol;Acc:HGNC:11456] |
ENSG00000184371 | ENSG00000184371 | colony stimulating factor 1 [Source:HGNC Symbol;Acc:HGNC:2432] |
ENSG00000149635 | ENSG00000149635 | osteoclast stimulatory transmembrane protein [Source:HGNC Symbol;Acc:HGNC:16116] |
ENSG00000204644 | ENSG00000204644 | ZFP57 zinc finger protein [Source:HGNC Symbol;Acc:HGNC:18791] |
ENSG00000050165 | ENSG00000050165 | dickkopf WNT signaling pathway inhibitor 3 [Source:HGNC Symbol;Acc:HGNC:2893] |
ENSG00000150510 | ENSG00000150510 | family with sequence similarity 124 member A [Source:HGNC Symbol;Acc:HGNC:26413] |
ENSG00000182871 | ENSG00000182871 | collagen type XVIII alpha 1 chain [Source:HGNC Symbol;Acc:HGNC:2195] |
ENSG00000134668 | ENSG00000134668 | SPOC domain containing 1 [Source:HGNC Symbol;Acc:HGNC:26338] |
ENSG00000164530 | ENSG00000164530 | peptidase inhibitor 16 [Source:HGNC Symbol;Acc:HGNC:21245] |
ENSG00000178860 | ENSG00000178860 | musculin [Source:HGNC Symbol;Acc:HGNC:7321] |
ENSG00000109471 | ENSG00000109471 | interleukin 2 [Source:HGNC Symbol;Acc:HGNC:6001] |
ENSG00000197121 | ENSG00000197121 | post-GPI attachment to proteins 1 [Source:HGNC Symbol;Acc:HGNC:25712] |
ENSG00000164400 | ENSG00000164400 | colony stimulating factor 2 [Source:HGNC Symbol;Acc:HGNC:2434] |
ENSG00000180616 | ENSG00000180616 | somatostatin receptor 2 [Source:HGNC Symbol;Acc:HGNC:11331] |
ENSG00000069482 | ENSG00000069482 | galanin and GMAP prepropeptide [Source:HGNC Symbol;Acc:HGNC:4114] |
ENSG00000133687 | ENSG00000133687 | transmembrane and tetratricopeptide repeat containing 1 [Source:HGNC Symbol;Acc:HGNC:24099] |
ENSG00000134716 | ENSG00000134716 | cytochrome P450 family 2 subfamily J member 2 [Source:HGNC Symbol;Acc:HGNC:2634] |
ENSG00000115009 | ENSG00000115009 | C-C motif chemokine ligand 20 [Source:HGNC Symbol;Acc:HGNC:10619] |
ENSG00000197822 | ENSG00000197822 | occludin [Source:HGNC Symbol;Acc:HGNC:8104] |
ENSG00000172724 | ENSG00000172724 | C-C motif chemokine ligand 19 [Source:HGNC Symbol;Acc:HGNC:10617] |
ENSG00000181634 | ENSG00000181634 | tumor necrosis factor superfamily member 15 [Source:HGNC Symbol;Acc:HGNC:11931] |
comp2_genes <- unlist(strsplit(x=as.character(scores$highest[, 2]), split=":"))
idx <- grep(pattern="^ENSG", x=comp2_genes)
comp2_genes <- comp2_genes[idx]
comp2_high_scores <- exclude_genes_expt(hs_inf_filt, ids=comp2_genes, method="keep")
## Before removal, there were 11944 entries.
## Now there are 50 entries.
## Percent kept: 0.097, 0.087, 0.077, 0.094, 0.077, 0.088, 0.053, 0.050, 0.059, 0.055, 0.063, 0.057, 0.024, 0.017, 0.015, 0.021, 0.017, 0.017
## Percent removed: 99.903, 99.913, 99.923, 99.906, 99.923, 99.912, 99.947, 99.950, 99.941, 99.945, 99.937, 99.943, 99.976, 99.983, 99.985, 99.979, 99.983, 99.983
plot_sample_heatmap(
data=sm(normalize_expt(comp2_high_scores, transform="log2", convert="cpm", norm="quant")),
row_label=NULL)
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
ensembl_gene_id | description | |
---|---|---|
ENSG00000100079 | ENSG00000100079 | galectin 2 [Source:HGNC Symbol;Acc:HGNC:6562] |
ENSG00000128283 | ENSG00000128283 | CDC42 effector protein 1 [Source:HGNC Symbol;Acc:HGNC:17014] |
ENSG00000170439 | ENSG00000170439 | methyltransferase like 7B [Source:HGNC Symbol;Acc:HGNC:28276] |
ENSG00000156113 | ENSG00000156113 | potassium calcium-activated channel subfamily M alpha 1 [Source:HGNC Symbol;Acc:HGNC:6284] |
ENSG00000007350 | ENSG00000007350 | transketolase like 1 [Source:HGNC Symbol;Acc:HGNC:11835] |
ENSG00000123219 | ENSG00000123219 | centromere protein K [Source:HGNC Symbol;Acc:HGNC:29479] |
ENSG00000205809 | ENSG00000205809 | killer cell lectin like receptor C2 [Source:HGNC Symbol;Acc:HGNC:6375] |
ENSG00000165810 | ENSG00000165810 | butyrophilin like 9 [Source:HGNC Symbol;Acc:HGNC:24176] |
ENSG00000150687 | ENSG00000150687 | protease, serine 23 [Source:HGNC Symbol;Acc:HGNC:14370] |
ENSG00000167880 | ENSG00000167880 | envoplakin [Source:HGNC Symbol;Acc:HGNC:3503] |
ENSG00000150510 | ENSG00000150510 | family with sequence similarity 124 member A [Source:HGNC Symbol;Acc:HGNC:26413] |
ENSG00000167094 | ENSG00000167094 | tetratricopeptide repeat domain 16 [Source:HGNC Symbol;Acc:HGNC:26536] |
ENSG00000146918 | ENSG00000146918 | non-SMC condensin II complex subunit G2 [Source:HGNC Symbol;Acc:HGNC:21904] |
ENSG00000173372 | ENSG00000173372 | complement C1q A chain [Source:HGNC Symbol;Acc:HGNC:1241] |
ENSG00000214944 | ENSG00000214944 | Rho guanine nucleotide exchange factor 28 [Source:HGNC Symbol;Acc:HGNC:30322] |
ENSG00000174945 | ENSG00000174945 | archaelysin family metallopeptidase 1 [Source:HGNC Symbol;Acc:HGNC:22231] |
ENSG00000240247 | ENSG00000240247 | defensin alpha 1B [Source:HGNC Symbol;Acc:HGNC:33596] |
ENSG00000205810 | ENSG00000205810 | killer cell lectin like receptor C3 [Source:HGNC Symbol;Acc:HGNC:6376] |
ENSG00000149635 | ENSG00000149635 | osteoclast stimulatory transmembrane protein [Source:HGNC Symbol;Acc:HGNC:16116] |
ENSG00000102837 | ENSG00000102837 | olfactomedin 4 [Source:HGNC Symbol;Acc:HGNC:17190] |
ENSG00000109956 | ENSG00000109956 | beta-1,3-glucuronyltransferase 1 [Source:HGNC Symbol;Acc:HGNC:921] |
ENSG00000061337 | ENSG00000061337 | leucine zipper tumor suppressor 1 [Source:HGNC Symbol;Acc:HGNC:13861] |
ENSG00000206047 | ENSG00000206047 | defensin alpha 1 [Source:HGNC Symbol;Acc:HGNC:2761] |
ENSG00000205846 | ENSG00000205846 | C-type lectin domain family 6 member A [Source:HGNC Symbol;Acc:HGNC:14556] |
ENSG00000183542 | ENSG00000183542 | killer cell lectin like receptor C4 [Source:HGNC Symbol;Acc:HGNC:6377] |
ENSG00000113721 | ENSG00000113721 | platelet derived growth factor receptor beta [Source:HGNC Symbol;Acc:HGNC:8804] |
ENSG00000101188 | ENSG00000101188 | neurotensin receptor 1 [Source:HGNC Symbol;Acc:HGNC:8039] |
ENSG00000113657 | ENSG00000113657 | dihydropyrimidinase like 3 [Source:HGNC Symbol;Acc:HGNC:3015] |
ENSG00000196421 | ENSG00000196421 | long intergenic non-protein coding RNA 176 [Source:HGNC Symbol;Acc:HGNC:27655] |
ENSG00000145708 | ENSG00000145708 | corticotropin releasing hormone binding protein [Source:HGNC Symbol;Acc:HGNC:2356] |
ENSG00000133962 | ENSG00000133962 | cation channel sperm associated auxiliary subunit beta [Source:HGNC Symbol;Acc:HGNC:20500] |
ENSG00000148488 | ENSG00000148488 | ST8 alpha-N-acetyl-neuraminide alpha-2,8-sialyltransferase 6 [Source:HGNC Symbol;Acc:HGNC:23317] |
ENSG00000125780 | ENSG00000125780 | transglutaminase 3 [Source:HGNC Symbol;Acc:HGNC:11779] |
ENSG00000243772 | ENSG00000243772 | killer cell immunoglobulin like receptor, two Ig domains and long cytoplasmic tail 3 [Source:HGNC Symbol;Acc:HGNC:6331] |
ENSG00000170160 | ENSG00000170160 | coiled-coil domain containing 144A [Source:HGNC Symbol;Acc:HGNC:29072] |
ENSG00000162383 | ENSG00000162383 | solute carrier family 1 member 7 [Source:HGNC Symbol;Acc:HGNC:10945] |
ENSG00000108370 | ENSG00000108370 | regulator of G-protein signaling 9 [Source:HGNC Symbol;Acc:HGNC:10004] |
ENSG00000054793 | ENSG00000054793 | ATPase phospholipid transporting 9A (putative) [Source:HGNC Symbol;Acc:HGNC:13540] |
ENSG00000168329 | ENSG00000168329 | C-X3-C motif chemokine receptor 1 [Source:HGNC Symbol;Acc:HGNC:2558] |
ENSG00000184349 | ENSG00000184349 | ephrin A5 [Source:HGNC Symbol;Acc:HGNC:3225] |
ENSG00000100450 | ENSG00000100450 | granzyme H [Source:HGNC Symbol;Acc:HGNC:4710] |
ENSG00000173239 | ENSG00000173239 | lipase family member M [Source:HGNC Symbol;Acc:HGNC:23455] |
ENSG00000086548 | ENSG00000086548 | carcinoembryonic antigen related cell adhesion molecule 6 [Source:HGNC Symbol;Acc:HGNC:1818] |
ENSG00000111199 | ENSG00000111199 | transient receptor potential cation channel subfamily V member 4 [Source:HGNC Symbol;Acc:HGNC:18083] |
ENSG00000239839 | ENSG00000239839 | defensin alpha 3 [Source:HGNC Symbol;Acc:HGNC:2762] |
ENSG00000139567 | ENSG00000139567 | activin A receptor like type 1 [Source:HGNC Symbol;Acc:HGNC:175] |
ENSG00000188820 | ENSG00000188820 | family with sequence similarity 26 member F [Source:HGNC Symbol;Acc:HGNC:33391] |
ENSG00000101425 | ENSG00000101425 | bactericidal/permeability-increasing protein [Source:HGNC Symbol;Acc:HGNC:1095] |
ENSG00000176083 | ENSG00000176083 | zinc finger protein 683 [Source:HGNC Symbol;Acc:HGNC:28495] |
ENSG00000169908 | ENSG00000169908 | transmembrane 4 L six family member 1 [Source:HGNC Symbol;Acc:HGNC:11853] |
comp3_genes <- unlist(strsplit(x=as.character(scores$highest[, 3]), split=":"))
idx <- grep(pattern="^ENSG", x=comp3_genes)
comp3_genes <- comp3_genes[idx]
comp3_high_scores <- exclude_genes_expt(hs_inf_filt, ids=comp3_genes, method="keep")
## Before removal, there were 11944 entries.
## Now there are 50 entries.
## Percent kept: 0.401, 0.316, 0.140, 0.384, 0.155, 0.370, 0.275, 0.242, 0.094, 0.218, 0.125, 0.158, 0.486, 0.395, 0.215, 0.493, 0.271, 0.316
## Percent removed: 99.599, 99.684, 99.860, 99.616, 99.845, 99.630, 99.725, 99.758, 99.906, 99.782, 99.875, 99.842, 99.514, 99.605, 99.785, 99.507, 99.729, 99.684
plot_sample_heatmap(
data=sm(normalize_expt(comp3_high_scores, transform="log2", convert="cpm", norm="quant")),
row_label=NULL)
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
ensembl_gene_id | description | |
---|---|---|
ENSG00000137673 | ENSG00000137673 | matrix metallopeptidase 7 [Source:HGNC Symbol;Acc:HGNC:7174] |
ENSG00000113302 | ENSG00000113302 | interleukin 12B [Source:HGNC Symbol;Acc:HGNC:5970] |
ENSG00000169908 | ENSG00000169908 | transmembrane 4 L six family member 1 [Source:HGNC Symbol;Acc:HGNC:11853] |
ENSG00000074410 | ENSG00000074410 | carbonic anhydrase 12 [Source:HGNC Symbol;Acc:HGNC:1371] |
ENSG00000122641 | ENSG00000122641 | inhibin beta A subunit [Source:HGNC Symbol;Acc:HGNC:6066] |
ENSG00000166923 | ENSG00000166923 | gremlin 1, DAN family BMP antagonist [Source:HGNC Symbol;Acc:HGNC:2001] |
ENSG00000106178 | ENSG00000106178 | C-C motif chemokine ligand 24 [Source:HGNC Symbol;Acc:HGNC:10623] |
ENSG00000125144 | ENSG00000125144 | metallothionein 1G [Source:HGNC Symbol;Acc:HGNC:7399] |
ENSG00000197506 | ENSG00000197506 | solute carrier family 28 member 3 [Source:HGNC Symbol;Acc:HGNC:16484] |
ENSG00000136695 | ENSG00000136695 | interleukin 36 receptor antagonist [Source:HGNC Symbol;Acc:HGNC:15561] |
ENSG00000205362 | ENSG00000205362 | metallothionein 1A [Source:HGNC Symbol;Acc:HGNC:7393] |
ENSG00000187848 | ENSG00000187848 | purinergic receptor P2X 2 [Source:HGNC Symbol;Acc:HGNC:15459] |
ENSG00000111052 | ENSG00000111052 | lin-7 homolog A, crumbs cell polarity complex component [Source:HGNC Symbol;Acc:HGNC:17787] |
ENSG00000069482 | ENSG00000069482 | galanin and GMAP prepropeptide [Source:HGNC Symbol;Acc:HGNC:4114] |
ENSG00000113657 | ENSG00000113657 | dihydropyrimidinase like 3 [Source:HGNC Symbol;Acc:HGNC:3015] |
ENSG00000186081 | ENSG00000186081 | keratin 5 [Source:HGNC Symbol;Acc:HGNC:6442] |
ENSG00000148677 | ENSG00000148677 | ankyrin repeat domain 1 [Source:HGNC Symbol;Acc:HGNC:15819] |
ENSG00000187134 | ENSG00000187134 | aldo-keto reductase family 1 member C1 [Source:HGNC Symbol;Acc:HGNC:384] |
ENSG00000105976 | ENSG00000105976 | MET proto-oncogene, receptor tyrosine kinase [Source:HGNC Symbol;Acc:HGNC:7029] |
ENSG00000105855 | ENSG00000105855 | integrin subunit beta 8 [Source:HGNC Symbol;Acc:HGNC:6163] |
ENSG00000172724 | ENSG00000172724 | C-C motif chemokine ligand 19 [Source:HGNC Symbol;Acc:HGNC:10617] |
ENSG00000133048 | ENSG00000133048 | chitinase 3 like 1 [Source:HGNC Symbol;Acc:HGNC:1932] |
ENSG00000100234 | ENSG00000100234 | TIMP metallopeptidase inhibitor 3 [Source:HGNC Symbol;Acc:HGNC:11822] |
ENSG00000138772 | ENSG00000138772 | annexin A3 [Source:HGNC Symbol;Acc:HGNC:541] |
ENSG00000108688 | ENSG00000108688 | C-C motif chemokine ligand 7 [Source:HGNC Symbol;Acc:HGNC:10634] |
ENSG00000151012 | ENSG00000151012 | solute carrier family 7 member 11 [Source:HGNC Symbol;Acc:HGNC:11059] |
ENSG00000117594 | ENSG00000117594 | hydroxysteroid 11-beta dehydrogenase 1 [Source:HGNC Symbol;Acc:HGNC:5208] |
ENSG00000039560 | ENSG00000039560 | retinoic acid induced 14 [Source:HGNC Symbol;Acc:HGNC:14873] |
ENSG00000118785 | ENSG00000118785 | secreted phosphoprotein 1 [Source:HGNC Symbol;Acc:HGNC:11255] |
ENSG00000154856 | ENSG00000154856 | APC down-regulated 1 [Source:HGNC Symbol;Acc:HGNC:15718] |
ENSG00000110436 | ENSG00000110436 | solute carrier family 1 member 2 [Source:HGNC Symbol;Acc:HGNC:10940] |
ENSG00000050730 | ENSG00000050730 | TNFAIP3 interacting protein 3 [Source:HGNC Symbol;Acc:HGNC:19315] |
ENSG00000176076 | ENSG00000176076 | potassium voltage-gated channel subfamily E regulatory subunit 5 [Source:HGNC Symbol;Acc:HGNC:6241] |
ENSG00000205358 | ENSG00000205358 | metallothionein 1H [Source:HGNC Symbol;Acc:HGNC:7400] |
ENSG00000053747 | ENSG00000053747 | laminin subunit alpha 3 [Source:HGNC Symbol;Acc:HGNC:6483] |
ENSG00000136689 | ENSG00000136689 | interleukin 1 receptor antagonist [Source:HGNC Symbol;Acc:HGNC:6000] |
ENSG00000166165 | ENSG00000166165 | creatine kinase B [Source:HGNC Symbol;Acc:HGNC:1991] |
ENSG00000102837 | ENSG00000102837 | olfactomedin 4 [Source:HGNC Symbol;Acc:HGNC:17190] |
ENSG00000112139 | ENSG00000112139 | MAM domain containing glycosylphosphatidylinositol anchor 1 [Source:HGNC Symbol;Acc:HGNC:19267] |
ENSG00000163687 | ENSG00000163687 | deoxyribonuclease 1 like 3 [Source:HGNC Symbol;Acc:HGNC:2959] |
ENSG00000205364 | ENSG00000205364 | metallothionein 1M [Source:HGNC Symbol;Acc:HGNC:14296] |
ENSG00000169429 | ENSG00000169429 | C-X-C motif chemokine ligand 8 [Source:HGNC Symbol;Acc:HGNC:6025] |
ENSG00000255398 | ENSG00000255398 | hydroxycarboxylic acid receptor 3 [Source:HGNC Symbol;Acc:HGNC:16824] |
ENSG00000178726 | ENSG00000178726 | thrombomodulin [Source:HGNC Symbol;Acc:HGNC:11784] |
ENSG00000066056 | ENSG00000066056 | tyrosine kinase with immunoglobulin like and EGF like domains 1 [Source:HGNC Symbol;Acc:HGNC:11809] |
ENSG00000136052 | ENSG00000136052 | solute carrier family 41 member 2 [Source:HGNC Symbol;Acc:HGNC:31045] |
ENSG00000125780 | ENSG00000125780 | transglutaminase 3 [Source:HGNC Symbol;Acc:HGNC:11779] |
ENSG00000178860 | ENSG00000178860 | musculin [Source:HGNC Symbol;Acc:HGNC:7321] |
ENSG00000197632 | ENSG00000197632 | serpin family B member 2 [Source:HGNC Symbol;Acc:HGNC:8584] |
ENSG00000162433 | ENSG00000162433 | adenylate kinase 4 [Source:HGNC Symbol;Acc:HGNC:363] |
test <- create_expt(metadata="sample_sheets/pbmc_samples-manual_switch.xlsx",
file_column="humanfile")
## Reading the sample metadata.
## The sample definitions comprises: 21 rows(samples) and 33 columns(metadata fields).
## Reading count tables.
## Reading count tables with read.table().
## /cbcb/nelsayed-scratch/atb/rnaseq/lpanamensis_2016/preprocessing/hpgl0630/outputs/tophat_hsapiens/accepted_paired.count.xz contains 51046 rows.
## /cbcb/nelsayed-scratch/atb/rnaseq/lpanamensis_2016/preprocessing/hpgl0631/outputs/tophat_hsapiens/accepted_paired.count.xz contains 51046 rows and merges to 51046 rows.
## /cbcb/nelsayed-scratch/atb/rnaseq/lpanamensis_2016/preprocessing/hpgl0632/outputs/tophat_hsapiens/accepted_paired.count.xz contains 51046 rows and merges to 51046 rows.
## /cbcb/nelsayed-scratch/atb/rnaseq/lpanamensis_2016/preprocessing/hpgl0633/outputs/tophat_hsapiens/accepted_paired.count.xz contains 51046 rows and merges to 51046 rows.
## /cbcb/nelsayed-scratch/atb/rnaseq/lpanamensis_2016/preprocessing/hpgl0634/outputs/tophat_hsapiens/accepted_paired.count.xz contains 51046 rows and merges to 51046 rows.
## /cbcb/nelsayed-scratch/atb/rnaseq/lpanamensis_2016/preprocessing/hpgl0635/outputs/tophat_hsapiens/accepted_paired.count.xz contains 51046 rows and merges to 51046 rows.
## /cbcb/nelsayed-scratch/atb/rnaseq/lpanamensis_2016/preprocessing/hpgl0636/outputs/tophat_hsapiens/accepted_paired.count.xz contains 51046 rows and merges to 51046 rows.
## /cbcb/nelsayed-scratch/atb/rnaseq/lpanamensis_2016/preprocessing/hpgl0650/outputs/tophat_hsapiens/accepted_paired.count.xz contains 51046 rows and merges to 51046 rows.
## /cbcb/nelsayed-scratch/atb/rnaseq/lpanamensis_2016/preprocessing/hpgl0651/outputs/tophat_hsapiens/accepted_paired.count.xz contains 51046 rows and merges to 51046 rows.
## /cbcb/nelsayed-scratch/atb/rnaseq/lpanamensis_2016/preprocessing/hpgl0652/outputs/tophat_hsapiens/accepted_paired.count.xz contains 51046 rows and merges to 51046 rows.
## /cbcb/nelsayed-scratch/atb/rnaseq/lpanamensis_2016/preprocessing/hpgl0653/outputs/tophat_hsapiens/accepted_paired.count.xz contains 51046 rows and merges to 51046 rows.
## /cbcb/nelsayed-scratch/atb/rnaseq/lpanamensis_2016/preprocessing/hpgl0654/outputs/tophat_hsapiens/accepted_paired.count.xz contains 51046 rows and merges to 51046 rows.
## /cbcb/nelsayed-scratch/atb/rnaseq/lpanamensis_2016/preprocessing/hpgl0655/outputs/tophat_hsapiens/accepted_paired.count.xz contains 51046 rows and merges to 51046 rows.
## /cbcb/nelsayed-scratch/atb/rnaseq/lpanamensis_2016/preprocessing/hpgl0656/outputs/tophat_hsapiens/accepted_paired.count.xz contains 51046 rows and merges to 51046 rows.
## /cbcb/nelsayed-scratch/atb/rnaseq/lpanamensis_2016/preprocessing/hpgl0657/outputs/tophat_hsapiens/accepted_paired.count.xz contains 51046 rows and merges to 51046 rows.
## /cbcb/nelsayed-scratch/atb/rnaseq/lpanamensis_2016/preprocessing/hpgl0658/outputs/tophat_hsapiens/accepted_paired.count.xz contains 51046 rows and merges to 51046 rows.
## /cbcb/nelsayed-scratch/atb/rnaseq/lpanamensis_2016/preprocessing/hpgl0659/outputs/tophat_hsapiens/accepted_paired.count.xz contains 51046 rows and merges to 51046 rows.
## /cbcb/nelsayed-scratch/atb/rnaseq/lpanamensis_2016/preprocessing/hpgl0660/outputs/tophat_hsapiens/accepted_paired.count.xz contains 51046 rows and merges to 51046 rows.
## /cbcb/nelsayed-scratch/atb/rnaseq/lpanamensis_2016/preprocessing/hpgl0661/outputs/tophat_hsapiens/accepted_paired.count.xz contains 51046 rows and merges to 51046 rows.
## /cbcb/nelsayed-scratch/atb/rnaseq/lpanamensis_2016/preprocessing/hpgl0662/outputs/tophat_hsapiens/accepted_paired.count.xz contains 51046 rows and merges to 51046 rows.
## /cbcb/nelsayed-scratch/atb/rnaseq/lpanamensis_2016/preprocessing/hpgl0663/outputs/tophat_hsapiens/accepted_paired.count.xz contains 51046 rows and merges to 51046 rows.
## Finished reading count tables.
## Matched 51041 annotations and counts.
## Bringing together the count matrix and gene information.
## There were 21, now there are 18 samples.
## This function will replace the expt$expressionset slot with:
## log2(cpm(quant(hpgl(data))))
## It backs up the current data into a slot named:
## expt$backup_expressionset. It will also save copies of each step along the way
## in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
## when invoking limma. The appropriate libsize is the 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: hpgl
## Removing 37484 low-count genes (13557 remaining).
## Step 2: normalizing the data with quant.
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
## Error in plot_pca(test_norm): object 'test_norm' not found
test_batch <- normalize_expt(test, transform="log2", convert="cpm", filter=TRUE,
norm="quant", batch="ruvg")
## This function will replace the expt$expressionset slot with:
## log2(ruvg(cpm(quant(hpgl(data)))))
## It backs up the current data into a slot named:
## expt$backup_expressionset. It will also save copies of each step along the way
## in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
## when invoking limma. The appropriate libsize is the 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
## Step 1: performing count filter with option: hpgl
## Removing 37484 low-count genes (13557 remaining).
## Step 2: normalizing the data with quant.
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
## Error in plot_pca(test_batch): object 'test_batch' not found
## Error in pca_highscores(test_batch, n = 50): object 'test_batch' not found
## Comp.1 Comp.2
## [1,] "31.7:ENSG00000129824" "22.16:ENSG00000100079"
## [2,] "30.31:ENSG00000012817" "16.86:ENSG00000128283"
## [3,] "30:ENSG00000185736" "15.92:ENSG00000170439"
## [4,] "29.65:ENSG00000067048" "13.44:ENSG00000156113"
## [5,] "27.68:ENSG00000114374" "13.39:ENSG00000007350"
## [6,] "25.19:ENSG00000183878" "13.25:ENSG00000123219"
## [7,] "21.54:ENSG00000198692" "13.11:ENSG00000205809"
## [8,] "20.86:ENSG00000067646" "13.01:ENSG00000165810"
## [9,] "18.84:ENSG00000074410" "12.6:ENSG00000150687"
## [10,] "16.59:ENSG00000179344" "12.6:ENSG00000167880"
## [11,] "15.01:ENSG00000134184" "12.36:ENSG00000150510"
## [12,] "14.9:ENSG00000196735" "12.1:ENSG00000167094"
## [13,] "14.46:ENSG00000122641" "12.09:ENSG00000146918"
## [14,] "14.37:ENSG00000121380" "11.77:ENSG00000173372"
## [15,] "14.37:ENSG00000161055" "11.67:ENSG00000214944"
## [16,] "13.88:ENSG00000154620" "11.34:ENSG00000174945"
## [17,] "13.76:ENSG00000166923" "11.01:ENSG00000240247"
## [18,] "13.5:ENSG00000123453" "10.99:ENSG00000205810"
## [19,] "13.06:ENSG00000123454" "10.95:ENSG00000149635"
## [20,] "13.02:ENSG00000108688" "10.72:ENSG00000102837"
## [21,] "12.49:ENSG00000131435" "10.69:ENSG00000109956"
## [22,] "12.47:ENSG00000136449" "10.68:ENSG00000061337"
## [23,] "12.22:ENSG00000110092" "10.6:ENSG00000206047"
## [24,] "12.05:ENSG00000120645" "10.51:ENSG00000205846"
## [25,] "11.83:ENSG00000205057" "10.5:ENSG00000183542"
## [26,] "11.54:ENSG00000123689" "10.04:ENSG00000113721"
## [27,] "11.34:ENSG00000170160" "9.932:ENSG00000101188"
## [28,] "11.32:ENSG00000113302" "9.836:ENSG00000113657"
## [29,] "11.2:ENSG00000102970" "9.783:ENSG00000196421"
## [30,] "11.08:ENSG00000198203" "9.572:ENSG00000145708"
## [31,] "11.08:ENSG00000184371" "9.241:ENSG00000133962"
## [32,] "11.07:ENSG00000149635" "9.213:ENSG00000148488"
## [33,] "10.94:ENSG00000204644" "8.953:ENSG00000125780"
## [34,] "10.83:ENSG00000050165" "8.802:ENSG00000243772"
## [35,] "10.66:ENSG00000150510" "8.796:ENSG00000170160"
## [36,] "10.47:ENSG00000182871" "8.76:ENSG00000162383"
## [37,] "10.42:ENSG00000134668" "8.719:ENSG00000108370"
## [38,] "10.23:ENSG00000164530" "8.705:ENSG00000054793"
## [39,] "10.22:ENSG00000178860" "8.564:ENSG00000168329"
## [40,] "10.19:ENSG00000109471" "8.503:ENSG00000184349"
## [41,] "10.17:ENSG00000197121" "8.377:ENSG00000100450"
## [42,] "10.15:ENSG00000164400" "8.224:ENSG00000173239"
## [43,] "9.967:ENSG00000180616" "8.088:ENSG00000086548"
## [44,] "9.874:ENSG00000069482" "8.001:ENSG00000111199"
## [45,] "9.8:ENSG00000133687" "7.934:ENSG00000239839"
## [46,] "9.786:ENSG00000134716" "7.928:ENSG00000139567"
## [47,] "9.778:ENSG00000115009" "7.858:ENSG00000188820"
## [48,] "9.737:ENSG00000197822" "7.744:ENSG00000101425"
## [49,] "9.705:ENSG00000172724" "7.624:ENSG00000176083"
## [50,] "9.686:ENSG00000181634" "7.616:ENSG00000169908"
annot <- fData(hs_inf_filt)
comp1_genes <- unlist(strsplit(x=as.character(scores$highest[, 1]), split=":"))
idx <- grep(pattern="^ENSG", x=comp1_genes)
comp1_genes <- comp1_genes[idx]
comp1_high_scores <- exclude_genes_expt(test_batch, ids=comp1_genes, method="keep")
## Error in exclude_genes_expt(test_batch, ids = comp1_genes, method = "keep"): object 'test_batch' not found
plot_sample_heatmap(
data=sm(normalize_expt(comp1_high_scores, transform="log2", convert="cpm", norm="quant")),
row_label=NULL)
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
ensembl_gene_id | description | |
---|---|---|
ENSG00000129824 | ENSG00000129824 | ribosomal protein S4, Y-linked 1 [Source:HGNC Symbol;Acc:HGNC:10425] |
ENSG00000012817 | ENSG00000012817 | lysine demethylase 5D [Source:HGNC Symbol;Acc:HGNC:11115] |
ENSG00000185736 | ENSG00000185736 | adenosine deaminase, RNA specific B2 (inactive) [Source:HGNC Symbol;Acc:HGNC:227] |
ENSG00000067048 | ENSG00000067048 | DEAD-box helicase 3, Y-linked [Source:HGNC Symbol;Acc:HGNC:2699] |
ENSG00000114374 | ENSG00000114374 | ubiquitin specific peptidase 9, Y-linked [Source:HGNC Symbol;Acc:HGNC:12633] |
ENSG00000183878 | ENSG00000183878 | ubiquitously transcribed tetratricopeptide repeat containing, Y-linked [Source:HGNC Symbol;Acc:HGNC:12638] |
ENSG00000198692 | ENSG00000198692 | eukaryotic translation initiation factor 1A, Y-linked [Source:HGNC Symbol;Acc:HGNC:3252] |
ENSG00000067646 | ENSG00000067646 | zinc finger protein, Y-linked [Source:HGNC Symbol;Acc:HGNC:12870] |
ENSG00000074410 | ENSG00000074410 | carbonic anhydrase 12 [Source:HGNC Symbol;Acc:HGNC:1371] |
ENSG00000179344 | ENSG00000179344 | major histocompatibility complex, class II, DQ beta 1 [Source:HGNC Symbol;Acc:HGNC:4944] |
ENSG00000134184 | ENSG00000134184 | glutathione S-transferase mu 1 [Source:HGNC Symbol;Acc:HGNC:4632] |
ENSG00000196735 | ENSG00000196735 | major histocompatibility complex, class II, DQ alpha 1 [Source:HGNC Symbol;Acc:HGNC:4942] |
ENSG00000122641 | ENSG00000122641 | inhibin beta A subunit [Source:HGNC Symbol;Acc:HGNC:6066] |
ENSG00000121380 | ENSG00000121380 | BCL2 like 14 [Source:HGNC Symbol;Acc:HGNC:16657] |
ENSG00000161055 | ENSG00000161055 | secretoglobin family 3A member 1 [Source:HGNC Symbol;Acc:HGNC:18384] |
ENSG00000154620 | ENSG00000154620 | thymosin beta 4, Y-linked [Source:HGNC Symbol;Acc:HGNC:11882] |
ENSG00000166923 | ENSG00000166923 | gremlin 1, DAN family BMP antagonist [Source:HGNC Symbol;Acc:HGNC:2001] |
ENSG00000123453 | ENSG00000123453 | sarcosine dehydrogenase [Source:HGNC Symbol;Acc:HGNC:10536] |
ENSG00000123454 | ENSG00000123454 | dopamine beta-hydroxylase [Source:HGNC Symbol;Acc:HGNC:2689] |
ENSG00000108688 | ENSG00000108688 | C-C motif chemokine ligand 7 [Source:HGNC Symbol;Acc:HGNC:10634] |
ENSG00000131435 | ENSG00000131435 | PDZ and LIM domain 4 [Source:HGNC Symbol;Acc:HGNC:16501] |
ENSG00000136449 | ENSG00000136449 | MYCBP associated protein [Source:HGNC Symbol;Acc:HGNC:19677] |
ENSG00000110092 | ENSG00000110092 | cyclin D1 [Source:HGNC Symbol;Acc:HGNC:1582] |
ENSG00000120645 | ENSG00000120645 | IQ motif and Sec7 domain 3 [Source:HGNC Symbol;Acc:HGNC:29193] |
ENSG00000205057 | ENSG00000205057 | chronic lymphocytic leukemia up-regulated 1 opposite strand [Source:HGNC Symbol;Acc:HGNC:24070] |
ENSG00000123689 | ENSG00000123689 | G0/G1 switch 2 [Source:HGNC Symbol;Acc:HGNC:30229] |
ENSG00000170160 | ENSG00000170160 | coiled-coil domain containing 144A [Source:HGNC Symbol;Acc:HGNC:29072] |
ENSG00000113302 | ENSG00000113302 | interleukin 12B [Source:HGNC Symbol;Acc:HGNC:5970] |
ENSG00000102970 | ENSG00000102970 | C-C motif chemokine ligand 17 [Source:HGNC Symbol;Acc:HGNC:10615] |
ENSG00000198203 | ENSG00000198203 | sulfotransferase family 1C member 2 [Source:HGNC Symbol;Acc:HGNC:11456] |
ENSG00000184371 | ENSG00000184371 | colony stimulating factor 1 [Source:HGNC Symbol;Acc:HGNC:2432] |
ENSG00000149635 | ENSG00000149635 | osteoclast stimulatory transmembrane protein [Source:HGNC Symbol;Acc:HGNC:16116] |
ENSG00000204644 | ENSG00000204644 | ZFP57 zinc finger protein [Source:HGNC Symbol;Acc:HGNC:18791] |
ENSG00000050165 | ENSG00000050165 | dickkopf WNT signaling pathway inhibitor 3 [Source:HGNC Symbol;Acc:HGNC:2893] |
ENSG00000150510 | ENSG00000150510 | family with sequence similarity 124 member A [Source:HGNC Symbol;Acc:HGNC:26413] |
ENSG00000182871 | ENSG00000182871 | collagen type XVIII alpha 1 chain [Source:HGNC Symbol;Acc:HGNC:2195] |
ENSG00000134668 | ENSG00000134668 | SPOC domain containing 1 [Source:HGNC Symbol;Acc:HGNC:26338] |
ENSG00000164530 | ENSG00000164530 | peptidase inhibitor 16 [Source:HGNC Symbol;Acc:HGNC:21245] |
ENSG00000178860 | ENSG00000178860 | musculin [Source:HGNC Symbol;Acc:HGNC:7321] |
ENSG00000109471 | ENSG00000109471 | interleukin 2 [Source:HGNC Symbol;Acc:HGNC:6001] |
ENSG00000197121 | ENSG00000197121 | post-GPI attachment to proteins 1 [Source:HGNC Symbol;Acc:HGNC:25712] |
ENSG00000164400 | ENSG00000164400 | colony stimulating factor 2 [Source:HGNC Symbol;Acc:HGNC:2434] |
ENSG00000180616 | ENSG00000180616 | somatostatin receptor 2 [Source:HGNC Symbol;Acc:HGNC:11331] |
ENSG00000069482 | ENSG00000069482 | galanin and GMAP prepropeptide [Source:HGNC Symbol;Acc:HGNC:4114] |
ENSG00000133687 | ENSG00000133687 | transmembrane and tetratricopeptide repeat containing 1 [Source:HGNC Symbol;Acc:HGNC:24099] |
ENSG00000134716 | ENSG00000134716 | cytochrome P450 family 2 subfamily J member 2 [Source:HGNC Symbol;Acc:HGNC:2634] |
ENSG00000115009 | ENSG00000115009 | C-C motif chemokine ligand 20 [Source:HGNC Symbol;Acc:HGNC:10619] |
ENSG00000197822 | ENSG00000197822 | occludin [Source:HGNC Symbol;Acc:HGNC:8104] |
ENSG00000172724 | ENSG00000172724 | C-C motif chemokine ligand 19 [Source:HGNC Symbol;Acc:HGNC:10617] |
ENSG00000181634 | ENSG00000181634 | tumor necrosis factor superfamily member 15 [Source:HGNC Symbol;Acc:HGNC:11931] |
## This function will replace the expt$expressionset slot with:
## hpgl(data)
## It backs up the current data into a slot named:
## expt$backup_expressionset. It will also save copies of each step along the way
## in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
## when invoking limma. The appropriate libsize is the non-log(cpm(normalized)).
## This is most likely kept at:
## 'new_expt$normalized$intermediate_counts$normalization$libsizes'
## A copy of this may also be found at:
## new_expt$best_libsize
## Leaving the data in its current base format, keep in mind that
## some metrics are easier to see when the data is log2 transformed, but
## EdgeR/DESeq do not accept transformed data.
## Leaving the data unconverted. It is often advisable to cpm/rpkm
## the data to normalize for sampling differences, keep in mind though that rpkm
## has some annoying biases, and voom() by default does a cpm (though hpgl_voom()
## will try to detect this).
## Leaving the data unnormalized. This is necessary for DESeq, but
## EdgeR/limma might benefit from normalization. Good choices include quantile,
## size-factor, tmm, etc.
## 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: hpgl
## Removing 37484 low-count genes (13557 remaining).
## Step 2: not normalizing the data.
## Step 3: not converting the data.
## Step 4: not transforming the data.
## Step 5: not doing batch correction.
excel_file <- glue::glue("excel/{rundate}_hs_infect_switchone_contr-v{ver}.xlsx")
test_de <- all_pairwise(test_filt, model_batch=TRUE, parallel=FALSE)
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
## Error in extract_de_plots(test_de): object 'test_de' not found
## Error in combine_de_tables(test_de, excel = excel_file): object 'test_de' not found
excel_file <- glue::glue("excel/{rundate}_hs_infect_switchone_sig-v{ver}.xlsx")
test_sig <- extract_significant_genes(test_table, according_to="deseq", p=0.05, excel=excel_file)
## Error in extract_significant_genes(test_table, according_to = "deseq", : object 'test_table' not found
## Error in head(test_sig$deseq$ups[[1]]): object 'test_sig' not found
## Error in extract_de_plots(test_de, p = 0.1): object 'test_de' not found
## Error in eval(expr, envir, enclos): object 'test_ma' not found
## Error in plot_histogram(test_table[["data"]][["sh_vs_chr"]][, c("deseq_p")]): object 'test_table' not found
## Error in eval(expr, envir, enclos): object 'test_table' not found
In the following block, I will take the comparisons performed without any batch in the model/adjustment and use them to search for shared/unique genes among the self-healing vs. uninfected and the chronic vs. uninfected.
hs_pairwise_nobatch_pca <- sm(plot_pca(hs_inf_filt, convert="cpm", transform="log2"))
hs_pairwise_nobatch_pca$plot
hs_pairwise_nobatch <- sm(all_pairwise(hs_uninf_filt, model_batch=FALSE, parallel=FALSE,
do_ebseq=FALSE))
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
excel_file <- glue::glue("excel/{rundate}_hs_infect_nobatch_contr-v{ver}.xlsx")
hs_combined_nobatch <- sm(combine_de_tables(
hs_pairwise_nobatch,
excel=excel_file,
keepers=keepers))
## Error in combine_de_tables(hs_pairwise_nobatch, excel = excel_file, keepers = keepers): object 'hs_pairwise_nobatch' not found
excel_file <- glue::glue("excel/{rundate}_hs_infect_nobatch_sig-v{ver}.xlsx")
hs_sig_nobatch <- sm(extract_significant_genes(
hs_combined_nobatch,
excel=excel_file))
## Error in extract_significant_genes(hs_combined_nobatch, excel = excel_file): object 'hs_combined_nobatch' not found
## Error in eval(expr, envir, enclos): object 'hs_sig_nobatch' not found
limma_deseq_order <- rank_order_scatter(
hs_pairwise_nobatch,
first_type="limma", second_type="deseq",
first_table=1, second_table=1)
## Error in rank_order_scatter(hs_pairwise_nobatch, first_type = "limma", : object 'hs_pairwise_nobatch' not found
## Error in eval(expr, envir, enclos): object 'limma_deseq_order' not found
up_lst <- list(
"sh_up" = rownames(hs_sig_nobatch[["deseq"]][["ups"]][["sh_nil"]]),
"ch_up" = rownames(hs_sig_nobatch[["deseq"]][["ups"]][["ch_nil"]]))
## Error in rownames(hs_sig_nobatch[["deseq"]][["ups"]][["sh_nil"]]): object 'hs_sig_nobatch' not found
## Error in Vennerable::Venn(Sets = up_lst): object 'up_lst' not found
## Error in Vennerable::plot(nobatch_up_venn, doWeights = FALSE): object 'nobatch_up_venn' not found
## Maria-Adelaida and Najib asked about comparing the following pair of data against
## The intersection of a bunch of stuff later... So don't forget this!
## This gives me the genes only up in self vs. nil
nobatch_up_sh_solo <- nobatch_up_venn@IntersectionSets[["10"]]
## Error in eval(expr, envir, enclos): object 'nobatch_up_venn' not found
## Error in eval(expr, envir, enclos): object 'nobatch_up_venn' not found
down_lst <- list(
"sh_down" = rownames(hs_sig_nobatch[["deseq"]][["downs"]][["sh_nil"]]),
"ch_down" = rownames(hs_sig_nobatch[["deseq"]][["downs"]][["ch_nil"]]))
## Error in rownames(hs_sig_nobatch[["deseq"]][["downs"]][["sh_nil"]]): object 'hs_sig_nobatch' not found
## Error in Vennerable::Venn(Sets = down_lst): object 'down_lst' not found
## Error in Vennerable::plot(nobatch_down_venn, doWeights = FALSE): object 'nobatch_down_venn' not found
Repeat the previous set of analyses with d107/108/110 in the model.
hs_pairwise_batch_pca <- sm(plot_pca(hs_inf_filt, batch="limma", convert="cpm", transform="log2"))
hs_pairwise_batch_pca$plot
hs_pairwise_batch <- sm(all_pairwise(hs_uninf_filt, parallel=FALSE,
model_batch=TRUE, do_ebseq=FALSE))
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
excel_file <- glue::glue("excel/{rundate}_hs_infect_patbatch_contr-v{ver}.xlsx")
hs_combined_batch <- sm(combine_de_tables(
hs_pairwise_batch,
excel=excel_file,
keepers=keepers))
## Error in combine_de_tables(hs_pairwise_batch, excel = excel_file, keepers = keepers): object 'hs_pairwise_batch' not found
excel_file <- glue::glue("excel/{rundate}_hs_infect_patbatch_sig-v{ver}.xlsx")
hs_sig_batch <- sm(extract_significant_genes(
hs_combined_batch,
excel=excel_file))
## Error in extract_significant_genes(hs_combined_batch, excel = excel_file): object 'hs_combined_batch' not found
## Error in eval(expr, envir, enclos): object 'hs_sig_batch' not found
up_lst <- list(
"sh_up" = rownames(hs_sig_batch[["deseq"]][["ups"]][["sh_nil"]]),
"ch_up" = rownames(hs_sig_batch[["deseq"]][["ups"]][["ch_nil"]]))
## Error in rownames(hs_sig_batch[["deseq"]][["ups"]][["sh_nil"]]): object 'hs_sig_batch' not found
## Error in Vennerable::Venn(Sets = up_lst): object 'up_lst' not found
## Error in summary(batch_up_venn@IntersectionSets): object 'batch_up_venn' not found
## Error in Vennerable::plot(batch_up_venn, doWeights = FALSE): object 'batch_up_venn' not found
down_lst <- list(
"sh_down" = rownames(hs_sig_batch[["deseq"]][["downs"]][["sh_nil"]]),
"ch_down" = rownames(hs_sig_batch[["deseq"]][["downs"]][["ch_nil"]]))
## Error in rownames(hs_sig_batch[["deseq"]][["downs"]][["sh_nil"]]): object 'hs_sig_batch' not found
## Error in Vennerable::Venn(Sets = down_lst): object 'down_lst' not found
## Error in Vennerable::plot(batch_down_venn, doWeights = FALSE): object 'batch_down_venn' not found
## Error in summary(batch_down_venn@IntersectionSets): object 'batch_down_venn' not found
## Testing method: limma.
## Error in compare_de_results(hs_combined_nobatch, hs_combined_batch, cor_method = "spearman"): object 'hs_combined_nobatch' not found
## Error in eval(expr, envir, enclos): object 'similar' not found
The following block writes out the unique/shared genes observed among the contrasts which included donor in the model.
kept_columns <- c("ensembltranscriptid", "ensemblgeneid", "description",
"deseq_logfc", "deseq_adjp")
start_data <- hs_sig_batch[["deseq"]]
## Error in eval(expr, envir, enclos): object 'hs_sig_batch' not found
## Error in eval(expr, envir, enclos): object 'batch_up_venn' not found
## Error in eval(expr, envir, enclos): object 'batch_up_venn' not found
## Error in eval(expr, envir, enclos): object 'batch_up_venn' not found
## Error in eval(expr, envir, enclos): object 'batch_down_venn' not found
## Error in eval(expr, envir, enclos): object 'batch_down_venn' not found
## Error in eval(expr, envir, enclos): object 'batch_down_venn' not found
xls_result <- write_xls(
data=start_data[["ups"]][["sh_nil"]][sh_up_solo_genes, kept_columns],
sheet="sh_up_solo")
## Error in write_xls(data = start_data[["ups"]][["sh_nil"]][sh_up_solo_genes, : object 'start_data' not found
xls_result <- write_xls(
data=start_data[["ups"]][["ch_nil"]][ch_up_solo_genes, kept_columns],
sheet="ch_up_solo", wb=xls_result[["workbook"]])
## Error in write_xls(data = start_data[["ups"]][["ch_nil"]][ch_up_solo_genes, : object 'start_data' not found
xls_result <- write_xls(
data=start_data[["ups"]][["ch_nil"]][shch_up_shared_genes, kept_columns],
sheet="shch_up_shared", wb=xls_result[["workbook"]])
## Error in write_xls(data = start_data[["ups"]][["ch_nil"]][shch_up_shared_genes, : object 'start_data' not found
xls_result <- write_xls(
data=start_data[["downs"]][["sh_nil"]][sh_down_solo_genes, kept_columns],
sheet="sh_down_solo")
## Error in write_xls(data = start_data[["downs"]][["sh_nil"]][sh_down_solo_genes, : object 'start_data' not found
xls_result <- write_xls(
data=start_data[["downs"]][["ch_nil"]][ch_down_solo_genes, kept_columns],
sheet="ch_down_solo", wb=xls_result[["workbook"]])
## Error in write_xls(data = start_data[["downs"]][["ch_nil"]][ch_down_solo_genes, : object 'start_data' not found
xls_result <- write_xls(
data=start_data[["downs"]][["ch_nil"]][shch_down_shared_genes, kept_columns],
sheet="shch_down_shared", wb=xls_result[["workbook"]],
excel=paste0("excel/figure_5a_stuff-v", ver, ".xlsx"))
## Error in write_xls(data = start_data[["downs"]][["ch_nil"]][shch_down_shared_genes, : object 'start_data' not found
Repeat, this time attmepting to tamp down the variance by person.
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
## Error in eval(expr, envir, enclos): object 'hs_pairwise_ssva_pca' not found
hs_pairwise_ssva <- sm(all_pairwise(hs_uninf_filt, model_batch="ssva",
parallel=FALSE, do_ebseq=FALSE))
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
excel_file <- glue::glue("excel/{rundate}_hs_infect_fsva_contr-v{ver}.xlsx")
hs_combined_ssva <- sm(combine_de_tables(
hs_pairwise_ssva,
excel=excel_file,
keepers=keepers))
## Error in combine_de_tables(hs_pairwise_ssva, excel = excel_file, keepers = keepers): object 'hs_pairwise_ssva' not found
excel_file <- glue::glue("excel/{rundate}_hs_infect_fsva_sig-v{ver}.xlsx")
hs_sig_ssva <- sm(extract_significant_genes(
hs_combined_ssva,
excel=excel_file))
## Error in extract_significant_genes(hs_combined_ssva, excel = excel_file): object 'hs_combined_ssva' not found
## Error in eval(expr, envir, enclos): object 'hs_sig_ssva' not found
up_lst <- list(
"sh_up" = rownames(hs_sig_ssva[["deseq"]][["ups"]][["sh_nil"]]),
"ch_up" = rownames(hs_sig_ssva[["deseq"]][["ups"]][["ch_nil"]]))
## Error in rownames(hs_sig_ssva[["deseq"]][["ups"]][["sh_nil"]]): object 'hs_sig_ssva' not found
## Error in Vennerable::Venn(Sets = up_lst): object 'up_lst' not found
## Error in summary(ssva_up_venn@IntersectionSets): object 'ssva_up_venn' not found
## Error in Vennerable::plot(ssva_up_venn, doWeights = FALSE): object 'ssva_up_venn' not found
down_lst <- list(
"sh_down" = rownames(hs_sig_ssva[["deseq"]][["downs"]][["sh_nil"]]),
"ch_down" = rownames(hs_sig_ssva[["deseq"]][["downs"]][["ch_nil"]]))
## Error in rownames(hs_sig_ssva[["deseq"]][["downs"]][["sh_nil"]]): object 'hs_sig_ssva' not found
## Error in Vennerable::Venn(Sets = down_lst): object 'down_lst' not found
## Error in Vennerable::plot(ssva_down_venn, doWeights = FALSE): object 'ssva_down_venn' not found
## Error in summary(ssva_down_venn@IntersectionSets): object 'ssva_down_venn' not found
## Error in compare_de_results(hs_combined_nobatch, hs_combined_ssva, cor_method = "spearman"): object 'hs_combined_nobatch' not found
## Error in eval(expr, envir, enclos): object 'similar' not found
## Error in compare_de_results(hs_combined_batch, hs_combined_ssva, cor_method = "spearman"): object 'hs_combined_batch' not found
## Error in eval(expr, envir, enclos): object 'similar' not found
Repeat again using fsva.
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
## Error in eval(expr, envir, enclos): object 'hs_pairwise_fsva_pca' not found
hs_pairwise_fsva <- sm(all_pairwise(hs_uninf_filt, model_batch="fsva",
parallel=FALSE, do_ebseq=FALSE))
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
excel_file <- glue::glue("excel/{rundate}_hs_infect_fsva_contr-v{ver}.xlsx")
hs_combined_fsva <- sm(combine_de_tables(
hs_pairwise_fsva,
excel=excel_file,
keepers=keepers))
## Error in combine_de_tables(hs_pairwise_fsva, excel = excel_file, keepers = keepers): object 'hs_pairwise_fsva' not found
excel_file <- glue::glue("excel/{rundate}_hs_infect_fsva_sig-v{ver}.xlsx")
hs_sig_fsva <- sm(extract_significant_genes(
hs_combined_fsva,
excel=excel_file))
## Error in extract_significant_genes(hs_combined_fsva, excel = excel_file): object 'hs_combined_fsva' not found
up_lst <- list(
"sh_up" = rownames(hs_sig_fsva[["deseq"]][["ups"]][["sh_nil"]]),
"ch_up" = rownames(hs_sig_fsva[["deseq"]][["ups"]][["ch_nil"]]))
## Error in rownames(hs_sig_fsva[["deseq"]][["ups"]][["sh_nil"]]): object 'hs_sig_fsva' not found
## Error in Vennerable::Venn(Sets = up_lst): object 'up_lst' not found
## Error in summary(fsva_up_venn@IntersectionSets): object 'fsva_up_venn' not found
## Error in Vennerable::plot(fsva_up_venn, doWeights = FALSE): object 'fsva_up_venn' not found
down_lst <- list(
"sh_down" = rownames(hs_sig_fsva[["deseq"]][["downs"]][["sh_nil"]]),
"ch_down" = rownames(hs_sig_fsva[["deseq"]][["downs"]][["ch_nil"]]))
## Error in rownames(hs_sig_fsva[["deseq"]][["downs"]][["sh_nil"]]): object 'hs_sig_fsva' not found
## Error in Vennerable::Venn(Sets = down_lst): object 'down_lst' not found
## Error in Vennerable::plot(fsva_down_venn, doWeights = FALSE): object 'fsva_down_venn' not found
## Error in summary(fsva_down_venn@IntersectionSets): object 'fsva_down_venn' not found
##hs_sig_fsva$deseq$counts$common_solos_fsva <- sm(subset_significants(hs_sig_fsva))
similar <- sm(compare_de_results(hs_combined_nobatch, hs_combined_fsva,
cor_method="spearman"))
## Error in compare_de_results(hs_combined_nobatch, hs_combined_fsva, cor_method = "spearman"): object 'hs_combined_nobatch' not found
## Error in eval(expr, envir, enclos): object 'similar' not found
## Error in compare_de_results(hs_combined_ssva, hs_combined_fsva, cor_method = "spearman"): object 'hs_combined_ssva' not found
## Error in eval(expr, envir, enclos): object 'similar' not found
Repeat once again, this time using combat to try to limit the contribution of the strain to the data. I do not think we will ever use this set of contrasts, so I will deactivate it but leave it here if it is required later.
old_condition <- hs_uninf$design$condition
names(old_condition) <- hs_uninf$design$sampleid
new_condition <- paste0(hs_uninf$design$state, '_', hs_uninf$design$donor)
hs_uninf_recond <- set_expt_factors(hs_uninf_filt,
batch="pathogenstrain",
condition=new_condition)
combat_input <- normalize_expt(hs_uninf_recond, batch="combat_scale", filter=TRUE)
## This function will replace the expt$expressionset slot with:
## combat_scale(hpgl(data))
## It backs up the current data into a slot named:
## expt$backup_expressionset. It will also save copies of each step along the way
## in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
## when invoking limma. The appropriate libsize is the non-log(cpm(normalized)).
## This is most likely kept at:
## 'new_expt$normalized$intermediate_counts$normalization$libsizes'
## A copy of this may also be found at:
## new_expt$best_libsize
## Leaving the data in its current base format, keep in mind that
## some metrics are easier to see when the data is log2 transformed, but
## EdgeR/DESeq do not accept transformed data.
## Leaving the data unconverted. It is often advisable to cpm/rpkm
## the data to normalize for sampling differences, keep in mind though that rpkm
## has some annoying biases, and voom() by default does a cpm (though hpgl_voom()
## will try to detect this).
## Leaving the data unnormalized. This is necessary for DESeq, but
## EdgeR/limma might benefit from normalization. Good choices include quantile,
## size-factor, tmm, etc.
## Step 1: performing count filter with option: hpgl
## Removing 0 low-count genes (12086 remaining).
## Step 2: not normalizing the data.
## Step 3: not converting the data.
## Step 4: not transforming the data.
## Step 5: doing batch correction with combat_scale.
## 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.
## batch_counts: Before batch correction, 320 entries are >= 0.
## Passing off to all_adjusters.
## batch_counts: Before batch/surrogate estimation, 320 entries are x<=0.
## The be method chose 2 surrogate variable(s).
## batch_counts: Using combat with a prior and with scaling.
## The number of elements which are < 0 after batch correction is: 6972
## The variable low_to_zero sets whether to change <0 values to 0 and is: FALSE
combat_input <- set_expt_conditions(combat_input, fact=old_condition)
combat_input <- set_expt_colors(combat_input, colors=c("green","blue","red"))
## The new colors are a character, changing according to condition.
hs_pairwise_combatpath <- sm(all_pairwise(combat_input, model_batch=FALSE, parallel=FALSE,
force=TRUE))
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
excel_file <- glue::glue("excel/{rundate}_hs_infect_combatpath_contr-v{ver}.xlsx")
hs_combined_combatpath <- sm(combine_de_tables(
hs_pairwise_combatpath,
excel=excel_file,
keepers=keepers))
## Error in combine_de_tables(hs_pairwise_combatpath, excel = excel_file, : object 'hs_pairwise_combatpath' not found
excel_file <- glue::glue("excel/{rundate}_hs_infect_combatpath_sig-v{ver}.xlsx")
hs_sig_combatpath <- sm(extract_significant_genes(
hs_combined_combatpath,
excel=excel_file))
## Error in extract_significant_genes(hs_combined_combatpath, excel = excel_file): object 'hs_combined_combatpath' not found
## Error in eval(expr, envir, enclos): object 'hs_sig_combatpath' not found
up_lst <- list(
"sh_up" = rownames(hs_sig_combatpath[["deseq"]][["ups"]][["sh_nil"]]),
"ch_up" = rownames(hs_sig_combatpath[["deseq"]][["ups"]][["ch_nil"]]))
## Error in rownames(hs_sig_combatpath[["deseq"]][["ups"]][["sh_nil"]]): object 'hs_sig_combatpath' not found
## Error in Vennerable::Venn(Sets = up_lst): object 'up_lst' not found
## Error in summary(combatpath_up_venn@IntersectionSets): object 'combatpath_up_venn' not found
## Error in Vennerable::plot(combatpath_up_venn, doWeights = FALSE): object 'combatpath_up_venn' not found
down_lst <- list(
"sh_down" = rownames(hs_sig_combatpath[["deseq"]][["downs"]][["sh_nil"]]),
"ch_down" = rownames(hs_sig_combatpath[["deseq"]][["downs"]][["ch_nil"]]))
## Error in rownames(hs_sig_combatpath[["deseq"]][["downs"]][["sh_nil"]]): object 'hs_sig_combatpath' not found
## Error in Vennerable::Venn(Sets = down_lst): object 'down_lst' not found
## Error in Vennerable::plot(combatpath_down_venn, doWeights = FALSE): object 'combatpath_down_venn' not found
## Error in summary(combatpath_down_venn@IntersectionSets): object 'combatpath_down_venn' not found
similar <- sm(compare_de_results(hs_combined_nobatch, hs_combined_combatpath,
cor_method="spearman"))
## Error in compare_de_results(hs_combined_nobatch, hs_combined_combatpath, : object 'hs_combined_nobatch' not found
## Error in eval(expr, envir, enclos): object 'similar' not found
## Error in compare_de_results(hs_combined_fsva, hs_combined_combatpath, : object 'hs_combined_fsva' not found
## Error in eval(expr, envir, enclos): object 'similar' not found
For each of the following, perform a simple DE and see what happens:
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.
## chr_5430_d108 chr_5397_d108 chr_2504_d108 sh_2272_d108 sh_1022_d108
## chr chr chr sh sh
## sh_2189_d108 chr_5430_d110 chr_5397_d110 chr_2504_d110 sh_2272_d110
## sh chr chr chr sh
## sh_1022_d110 sh_2189_d110 chr_5430_d107 chr_5397_d107 chr_2504_d107
## sh sh chr chr chr
## sh_2272_d107 sh_1022_d107 sh_2189_d107
## sh sh sh
## Levels: sh chr
## 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, parallel=FALSE,
model_batch=FALSE, do_ebseq=FALSE))
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
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")))
## Error in combine_de_tables(hs_uninf_filtv2_pairwise, keepers = uninf_patient_keepers, : object 'hs_uninf_filtv2_pairwise' not found
hs_uninf_filtv2_sig <- sm(extract_significant_genes(
hs_uninf_filtv2_combined, according_to="deseq",
excel=paste0("excel/", rundate, "hs_infect_patient_nobatch_sig-v", ver, ".xlsx")))
## Error in extract_significant_genes(hs_uninf_filtv2_combined, according_to = "deseq", : object 'hs_uninf_filtv2_combined' not found
##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")))
Now we want to look at intersections from the perspective of contrasts performed comparing the self-healing/chronic vs. uninfected for the three donors separately.
This time for each of the three donors: self-healing up vs. uninfected.
## Error in eval(expr, envir, enclos): object 'hs_uninf_filtv2_sig' not found
## Do this in a way which is not stupid...
comp_lst <- list(
"donor_107" = rownames(up_sig[["d107_shun"]]),
"donor_108" = rownames(up_sig[["d108_shun"]]),
"donor_110" = rownames(up_sig[["d110_shun"]]))
## Error in rownames(up_sig[["d107_shun"]]): object 'up_sig' not found
## Error in Vennerable::Venn(Sets = comp_lst): object 'comp_lst' not found
## Error in Vennerable::plot(comp_shun_up, doWeights = FALSE): object 'comp_shun_up' not found
## Error in eval(expr, envir, enclos): object 'comp_shun_up' not found
## Error in eval(expr, envir, enclos): object 'hs_uninf_filtv2_combined' not found
## Error in eval(expr, envir, enclos): object 'hs_uninf_filtv2_combined' not found
## Error in eval(expr, envir, enclos): object 'hs_uninf_filtv2_combined' not found
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")
## Error in merge(de_table_shared_up_sh_first[, c("description", "deseq_logfc", : object 'de_table_shared_up_sh_first' not found
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")
## Error in merge(de_table_shared_up_sh_all, de_table_shared_up_sh_third[, : object 'de_table_shared_up_sh_all' not found
## Error in eval(expr, envir, enclos): object 'de_table_shared_up_sh_all' not found
## Error in eval(expr, envir, enclos): object 'de_table_shared_up_sh_all' not found
colnames(de_table_shared_up_sh_all) <- c("description", "logfc_107", "adjp_107",
"logfc_108", "adjp_108", "logfc_110", "adjp_110")
## Error in colnames(de_table_shared_up_sh_all) <- c("description", "logfc_107", : object 'de_table_shared_up_sh_all' not found
write.csv(de_table_shared_up_sh_all,
file=paste0("images/", rundate, "_de_table_shared_up_sh_all.csv"))
## Error in is.data.frame(x): object 'de_table_shared_up_sh_all' not found
This time for each of the three donors: chronic up vs. uninfected.
## Error in eval(expr, envir, enclos): object 'hs_uninf_filtv2_sig' not found
comp_lst <- list(
"donor_107" = rownames(up_sig[["d107_chun"]]),
"donor_108" = rownames(up_sig[["d108_chun"]]),
"donor_110" = rownames(up_sig[["d110_chun"]]))
## Error in rownames(up_sig[["d107_chun"]]): object 'up_sig' not found
## Error in Vennerable::Venn(Sets = comp_lst): object 'comp_lst' not found
## Error in Vennerable::plot(comp_chun_up, doWeights = FALSE): object 'comp_chun_up' not found
## Error in eval(expr, envir, enclos): object 'comp_chun_up' not found
## Error in eval(expr, envir, enclos): object 'hs_uninf_filtv2_combined' not found
## Error in eval(expr, envir, enclos): object 'hs_uninf_filtv2_combined' not found
## Error in eval(expr, envir, enclos): object 'hs_uninf_filtv2_combined' not found
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")
## Error in merge(de_table_shared_up_ch_first[, c("description", "deseq_logfc", : object 'de_table_shared_up_ch_first' not found
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")
## Error in merge(de_table_shared_up_ch_all, de_table_shared_up_ch_third[, : object 'de_table_shared_up_ch_all' not found
## Error in eval(expr, envir, enclos): object 'de_table_shared_up_ch_all' not found
## Error in eval(expr, envir, enclos): object 'de_table_shared_up_ch_all' not found
colnames(de_table_shared_up_ch_all) <- c("description", "logfc_107", "adjp_107",
"logfc_108", "adjp_108", "logfc_110", "adjp_110")
## Error in colnames(de_table_shared_up_ch_all) <- c("description", "logfc_107", : object 'de_table_shared_up_ch_all' not found
write.csv(de_table_shared_up_ch_all,
file=paste0("images/", rundate, "_de_table_shared_up_ch_all.csv"))
## Error in is.data.frame(x): object 'de_table_shared_up_ch_all' not found
This time for each of the three donors: self-healing down vs. uninfected.
## Error in eval(expr, envir, enclos): object 'hs_uninf_filtv2_sig' not found
comp_lst <- list(
"donor_107" = rownames(down_sig[["d107_shun"]]),
"donor_108" = rownames(down_sig[["d108_shun"]]),
"donor_110" = rownames(down_sig[["d110_shun"]]))
## Error in rownames(down_sig[["d107_shun"]]): object 'down_sig' not found
## Error in Vennerable::Venn(Sets = comp_lst): object 'comp_lst' not found
## Error in Vennerable::plot(comp_shun_down, doWeights = FALSE): object 'comp_shun_down' not found
## Error in eval(expr, envir, enclos): object 'comp_shun_down' not found
## Error in eval(expr, envir, enclos): object 'hs_uninf_filtv2_combined' not found
## Error in eval(expr, envir, enclos): object 'hs_uninf_filtv2_combined' not found
## Error in eval(expr, envir, enclos): object 'hs_uninf_filtv2_combined' not found
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")
## Error in merge(de_table_shared_down_sh_first[, c("description", "deseq_logfc", : object 'de_table_shared_down_sh_first' not found
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")
## Error in merge(de_table_shared_down_sh_all, de_table_shared_down_sh_third[, : object 'de_table_shared_down_sh_all' not found
## Error in eval(expr, envir, enclos): object 'de_table_shared_down_sh_all' not found
## Error in eval(expr, envir, enclos): object 'de_table_shared_down_sh_all' not found
colnames(de_table_shared_down_sh_all) <- c("description", "logfc_107",
"adjp_107", "logfc_108", "adjp_108",
"logfc_110", "adjp_110")
## Error in colnames(de_table_shared_down_sh_all) <- c("description", "logfc_107", : object 'de_table_shared_down_sh_all' not found
write.csv(de_table_shared_down_sh_all,
file=paste0("images/", rundate, "_de_table_shared_down_sh_all.csv"))
## Error in is.data.frame(x): object 'de_table_shared_down_sh_all' not found
This time for each of the three donors: chronic down vs. uninfected.
## Error in eval(expr, envir, enclos): object 'hs_uninf_filtv2_sig' not found
comp_lst <- list(
"donor_107" = rownames(down_sig[["d107_chun"]]),
"donor_108" = rownames(down_sig[["d108_chun"]]),
"donor_110" = rownames(down_sig[["d110_chun"]]))
## Error in rownames(down_sig[["d107_chun"]]): object 'down_sig' not found
## Error in Vennerable::Venn(Sets = comp_lst): object 'comp_lst' not found
## Error in Vennerable::plot(comp_chun_down, doWeights = FALSE): object 'comp_chun_down' not found
## Error in eval(expr, envir, enclos): object 'comp_chun_down' not found
## Error in eval(expr, envir, enclos): object 'hs_uninf_filtv2_combined' not found
## Error in eval(expr, envir, enclos): object 'hs_uninf_filtv2_combined' not found
## Error in eval(expr, envir, enclos): object 'hs_uninf_filtv2_combined' not found
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")
## Error in merge(de_table_shared_down_ch_first[, c("description", "deseq_logfc", : object 'de_table_shared_down_ch_first' not found
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")
## Error in merge(de_table_shared_down_ch_all, de_table_shared_down_ch_third[, : object 'de_table_shared_down_ch_all' not found
## Error in eval(expr, envir, enclos): object 'de_table_shared_down_ch_all' not found
## Error in eval(expr, envir, enclos): object 'de_table_shared_down_ch_all' not found
colnames(de_table_shared_down_ch_all) <- c("description", "logfc_107",
"adjp_107", "logfc_108", "adjp_108",
"logfc_110", "adjp_110")
## Error in colnames(de_table_shared_down_ch_all) <- c("description", "logfc_107", : object 'de_table_shared_down_ch_all' not found
write.csv(de_table_shared_down_ch_all,
file=paste0("images/", rundate, "_de_table_shared_down_ch_all.csv"))
## Error in is.data.frame(x): object 'de_table_shared_down_ch_all' not found
At this point, we should have a set of genes which are up/down in the self/uninfected and chronic/uninfected, kept in variables with names like: ‘shared_shun_down’ and ‘shared_chun_down’
Now we want a sense of what genes are shared/unique among the self-healing vs. uninfected and the chronic vs. uninfected comparisons performed above. One would assume that the most interesting genes in these sets will prove to be the the ones which are not shared.
## Error in eval(expr, envir, enclos): object 'shared_shun_up' not found
## Error in Vennerable::Venn(Sets = shch_up_lst): object 'shch_up_lst' not found
## Error in Vennerable::plot(shared_up, doWeights = FALSE): object 'shared_up' not found
## Error in eval(expr, envir, enclos): object 'shared_shun_down' not found
## Error in Vennerable::Venn(Sets = shch_down_lst): object 'shch_down_lst' not found
## Error in Vennerable::plot(shared_down, doWeights = FALSE): object 'shared_down' not found
## Error in eval(expr, envir, enclos): object 'shared_up' not found
## Error in eval(expr, envir, enclos): object 'shared_up' not found
## Error in eval(expr, envir, enclos): object 'shared_up' not found
## Error in eval(expr, envir, enclos): object 'shared_down' not found
## Error in eval(expr, envir, enclos): object 'shared_down' not found
## Error in eval(expr, envir, enclos): object 'shared_down' not found
kept_columns <- c("ensembltranscriptid", "ensemblgeneid", "description",
"deseq_logfc", "deseq_adjp")
## Get the shared things in up_sh and up_ch, ergo can use either sheet
xls_result <- write_xls(data=up_sig[["d107_shun"]][upup_genes, kept_columns],
sheet="upsh_upch")
## Error in write_xls(data = up_sig[["d107_shun"]][upup_genes, kept_columns], : object 'up_sig' not found
## Get the only-up_sh
xls_result <- write_xls(data=up_sig[["d107_shun"]][upsh_notch, kept_columns],
sheet="upsh_noch",
wb=xls_result[["workbook"]])
## Error in write_xls(data = up_sig[["d107_shun"]][upsh_notch, kept_columns], : object 'up_sig' not found
## Get the only-up_ch
xls_result <- write_xls(data=up_sig[["d107_chun"]][upch_notsh, kept_columns],
sheet="upch_nosh",
wb=xls_result[["workbook"]])
## Error in write_xls(data = up_sig[["d107_chun"]][upch_notsh, kept_columns], : object 'up_sig' not found
## Now repeat for the down: down-shared
## Get shared down down_sh down_ch
xls_result <- write_xls(data=down_sig[["d107_shun"]][downdown_genes, kept_columns],
sheet="downsh_downch",
wb=xls_result[["workbook"]])
## Error in write_xls(data = down_sig[["d107_shun"]][downdown_genes, kept_columns], : object 'down_sig' not found
## The down_sh only
xls_result <- write_xls(data=down_sig[["d107_shun"]][downsh_notch, kept_columns],
sheet="downsh_noch",
wb=xls_result[["workbook"]])
## Error in write_xls(data = down_sig[["d107_shun"]][downsh_notch, kept_columns], : object 'down_sig' not found
## The down_ch
xls_result <- write_xls(data=down_sig[["d107_chun"]][downch_notsh, kept_columns],
sheet="downch_nosh",
wb=xls_result[["workbook"]],
excel=paste0("excel/", rundate, "_figure_5c_stuff-v", ver, ".xlsx"))
## Error in write_xls(data = down_sig[["d107_chun"]][downch_notsh, kept_columns], : object 'down_sig' not found
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_batches(expt=hs_cds_uninf, fact="pathogenstrain")
uninf_strainbatch_qcf <- set_expt_conditions(expt=uninf_strainbatch_qcf, fact="state")
withuninf_strainbatch_pairs_chsh <- all_pairwise(
parallel=FALSE, uninf_strainbatch_qcf,
model_batch=FALSE, force=TRUE, do_ebseq=FALSE)
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/", rundate, "_withuninf_strainbatch_chsh_pairwise-v", ver, ".xlsx")))
withuninf_strainbatch_sig_chsh <- extract_significant_genes(
withuninf_strainbatch_tables_chsh,
excel=paste0("excel/", rundate, "_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
Generate DE lists of each donor for all contrasts for PBMCs.
I renamed these plots and am now hopelessly confused as to which is which. I think I will not run this for now but instead generate the worksheet without them and then return to this in the hopes that I can do a better job.
pp(file=paste0("images/", rundate, "_fig_05a-sh_uninf_vs_chr_uninf_up.pdf"))
Vennerable::plot(common_solos_batch$up_venn, doWeights=FALSE)
dev.off()
pp(file=paste0("images/", rundate, "_fig_05b-sh_uninf_vs_chr_uninf_down.pdf"))
Vennerable::plot(common_solos_batch$down_venn, doWeights=FALSE)
dev.off()
pp(file=paste0("images/", rundate, "_fig_05c-sh_uninf_donors_up.pdf"))
Vennerable::plot(sh_up_venn, doWeights=FALSE)
dev.off()
pp(file=paste0("images/", rundate, "_fig_05d-sh_uninf_donors_down.pdf"))
Vennerable::plot(sh_down_venn, doWeights=FALSE)
dev.off()
pp(file=paste0("images/", rundate, "_fig_05e-chr_uninf_donors_up.pdf"))
Vennerable::plot(chr_up_venn, doWeights=FALSE)
dev.off()
pp(file=paste0("images/", rundate, "_fig_05f-chr_uninf_donors_down.pdf"))
Vennerable::plot(chr_down_venn, doWeights=FALSE)
dev.off()
pp(file=paste0("images/", rundate, "_fig_05g-up_common.pdf"))
Vennerable::plot(shared_up, doWeights=FALSE)
dev.off()
pp(file=paste0("images/", rundate, "_fig_05h-down_common.pdf"))
Vennerable::plot(shared_down, doWeights=FALSE)
dev.off()
lp_pairwise_nobatch <- sm(all_pairwise(lp_inf_filt, model_batch=FALSE,
parallel=FALSE, do_ebseq=FALSE))
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
excel_file <- glue::glue("excel/{rundate}_lp_infect_nobatch_contr-v{ver}.xlsx")
lp_combined_nobatch <- sm(combine_de_tables(lp_pairwise_nobatch, excel=excel_file))
## Error in combine_de_tables(lp_pairwise_nobatch, excel = excel_file): object 'lp_pairwise_nobatch' not found
excel_file <- glue::glue("excel/{rundate}_lp_infect_nobatch_sig-v{ver}.xlsx")
lp_sig_nobatch <- sm(extract_significant_genes(lp_combined_nobatch, excel=excel_file))
## Error in extract_significant_genes(lp_combined_nobatch, excel = excel_file): object 'lp_combined_nobatch' not found
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
excel_file <- glue::glue("excel/{rundate}_lp_infect_batch_contr-v{ver}.xlsx")
lp_combined_batch <- sm(combine_de_tables(lp_pairwise_batch, excel=excel_file))
## Error in combine_de_tables(lp_pairwise_batch, excel = excel_file): object 'lp_pairwise_batch' not found
excel_file <- glue::glue("excel/{rundate}_lp_infect_batch_sig-v{ver}.xlsx")
lp_sig_batch <- sm(extract_significant_genes(lp_combined_batch, excel=excel_file))
## Error in extract_significant_genes(lp_combined_batch, excel = excel_file): object 'lp_combined_batch' not found
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
excel_file <- glue::glue("excel/{rundate}_lp_infect_ssva_contr-v{ver}.xlsx")
lp_combined_ssva <- sm(combine_de_tables(lp_pairwise_ssva, excel=excel_file))
## Error in combine_de_tables(lp_pairwise_ssva, excel = excel_file): object 'lp_pairwise_ssva' not found
excel_file <- glue::glue("excel/{rundate}_lp_infect_ssva_sig-v{ver}.xlsx")
lp_sig_ssva <- sm(extract_significant_genes(lp_combined_ssva, excel=excel_file))
## Error in extract_significant_genes(lp_combined_ssva, excel = excel_file): object 'lp_combined_ssva' not found
## Error in preprocessCore::normalize.quantiles(as.matrix(count_table), copy = TRUE): ERROR; return code from pthread_create() is 22
excel_file <- glue::glue("excel/{rundate}_lp_infect_fsva_contr-v{ver}.xlsx")
lp_combined_fsva <- sm(combine_de_tables(lp_pairwise_fsva, excel=excel_file))
## Error in combine_de_tables(lp_pairwise_fsva, excel = excel_file): object 'lp_pairwise_fsva' not found
excel_file <- glue::glue("excel/{rundate}_lp_infect_fsva_sig-v{ver}.xlsx")
lp_sig_fsva <- sm(extract_significant_genes(lp_combined_fsva, excel=excel_file))
## Error in extract_significant_genes(lp_combined_fsva, excel = excel_file): object 'lp_combined_fsva' not found
R version 3.5.1 (2018-07-02)
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, stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: bindrcpp(v.0.2.2), ruv(v.0.9.7), variancePartition(v.1.12.0), scales(v.1.0.0), foreach(v.1.4.4), limma(v.3.38.3), ggplot2(v.3.1.0), hpgltools(v.2018.11), Biobase(v.2.42.0) and BiocGenerics(v.0.28.0)
loaded via a namespace (and not attached): tidyselect(v.0.2.5), lme4(v.1.1-19), RSQLite(v.2.1.1), AnnotationDbi(v.1.44.0), htmlwidgets(v.1.3), grid(v.3.5.1), BiocParallel(v.1.16.2), devtools(v.2.0.1), munsell(v.0.5.0), codetools(v.0.2-15), preprocessCore(v.1.44.0), units(v.0.6-2), withr(v.2.1.2), colorspace(v.1.3-2), GOSemSim(v.2.8.0), OrganismDbi(v.1.24.0), highr(v.0.7), knitr(v.1.20), rstudioapi(v.0.8), stats4(v.3.5.1), Vennerable(v.3.1.0.9000), DOSE(v.3.8.0), labeling(v.0.3), urltools(v.1.7.1), GenomeInfoDbData(v.1.2.0), bit64(v.0.9-7), farver(v.1.1.0), rprojroot(v.1.3-2), R6(v.2.3.0), doParallel(v.1.0.14), GenomeInfoDb(v.1.18.1), locfit(v.1.5-9.1), bitops(v.1.0-6), fgsea(v.1.8.0), gridGraphics(v.0.3-0), DelayedArray(v.0.8.0), assertthat(v.0.2.0), ggraph(v.1.0.2), nnet(v.7.3-12), enrichplot(v.1.2.0), gtable(v.0.2.0), sva(v.3.30.0), processx(v.3.2.1), rlang(v.0.3.0.1), genefilter(v.1.64.0), splines(v.3.5.1), rtracklayer(v.1.42.1), lazyeval(v.0.2.1), acepack(v.1.4.1), checkmate(v.1.8.5), europepmc(v.0.3), BiocManager(v.1.30.4), yaml(v.2.2.0), reshape2(v.1.4.3), GenomicFeatures(v.1.34.1), backports(v.1.1.2), qvalue(v.2.14.0), Hmisc(v.4.1-1), clusterProfiler(v.3.10.0), RBGL(v.1.58.1), tools(v.3.5.1), usethis(v.1.4.0), ggplotify(v.0.0.3), gplots(v.3.0.1), RColorBrewer(v.1.1-2), sessioninfo(v.1.1.1), ggridges(v.0.5.1), Rcpp(v.1.0.0), plyr(v.1.8.4), base64enc(v.0.1-3), progress(v.1.2.0), zlibbioc(v.1.28.0), purrr(v.0.2.5), RCurl(v.1.95-4.11), ps(v.1.2.1), prettyunits(v.1.0.2), rpart(v.4.1-13), viridis(v.0.5.1), cowplot(v.0.9.3), S4Vectors(v.0.20.1), SummarizedExperiment(v.1.12.0), ggrepel(v.0.8.0), cluster(v.2.0.7-1), colorRamps(v.2.3), fs(v.1.2.6), magrittr(v.1.5), data.table(v.1.11.8), openxlsx(v.4.1.0), DO.db(v.2.9), triebeard(v.0.3.0), packrat(v.0.5.0), matrixStats(v.0.54.0), pkgload(v.1.0.2), hms(v.0.4.2), evaluate(v.0.12), xtable(v.1.8-3), pbkrtest(v.0.4-7), XML(v.3.98-1.16), IRanges(v.2.16.0), gridExtra(v.2.3), testthat(v.2.0.1), compiler(v.3.5.1), biomaRt(v.2.38.0), tibble(v.1.4.2), KernSmooth(v.2.23-15), crayon(v.1.3.4), minqa(v.1.2.4), htmltools(v.0.3.6), corpcor(v.1.6.9), mgcv(v.1.8-26), Formula(v.1.2-3), geneplotter(v.1.60.0), tidyr(v.0.8.2), DBI(v.1.0.0), tweenr(v.1.0.0), MASS(v.7.3-51.1), Matrix(v.1.2-15), cli(v.1.0.1), quadprog(v.1.5-5), gdata(v.2.18.0), bindr(v.0.1.1), igraph(v.1.2.2), GenomicRanges(v.1.34.0), pkgconfig(v.2.0.2), rvcheck(v.0.1.3), GenomicAlignments(v.1.18.0), foreign(v.0.8-71), plotly(v.4.8.0), xml2(v.1.2.0), annotate(v.1.60.0), XVector(v.0.22.0), stringr(v.1.3.1), callr(v.3.0.0), digest(v.0.6.18), graph(v.1.60.0), Biostrings(v.2.50.1), rmarkdown(v.1.10), fastmatch(v.1.1-0), htmlTable(v.1.12), edgeR(v.3.24.1), directlabels(v.2018.05.22), Rsamtools(v.1.34.0), gtools(v.3.8.1), nloptr(v.1.2.1), nlme(v.3.1-137), jsonlite(v.1.6), desc(v.1.2.0), viridisLite(v.0.3.0), pillar(v.1.3.0), lattice(v.0.20-38), httr(v.1.3.1), pkgbuild(v.1.0.2), survival(v.2.43-3), GO.db(v.3.7.0), glue(v.1.3.0), remotes(v.2.0.2), zip(v.1.0.0), UpSetR(v.1.3.3), iterators(v.1.0.10), pander(v.0.6.3), bit(v.1.1-14), ggforce(v.0.1.3), stringi(v.1.2.4), blob(v.1.1.1), DESeq2(v.1.22.1), latticeExtra(v.0.6-28), caTools(v.1.17.1.1), memoise(v.1.1.0) and dplyr(v.0.7.8)
## If you wish to reproduce this exact build of hpgltools, invoke the following:
## > git clone http://github.com/abelew/hpgltools.git
## > git reset 242609d9ad73083ea13099760f8f0678a4befbf8
## This is hpgltools commit: Wed Dec 5 15:02:17 2018 -0500: 242609d9ad73083ea13099760f8f0678a4befbf8
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_20180822-v20180822.rda.xz
this_save