I want to perform a series of comparisons among the host cells: human and mouse. Thus I need to collect annotation data for both species and get the set of orthologs between them.
## The biomart annotations file already exists, loading from it.
rownames(hs_annot) <- make.names(
paste0(hs_annot[["ensembl_transcript_id"]], ".",
hs_annot[["transcript_version"]]),
unique=TRUE)
hs_tx_gene <- hs_annot[, c("ensembl_gene_id", "ensembl_transcript_id")]
hs_tx_gene[["id"]] <- rownames(hs_tx_gene)
hs_tx_gene <- hs_tx_gene[, c("id", "ensembl_gene_id")]
new_hs_annot <- hs_annot
rownames(new_hs_annot) <- make.names(hs_annot[["ensembl_gene_id"]], unique=TRUE)
The question is reasonably self-contained. I want to compare the uninfected human samples against any samples which were infected for 4 hours. So let us first pull those samples and then poke at them a bit.
sample_sheet <- "sample_sheets/leishmania_host_metasheet_20190401.xlsx"
hs_expt <- create_expt(sample_sheet,
file_column="hsapiensfile",
gene_info=new_hs_annot,
tx_gene_map=hs_tx_gene)
## Reading the sample metadata.
## The sample definitions comprises: 437 rows(samples) and 56 columns(metadata fields).
## Reading count tables.
## Using the transcript to gene mapping.
## Reading salmon data with tximport.
## Finished reading count tables.
## Matched 19629 annotations and counts.
## Bringing together the count matrix and gene information.
## The mapped IDs are not the rownames of your gene information, changing them now.
## Some annotations were lost in merging, setting them to 'undefined'.
## Saving the expressionset to 'expt.rda'.
## The final expressionset has 19629 rows and 267 columns.
## Using a subset expression.
## There were 267, now there are 247 samples.
## Using a subset expression.
## There were 247, now there are 64 samples.
hs_t4h_expt <- set_expt_conditions(hs_t4h_expt, fact="infectstate")
hs_t4h_expt <- set_expt_batches(hs_t4h_expt, fact="study")
table(hs_t4h_expt$conditions)
##
## bead no stim yes
## 3 18 35 8
##
## lps-timecourse m-gm-csf mbio
## 8 39 17
hs_t4h_norm <- normalize_expt(hs_t4h_expt, norm="quant", convert="cpm",
transform="log2", filter=TRUE, batch="svaseq")
## This function will replace the expt$expressionset slot with:
## log2(svaseq(cpm(quant(cbcb(data)))))
## It will save copies of each step along the way
## in expt$normalized with the corresponding libsizes. Keep libsizes in mind
## when invoking limma. The appropriate libsize is non-log(cpm(normalized)).
## This is most likely kept at:
## 'new_expt$normalized$intermediate_counts$normalization$libsizes'
## A copy of this may also be found at:
## new_expt$best_libsize
## Warning in normalize_expt(hs_t4h_expt, norm = "quant", convert = "cpm", :
## Quantile normalization and sva do not always play well together.
## Step 1: performing count filter with option: cbcb
## Removing 7360 low-count genes (12269 remaining).
## Step 2: normalizing the data with quant.
## Step 3: converting the data with cpm.
## Step 4: transforming the data with log2.
## transform_counts: Found 3822 values equal to 0, adding 1 to the matrix.
## Step 5: doing batch correction with svaseq.
## Note to self: If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.
## Passing off to all_adjusters.
## batch_counts: Before batch/surrogate estimation, 711267 entries are x>1: 91%.
## batch_counts: Before batch/surrogate estimation, 3822 entries are x==0: 0%.
## batch_counts: Before batch/surrogate estimation, 70127 entries are 0<x<1: 9%.
## The be method chose 12 surrogate variable(s).
## Attempting svaseq estimation with 12 surrogates.
## There are 1769 (0%) elements which are < 0 after batch correction.
## Setting low elements to zero.
## Using a subset expression.
## There were 64, now there are 29 samples.
hs_t4h_inf_norm <- normalize_expt(hs_t4h_inf, transform="log2", convert="cpm",
filter=TRUE, batch="svaseq")
## This function will replace the expt$expressionset slot with:
## log2(svaseq(cpm(cbcb(data))))
## It will save copies of each step along the way
## in expt$normalized with the corresponding libsizes. Keep libsizes in mind
## when invoking limma. The appropriate libsize is non-log(cpm(normalized)).
## This is most likely kept at:
## 'new_expt$normalized$intermediate_counts$normalization$libsizes'
## A copy of this may also be found at:
## new_expt$best_libsize
## Leaving the data unnormalized. This is necessary for DESeq, but
## EdgeR/limma might benefit from normalization. Good choices include quantile,
## size-factor, tmm, etc.
## Step 1: performing count filter with option: cbcb
## Removing 7759 low-count genes (11870 remaining).
## Step 2: not normalizing the data.
## Step 3: converting the data with cpm.
## Step 4: transforming the data with log2.
## transform_counts: Found 3276 values equal to 0, adding 1 to the matrix.
## Step 5: doing batch correction with svaseq.
## Note to self: If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.
## Passing off to all_adjusters.
## batch_counts: Before batch/surrogate estimation, 319013 entries are x>1: 93%.
## batch_counts: Before batch/surrogate estimation, 3276 entries are x==0: 1%.
## batch_counts: Before batch/surrogate estimation, 21941 entries are 0<x<1: 6%.
## The be method chose 6 surrogate variable(s).
## Attempting svaseq estimation with 6 surrogates.
## There are 557 (0%) elements which are < 0 after batch correction.
## Setting low elements to zero.
keepers <- list("infection" = c("yes", "no"))
hs_t4h_de <- all_pairwise(hs_t4h_inf, model_batch="svaseq", filter=TRUE, force=TRUE)
## batch_counts: Before batch/surrogate estimation, 339179 entries are x>1: 99%.
## batch_counts: Before batch/surrogate estimation, 3276 entries are x==0: 1%.
## batch_counts: Before batch/surrogate estimation, 216 entries are 0<x<1: 0%.
## The be method chose 4 surrogate variable(s).
## Attempting svaseq estimation with 4 surrogates.
## Plotting a PCA before surrogates/batch inclusion.
## Not putting labels on the plot.
## Using svaseq to visualize before/after batch inclusion.
## Performing a test normalization with: raw
## This function will replace the expt$expressionset slot with:
## log2(svaseq(cpm(cbcb(data))))
## It will save copies of each step along the way
## in expt$normalized with the corresponding libsizes. Keep libsizes in mind
## when invoking limma. The appropriate libsize is non-log(cpm(normalized)).
## This is most likely kept at:
## 'new_expt$normalized$intermediate_counts$normalization$libsizes'
## A copy of this may also be found at:
## new_expt$best_libsize
## Leaving the data unnormalized. This is necessary for DESeq, but
## EdgeR/limma might benefit from normalization. Good choices include quantile,
## size-factor, tmm, etc.
## Step 1: performing count filter with option: cbcb
## Removing 0 low-count genes (11870 remaining).
## Step 2: not normalizing the data.
## Step 3: converting the data with cpm.
## Step 4: transforming the data with log2.
## transform_counts: Found 3276 values equal to 0, adding 1 to the matrix.
## Step 5: doing batch correction with svaseq.
## Note to self: If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.
## Passing off to all_adjusters.
## batch_counts: Before batch/surrogate estimation, 319013 entries are x>1: 93%.
## batch_counts: Before batch/surrogate estimation, 3276 entries are x==0: 1%.
## batch_counts: Before batch/surrogate estimation, 21941 entries are 0<x<1: 6%.
## The be method chose 6 surrogate variable(s).
## Attempting svaseq estimation with 6 surrogates.
## There are 557 (0%) elements which are < 0 after batch correction.
## Setting low elements to zero.
## Not putting labels on the plot.
## Finished running DE analyses, collecting outputs.
## Comparing analyses.
hs_t4h_table <- combine_de_tables(hs_t4h_de, keepers=keepers,
excel=glue::glue("excel/hs_t4h_de_table-v{ver}.xlsx"))
## Writing a legend of columns.
## Printing a pca plot before/after surrogates/batch estimation.
## Working on 1/1: infection which is: yes/no.
## Found table with yes_vs_no
## Adding venn plots for infection.
## Limma expression coefficients for infection; R^2: 0.984; equation: y = 0.981x + 0.0934
## Deseq expression coefficients for infection; R^2: 0.982; equation: y = 0.96x + 0.37
## Edger expression coefficients for infection; R^2: 0.982; equation: y = 0.961x + 0.377
## Writing summary information, compare_plot is: TRUE.
## Performing save of excel/hs_t4h_de_table-v20190326.xlsx.
## Writing a legend of columns.
## Writing excel data according to limma for infection: 1/5.
## After (adj)p filter, the up genes table has 2813 genes.
## After (adj)p filter, the down genes table has 3963 genes.
## After fold change filter, the up genes table has 1032 genes.
## After fold change filter, the down genes table has 793 genes.
## Printing significant genes to the file: excel/significant_genes.xlsx
## 1/1: Creating significant table up_limma_infection
## Writing excel data according to edger for infection: 1/5.
## After (adj)p filter, the up genes table has 3277 genes.
## After (adj)p filter, the down genes table has 3261 genes.
## After fold change filter, the up genes table has 1179 genes.
## After fold change filter, the down genes table has 1022 genes.
## Printing significant genes to the file: excel/significant_genes.xlsx
## 1/1: Creating significant table up_edger_infection
## Writing excel data according to deseq for infection: 1/5.
## After (adj)p filter, the up genes table has 3454 genes.
## After (adj)p filter, the down genes table has 3607 genes.
## After fold change filter, the up genes table has 1165 genes.
## After fold change filter, the down genes table has 1040 genes.
## Printing significant genes to the file: excel/significant_genes.xlsx
## 1/1: Creating significant table up_deseq_infection
## Writing excel data according to ebseq for infection: 1/5.
## After (adj)p filter, the up genes table has 1893 genes.
## After (adj)p filter, the down genes table has 1940 genes.
## After fold change filter, the up genes table has 895 genes.
## After fold change filter, the down genes table has 1031 genes.
## Printing significant genes to the file: excel/significant_genes.xlsx
## 1/1: Creating significant table up_ebseq_infection
## Writing excel data according to basic for infection: 1/5.
## After (adj)p filter, the up genes table has 1335 genes.
## After (adj)p filter, the down genes table has 2592 genes.
## After fold change filter, the up genes table has 909 genes.
## After fold change filter, the down genes table has 1361 genes.
## Printing significant genes to the file: excel/significant_genes.xlsx
## 1/1: Creating significant table up_basic_infection
## Adding significance bar plots.
sig_up <- hs_t4h_sig$deseq$ups[[1]]
sig_down <- hs_t4h_sig$deseq$downs[[1]]
hs_t4h_gprofiler_up <- simple_gprofiler(sig_up, species="hsapiens")
## Performing gProfiler GO search of 1165 genes against hsapiens.
## GO search found 1047 hits.
## Performing gProfiler KEGG search of 1165 genes against hsapiens.
## KEGG search found 67 hits.
## Performing gProfiler REAC search of 1165 genes against hsapiens.
## REAC search found 57 hits.
## Performing gProfiler MI search of 1165 genes against hsapiens.
## MI search found 0 hits.
## Performing gProfiler TF search of 1165 genes against hsapiens.
## TF search found 405 hits.
## Performing gProfiler CORUM search of 1165 genes against hsapiens.
## CORUM search found 12 hits.
## Performing gProfiler HP search of 1165 genes against hsapiens.
## HP search found 3 hits.
## Performing gProfiler GO search of 1040 genes against hsapiens.
## GO search found 62 hits.
## Performing gProfiler KEGG search of 1040 genes against hsapiens.
## KEGG search found 6 hits.
## Performing gProfiler REAC search of 1040 genes against hsapiens.
## REAC search found 36 hits.
## Performing gProfiler MI search of 1040 genes against hsapiens.
## MI search found 0 hits.
## Performing gProfiler TF search of 1040 genes against hsapiens.
## TF search found 44 hits.
## Performing gProfiler CORUM search of 1040 genes against hsapiens.
## CORUM search found 2 hits.
## Performing gProfiler HP search of 1040 genes against hsapiens.
## HP search found 5 hits.
I want to perform a series of comparisons among the host cells: human and mouse. Thus I need to collect annotation data for both species and get the set of orthologs between them.
## The biomart annotations file already exists, loading from it.
rownames(mm_annot) <- make.names(
paste0(mm_annot[["ensembl_transcript_id"]], ".",
mm_annot[["transcript_version"]]),
unique=TRUE)
mm_tx_gene <- mm_annot[, c("ensembl_gene_id", "ensembl_transcript_id")]
mm_tx_gene[["id"]] <- rownames(mm_tx_gene)
mm_tx_gene <- mm_tx_gene[, c("id", "ensembl_gene_id")]
new_mm_annot <- mm_annot
rownames(new_mm_annot) <- make.names(mm_annot[["ensembl_gene_id"]], unique=TRUE)
The question is reasonably self-contained. I want to compare the uninfected human samples against any samples which were infected for 4 hours. So let us first pull those samples and then poke at them a bit.
mm_expt <- create_expt(sample_sheet,
file_column="mmusculusfile",
gene_info=new_mm_annot,
tx_gene_map=mm_tx_gene)
## Reading the sample metadata.
## The sample definitions comprises: 437 rows(samples) and 56 columns(metadata fields).
## Reading count tables.
## Using the transcript to gene mapping.
## Reading salmon data with tximport.
## Finished reading count tables.
## Matched 19660 annotations and counts.
## Bringing together the count matrix and gene information.
## The mapped IDs are not the rownames of your gene information, changing them now.
## Some annotations were lost in merging, setting them to 'undefined'.
## Saving the expressionset to 'expt.rda'.
## The final expressionset has 19660 rows and 105 columns.
## Using a subset expression.
## There were 105, now there are 41 samples.
##
## no stim yes
## 11 24 6
##
## undefined
## 41
mm_t4h_norm <- normalize_expt(mm_t4h_expt, norm="quant", convert="cpm",
transform="log2", filter=TRUE, batch="svaseq")
## This function will replace the expt$expressionset slot with:
## log2(svaseq(cpm(quant(cbcb(data)))))
## It will save copies of each step along the way
## in expt$normalized with the corresponding libsizes. Keep libsizes in mind
## when invoking limma. The appropriate libsize is non-log(cpm(normalized)).
## This is most likely kept at:
## 'new_expt$normalized$intermediate_counts$normalization$libsizes'
## A copy of this may also be found at:
## new_expt$best_libsize
## Warning in normalize_expt(mm_t4h_expt, norm = "quant", convert = "cpm", :
## Quantile normalization and sva do not always play well together.
## Step 1: performing count filter with option: cbcb
## Removing 9350 low-count genes (10310 remaining).
## Step 2: normalizing the data with quant.
## Step 3: converting the data with cpm.
## Step 4: transforming the data with log2.
## transform_counts: Found 38 values equal to 0, adding 1 to the matrix.
## Step 5: doing batch correction with svaseq.
## Note to self: If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.
## Passing off to all_adjusters.
## batch_counts: Before batch/surrogate estimation, 390851 entries are x>1: 92%.
## batch_counts: Before batch/surrogate estimation, 38 entries are x==0: 0%.
## batch_counts: Before batch/surrogate estimation, 31821 entries are 0<x<1: 8%.
## The be method chose 8 surrogate variable(s).
## Attempting svaseq estimation with 8 surrogates.
## There are 809 (0%) elements which are < 0 after batch correction.
## Setting low elements to zero.
mm_t4h_inf_norm <- normalize_expt(mm_t4h_expt, transform="log2", convert="cpm",
filter=TRUE, batch="svaseq")
## This function will replace the expt$expressionset slot with:
## log2(svaseq(cpm(cbcb(data))))
## It will save copies of each step along the way
## in expt$normalized with the corresponding libsizes. Keep libsizes in mind
## when invoking limma. The appropriate libsize is non-log(cpm(normalized)).
## This is most likely kept at:
## 'new_expt$normalized$intermediate_counts$normalization$libsizes'
## A copy of this may also be found at:
## new_expt$best_libsize
## Leaving the data unnormalized. This is necessary for DESeq, but
## EdgeR/limma might benefit from normalization. Good choices include quantile,
## size-factor, tmm, etc.
## Step 1: performing count filter with option: cbcb
## Removing 9350 low-count genes (10310 remaining).
## Step 2: not normalizing the data.
## Step 3: converting the data with cpm.
## Step 4: transforming the data with log2.
## transform_counts: Found 3253 values equal to 0, adding 1 to the matrix.
## Step 5: doing batch correction with svaseq.
## Note to self: If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.
## Passing off to all_adjusters.
## batch_counts: Before batch/surrogate estimation, 385615 entries are x>1: 91%.
## batch_counts: Before batch/surrogate estimation, 3253 entries are x==0: 1%.
## batch_counts: Before batch/surrogate estimation, 33842 entries are 0<x<1: 8%.
## The be method chose 7 surrogate variable(s).
## Attempting svaseq estimation with 7 surrogates.
## There are 1186 (0%) elements which are < 0 after batch correction.
## Setting low elements to zero.
## batch_counts: Before batch/surrogate estimation, 417208 entries are x>1: 99%.
## batch_counts: Before batch/surrogate estimation, 3253 entries are x==0: 1%.
## batch_counts: Before batch/surrogate estimation, 392 entries are 0<x<1: 0%.
## The be method chose 5 surrogate variable(s).
## Attempting svaseq estimation with 5 surrogates.
## Plotting a PCA before surrogates/batch inclusion.
## Not putting labels on the plot.
## Using svaseq to visualize before/after batch inclusion.
## Performing a test normalization with: raw
## This function will replace the expt$expressionset slot with:
## log2(svaseq(cpm(cbcb(data))))
## It will save copies of each step along the way
## in expt$normalized with the corresponding libsizes. Keep libsizes in mind
## when invoking limma. The appropriate libsize is non-log(cpm(normalized)).
## This is most likely kept at:
## 'new_expt$normalized$intermediate_counts$normalization$libsizes'
## A copy of this may also be found at:
## new_expt$best_libsize
## Leaving the data unnormalized. This is necessary for DESeq, but
## EdgeR/limma might benefit from normalization. Good choices include quantile,
## size-factor, tmm, etc.
## Step 1: performing count filter with option: cbcb
## Removing 0 low-count genes (10310 remaining).
## Step 2: not normalizing the data.
## Step 3: converting the data with cpm.
## Step 4: transforming the data with log2.
## transform_counts: Found 3253 values equal to 0, adding 1 to the matrix.
## Step 5: doing batch correction with svaseq.
## Note to self: If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.
## Passing off to all_adjusters.
## batch_counts: Before batch/surrogate estimation, 385615 entries are x>1: 91%.
## batch_counts: Before batch/surrogate estimation, 3253 entries are x==0: 1%.
## batch_counts: Before batch/surrogate estimation, 33842 entries are 0<x<1: 8%.
## The be method chose 7 surrogate variable(s).
## Attempting svaseq estimation with 7 surrogates.
## There are 1186 (0%) elements which are < 0 after batch correction.
## Setting low elements to zero.
## Not putting labels on the plot.
## Finished running DE analyses, collecting outputs.
## Comparing analyses.
mm_t4h_table <- combine_de_tables(mm_t4h_de, keepers=keepers,
excel=glue::glue("excel/mm_t4h_de_table-v{ver}.xlsx"))
## Deleting the file excel/mm_t4h_de_table-v20190326.xlsx before writing the tables.
## Writing a legend of columns.
## Printing a pca plot before/after surrogates/batch estimation.
## Working on 1/1: infection which is: yes/no.
## Found table with yes_vs_no
## Adding venn plots for infection.
## Limma expression coefficients for infection; R^2: 0.862; equation: y = 0.912x + 0.297
## Deseq expression coefficients for infection; R^2: 0.857; equation: y = 1.04x - 0.346
## Edger expression coefficients for infection; R^2: 0.858; equation: y = 1.04x - 0.373
## Writing summary information, compare_plot is: TRUE.
## Performing save of excel/mm_t4h_de_table-v20190326.xlsx.
## Writing a legend of columns.
## Writing excel data according to limma for infection: 1/5.
## After (adj)p filter, the up genes table has 2234 genes.
## After (adj)p filter, the down genes table has 2524 genes.
## After fold change filter, the up genes table has 825 genes.
## After fold change filter, the down genes table has 895 genes.
## Printing significant genes to the file: excel/significant_genes.xlsx
## 1/1: Creating significant table up_limma_infection
## Writing excel data according to edger for infection: 1/5.
## After (adj)p filter, the up genes table has 2153 genes.
## After (adj)p filter, the down genes table has 1967 genes.
## After fold change filter, the up genes table has 924 genes.
## After fold change filter, the down genes table has 911 genes.
## Printing significant genes to the file: excel/significant_genes.xlsx
## 1/1: Creating significant table up_edger_infection
## Writing excel data according to deseq for infection: 1/5.
## After (adj)p filter, the up genes table has 2474 genes.
## After (adj)p filter, the down genes table has 2404 genes.
## After fold change filter, the up genes table has 929 genes.
## After fold change filter, the down genes table has 932 genes.
## Printing significant genes to the file: excel/significant_genes.xlsx
## 1/1: Creating significant table up_deseq_infection
## Writing excel data according to ebseq for infection: 1/5.
## After (adj)p filter, the up genes table has 1149 genes.
## After (adj)p filter, the down genes table has 1310 genes.
## After fold change filter, the up genes table has 644 genes.
## After fold change filter, the down genes table has 710 genes.
## Printing significant genes to the file: excel/significant_genes.xlsx
## 1/1: Creating significant table up_ebseq_infection
## Writing excel data according to basic for infection: 1/5.
## After (adj)p filter, the up genes table has 2 genes.
## After (adj)p filter, the down genes table has 3 genes.
## After fold change filter, the up genes table has 2 genes.
## After fold change filter, the down genes table has 3 genes.
## Printing significant genes to the file: excel/significant_genes.xlsx
## 1/1: Creating significant table up_basic_infection
## Adding significance bar plots.
sig_up <- mm_t4h_sig$deseq$ups[[1]]
sig_down <- mm_t4h_sig$deseq$downs[[1]]
mm_t4h_gprofiler_up <- simple_gprofiler(sig_up, species="mmusculus")
## Performing gProfiler GO search of 929 genes against mmusculus.
## GO search found 841 hits.
## Performing gProfiler KEGG search of 929 genes against mmusculus.
## KEGG search found 36 hits.
## Performing gProfiler REAC search of 929 genes against mmusculus.
## REAC search found 30 hits.
## Performing gProfiler MI search of 929 genes against mmusculus.
## MI search found 0 hits.
## Performing gProfiler TF search of 929 genes against mmusculus.
## TF search found 439 hits.
## Performing gProfiler CORUM search of 929 genes against mmusculus.
## CORUM search found 4 hits.
## Performing gProfiler HP search of 929 genes against mmusculus.
## HP search found 7 hits.
## Performing gProfiler GO search of 932 genes against mmusculus.
## GO search found 51 hits.
## Performing gProfiler KEGG search of 932 genes against mmusculus.
## KEGG search found 8 hits.
## Performing gProfiler REAC search of 932 genes against mmusculus.
## REAC search found 9 hits.
## Performing gProfiler MI search of 932 genes against mmusculus.
## MI search found 0 hits.
## Performing gProfiler TF search of 932 genes against mmusculus.
## TF search found 331 hits.
## Performing gProfiler CORUM search of 932 genes against mmusculus.
## CORUM search found 2 hits.
## Performing gProfiler HP search of 932 genes against mmusculus.
## HP search found 2 hits.
Let us see if our human differential expression result is similar to that obtained in Table S2.
previous_hs <- readxl::read_excel("excel/inline-supplementary-material-5.xls", sheet=2)
previous_hs_lfc <- previous_hs[, c("ID", "Fold change")]
neg_idx <- previous_hs_lfc[[2]] < 0
previous_hs_lfc[neg_idx, 2] <- -1 * (1 / previous_hs_lfc[neg_idx, 2])
previous_hs_lfc[[2]] <- log2(previous_hs_lfc[[2]])
merged <- merge(previous_hs_lfc, hs_t4h_table$data[[1]], by.x="ID", by.y="row.names")
cor.test(merged[["limma_logfc"]], merged[["Fold change"]])
##
## Pearson's product-moment correlation
##
## data: merged[["limma_logfc"]] and merged[["Fold change"]]
## t = 148, df = 5082, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.8960 0.9063
## sample estimates:
## cor
## 0.9013
previous_mm <- readxl::read_excel("excel/12864_2015_2237_MOESM3_ESM.xls", sheet=2, skip=1)
previous_mm_lfc <- previous_mm[, c("ID", "Fold change")]
neg_idx <- previous_mm_lfc[[2]] < 0
previous_mm_lfc[neg_idx, 2] <- -1 * (1 / previous_mm_lfc[neg_idx, 2])
previous_mm_lfc[[2]] <- log2(previous_mm_lfc[[2]])
merged <- merge(previous_mm_lfc, mm_t4h_table$data[[1]], by.x="ID", by.y="row.names")
cor.test(merged[["limma_logfc"]], merged[["Fold change"]])
##
## Pearson's product-moment correlation
##
## data: merged[["limma_logfc"]] and merged[["Fold change"]]
## t = 223, df = 5718, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9445 0.9498
## sample estimates:
## cor
## 0.9472
First get the set of orthologs.
Side note: a different way of addressing this question resides in 20190220_host_comparisons.Rmd.
## The defaults of this function are suitable for mouse/human queries.
mm_hs_ortho <- load_biomart_orthologs()$all_linked_genes
mm_table <- mm_t4h_table$data[[1]]
hs_table <- hs_t4h_table$data[[1]]
mm_table <- merge(mm_hs_ortho, mm_table, by.x="mmusculus", by.y="row.names", all.y=TRUE)
hs_table <- merge(mm_hs_ortho, hs_table, by.x="hsapiens", by.y="row.names", all.y=TRUE)
both_table_hs <- merge(hs_table, mm_table, by.x="hsapiens", by.y="hsapiens")
both_table_mm <- merge(hs_table, mm_table, by.x="mmusculus", by.y="mmusculus")
cor.test(both_table_hs[["limma_logfc.x"]], both_table_hs[["limma_logfc.y"]])
##
## Pearson's product-moment correlation
##
## data: both_table_hs[["limma_logfc.x"]] and both_table_hs[["limma_logfc.y"]]
## t = 35, df = 13699, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2701 0.3008
## sample estimates:
## cor
## 0.2855
##
## Pearson's product-moment correlation
##
## data: both_table_mm[["limma_logfc.x"]] and both_table_mm[["limma_logfc.y"]]
## t = 35, df = 21859, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2160 0.2412
## sample estimates:
## cor
## 0.2286
Table S7 has a set of genes from human which are also up-regulated in mouse upon infection.
fig_s7_up <- readxl::read_excel("excel/inline-supplementary-material-7.xls", sheet=3)
fig_s7_hs_up <- unique(fig_s7_up[[5]])
fig_s7_down <- readxl::read_excel("excel/inline-supplementary-material-7.xls", sheet=4)
fig_s7_hs_down <- unique(fig_s7_down[[5]])
both_up_idx <- both_table_hs[["limma_logfc.x"]] >= 0.8 &
both_table_hs[["limma_logfc.y"]] >= 0.8 &
both_table_hs[["limma_adjp.x"]] <= 0.1 &
both_table_hs[["limma_adjp.y"]] <= 0.1
both_up_ids <- both_table_hs[both_up_idx, "hsapiens"]
up_venn <- Vennerable::Venn(Sets=list("figs7" = fig_s7_hs_up, "tables" = both_up_ids))
Vennerable::plot(up_venn)
both_down_idx <- both_table_hs[["limma_logfc.x"]] <= -0.8 &
both_table_hs[["limma_logfc.y"]] <= -0.8 &
both_table_hs[["limma_adjp.x"]] <= 0.1 &
both_table_hs[["limma_adjp.y"]] <= 0.1
both_down_ids <- both_table_hs[both_down_idx, "hsapiens"]
down_venn <- Vennerable::Venn(Sets=list("figs7" = fig_s7_hs_down, "tables" = both_down_ids))
Vennerable::plot(down_venn)
All of the neutrophil data is in mouse, apparently. This will make it more difficult, perhaps impossible to get an accurate answer.
So instead, look at infection vs. uninfected in mouse and then compare to the earliest Sacks’ timepoints in neutrophils.
subset <- "(hostcelltype=='PMN'&host=='mus_musculus'&expttime=='t12h') |
(hostcelltype=='macrophage'&host=='mus_musculus')"
neut_macr_mus <- subset_expt(mm_expt, subset=subset)
## Using a subset expression.
## There were 105, now there are 80 samples.
## Using a subset expression.
## There were 80, now there are 56 samples.
neut_macr_mus <- set_expt_conditions(neut_macr_mus, fact="infectstate")
neut_macr_mus <- set_expt_batches(neut_macr_mus, fact="hostcelltype")
neut_macr_mus_norm <- normalize_expt(neut_macr_mus, convert="cpm",
norm="quant", filter=TRUE)
## This function will replace the expt$expressionset slot with:
## cpm(quant(cbcb(data)))
## It will save copies of each step along the way
## in expt$normalized with the corresponding libsizes. Keep libsizes in mind
## when invoking limma. The appropriate libsize is non-log(cpm(normalized)).
## This is most likely kept at:
## 'new_expt$normalized$intermediate_counts$normalization$libsizes'
## A copy of this may also be found at:
## new_expt$best_libsize
## 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.
## Not correcting the count-data for batch effects. If batch is
## included in EdgerR/limma's model, then this is probably wise; but in extreme
## batch effects this is a good parameter to play with.
## Step 1: performing count filter with option: cbcb
## Removing 0 low-count genes (19660 remaining).
## Step 2: normalizing the data with quant.
## Step 3: converting the data with cpm.
## Step 4: not transforming the data.
## Step 5: not doing batch correction.
neut_macr_mus_normbatch <- normalize_expt(neut_macr_mus, convert="cpm",
norm="quant", filter=TRUE, batch="svaseq")
## This function will replace the expt$expressionset slot with:
## svaseq(cpm(quant(cbcb(data))))
## It will save copies of each step along the way
## in expt$normalized with the corresponding libsizes. Keep libsizes in mind
## when invoking limma. The appropriate libsize is non-log(cpm(normalized)).
## This is most likely kept at:
## 'new_expt$normalized$intermediate_counts$normalization$libsizes'
## A copy of this may also be found at:
## new_expt$best_libsize
## 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.
## Warning in normalize_expt(neut_macr_mus, convert = "cpm", norm = "quant", :
## Quantile normalization and sva do not always play well together.
## Step 1: performing count filter with option: cbcb
## Removing 0 low-count genes (19660 remaining).
## Step 2: normalizing the data with quant.
## Step 3: converting the data with cpm.
## Step 4: not transforming the data.
## Step 5: doing batch correction with svaseq.
## Note to self: If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.
## Passing off to all_adjusters.
## batch_counts: Before batch/surrogate estimation, 515332 entries are x>1: 47%.
## batch_counts: Before batch/surrogate estimation, 323375 entries are x==0: 29%.
## batch_counts: Before batch/surrogate estimation, 262253 entries are 0<x<1: 24%.
## The be method chose 10 surrogate variable(s).
## Attempting svaseq estimation with 10 surrogates.
## There are 30709 (3%) elements which are < 0 after batch correction.
## Setting low elements to zero.
## This function will replace the expt$expressionset slot with:
## simple(data)
## It will save copies of each step along the way
## in expt$normalized with the corresponding libsizes. Keep libsizes in mind
## when invoking limma. The appropriate libsize is non-log(cpm(normalized)).
## This is most likely kept at:
## 'new_expt$normalized$intermediate_counts$normalization$libsizes'
## A copy of this may also be found at:
## new_expt$best_libsize
## Leaving the data in its current base format, keep in mind that
## some metrics are easier to see when the data is log2 transformed, but
## EdgeR/DESeq do not accept transformed data.
## Leaving the data unconverted. It is often advisable to cpm/rpkm
## the data to normalize for sampling differences, keep in mind though that rpkm
## has some annoying biases, and voom() by default does a cpm (though hpgl_voom()
## will try to detect this).
## Leaving the data unnormalized. This is necessary for DESeq, but
## EdgeR/limma might benefit from normalization. Good choices include quantile,
## size-factor, tmm, etc.
## 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: simple
## Removing 3175 low-count genes (16485 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.
neut_macr_mus_de <- all_pairwise(neut_macr_mus_filt, parallel=FALSE,
force=TRUE, model_batch="svaseq")
## batch_counts: Before batch/surrogate estimation, 616059 entries are x>1: 67%.
## batch_counts: Before batch/surrogate estimation, 269846 entries are x==0: 29%.
## batch_counts: Before batch/surrogate estimation, 2347 entries are 0<x<1: 0%.
## The be method chose 3 surrogate variable(s).
## Attempting svaseq estimation with 3 surrogates.
## Plotting a PCA before surrogates/batch inclusion.
## Not putting labels on the plot.
## Using svaseq to visualize before/after batch inclusion.
## Performing a test normalization with: raw
## This function will replace the expt$expressionset slot with:
## log2(svaseq(cpm(cbcb(data))))
## It will save copies of each step along the way
## in expt$normalized with the corresponding libsizes. Keep libsizes in mind
## when invoking limma. The appropriate libsize is non-log(cpm(normalized)).
## This is most likely kept at:
## 'new_expt$normalized$intermediate_counts$normalization$libsizes'
## A copy of this may also be found at:
## new_expt$best_libsize
## Leaving the data unnormalized. This is necessary for DESeq, but
## EdgeR/limma might benefit from normalization. Good choices include quantile,
## size-factor, tmm, etc.
## Step 1: performing count filter with option: cbcb
## Removing 0 low-count genes (16485 remaining).
## Step 2: not normalizing the data.
## Step 3: converting the data with cpm.
## Step 4: transforming the data with log2.
## transform_counts: Found 269846 values equal to 0, adding 1 to the matrix.
## Step 5: doing batch correction with svaseq.
## Note to self: If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.
## Passing off to all_adjusters.
## batch_counts: Before batch/surrogate estimation, 511083 entries are x>1: 55%.
## batch_counts: Before batch/surrogate estimation, 269846 entries are x==0: 29%.
## batch_counts: Before batch/surrogate estimation, 142231 entries are 0<x<1: 15%.
## The be method chose 7 surrogate variable(s).
## Attempting svaseq estimation with 7 surrogates.
## There are 61691 (7%) elements which are < 0 after batch correction.
## Setting low elements to zero.
## Not putting labels on the plot.
## Starting basic_pairwise().
## Starting basic pairwise comparison.
## Leaving the data alone, regardless of normalization state.
## Basic step 0/3: Transforming data.
## Basic step 1/3: Creating mean and variance tables.
## Basic step 2/3: Performing 3 comparisons.
## Basic step 3/3: Creating faux DE Tables.
## Basic: Returning tables.
## Starting deseq_pairwise().
## Starting DESeq2 pairwise comparisons.
## About to round the data, this is a pretty terrible thing to do. But if you, like me, want to see what happens when you put non-standard data into deseq, then here you go.
## Warning in choose_binom_dataset(input, force = force): This data was
## inappropriately forced into integers.
## Including batch estimates from sva/ruv/pca in the model.
## Choosing the non-intercept containing model.
## DESeq2 step 1/5: Including a matrix of batch estimates in the deseq model.
## converting counts to integer mode
## DESeq2 step 2/5: Estimate size factors.
## DESeq2 step 3/5: Estimate dispersions.
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## Using a parametric fitting seems to have worked.
## DESeq2 step 4/5: nbinomWaldTest.
## Starting ebseq_pairwise().
## About to round the data, this is a pretty terrible thing to do. But if you, like me, want to see what happens when you put non-standard data into deseq, then here you go.
## Warning in choose_binom_dataset(input, force = force): This data was
## inappropriately forced into integers.
## Starting EBSeq pairwise subset.
## Choosing the non-intercept containing model.
## Starting EBTest of no vs. yes.
## Forcing out NA values by putting in the mean of all data.
## Copying ppee values as ajusted p-values until I figure out how to deal with them.
## Starting edger_pairwise().
## Starting edgeR pairwise comparisons.
## About to round the data, this is a pretty terrible thing to do. But if you, like me, want to see what happens when you put non-standard data into deseq, then here you go.
## Warning in choose_binom_dataset(input, force = force): This data was
## inappropriately forced into integers.
## Including batch estimates from sva/ruv/pca in the model.
## Choosing the non-intercept containing model.
## EdgeR step 1/9: Importing and normalizing data.
## EdgeR step 2/9: Estimating the common dispersion.
## EdgeR step 3/9: Estimating dispersion across genes.
## EdgeR step 4/9: Estimating GLM Common dispersion.
## EdgeR step 5/9: Estimating GLM Trended dispersion.
## EdgeR step 6/9: Estimating GLM Tagged dispersion.
## EdgeR step 7/9: Running glmFit, switch to glmQLFit by changing the argument 'edger_test'.
## EdgeR step 8/9: Making pairwise contrasts.
## Starting limma_pairwise().
## Starting limma pairwise comparison.
## Leaving the data alone, regardless of normalization state.
## libsize was not specified, this parameter has profound effects on limma's result.
## Using the libsize from expt$best_libsize.
## Limma step 1/6: choosing model.
## Including batch estimates from sva/ruv/pca in the model.
## Choosing the non-intercept containing model.
## Limma step 2/6: running limma::voom(), switch with the argument 'which_voom'.
## Using normalize.method=quantile for voom.
## Limma step 3/6: running lmFit with method: ls.
## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)
## Limma step 5/6: Running eBayes with robust=FALSE and trend=FALSE.
## Limma step 6/6: Writing limma outputs.
## Limma step 6/6: 1/1: Creating table: yes_vs_no. Adjust=BH
## Limma step 6/6: 1/2: Creating table: no. Adjust=BH
## Limma step 6/6: 2/2: Creating table: yes. Adjust=BH
## Comparing analyses.
## Writing a legend of columns.
## Working on 1/1: infection which is: yes/no.
## Found table with yes_vs_no
R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
locale: LC_CTYPE=en_US.UTF-8, LC_NUMERIC=C, LC_TIME=en_US.UTF-8, LC_COLLATE=en_US.UTF-8, LC_MONETARY=en_US.UTF-8, LC_MESSAGES=en_US.UTF-8, LC_PAPER=en_US.UTF-8, LC_NAME=C, LC_ADDRESS=C, LC_TELEPHONE=C, LC_MEASUREMENT=en_US.UTF-8 and LC_IDENTIFICATION=C
attached base packages: parallel, stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: edgeR(v.3.28.0), ruv(v.0.9.7.1), hpgltools(v.1.0), Biobase(v.2.46.0) and BiocGenerics(v.0.32.0)
loaded via a namespace (and not attached): rappdirs(v.0.3.1), rtracklayer(v.1.46.0), R.methodsS3(v.1.7.1), tidyr(v.1.0.2), ggplot2(v.3.2.1), acepack(v.1.4.1), bit64(v.0.9-7), knitr(v.1.28), R.utils(v.2.9.2), DelayedArray(v.0.12.2), data.table(v.1.12.8), rpart(v.4.1-15), RCurl(v.1.98-1.1), doParallel(v.1.0.15), GenomicFeatures(v.1.38.1), preprocessCore(v.1.48.0), callr(v.3.4.1), cowplot(v.1.0.0), usethis(v.1.5.1), RSQLite(v.2.2.0), europepmc(v.0.3), bit(v.1.1-15.2), enrichplot(v.1.6.1), xml2(v.1.2.2), httpuv(v.1.5.2), SummarizedExperiment(v.1.16.1), assertthat(v.0.2.1), viridis(v.0.5.1), xfun(v.0.12), tximport(v.1.14.0), hms(v.0.5.3), evaluate(v.0.14), promises(v.1.1.0), DEoptimR(v.1.0-8), fansi(v.0.4.1), progress(v.1.2.2), readxl(v.1.3.1), caTools(v.1.18.0), dbplyr(v.1.4.2), igraph(v.1.2.4.2), DBI(v.1.1.0), geneplotter(v.1.64.0), htmlwidgets(v.1.5.1), stats4(v.3.6.1), purrr(v.0.3.3), ellipsis(v.0.3.0), crosstalk(v.1.0.0), dplyr(v.0.8.4), backports(v.1.1.5), annotate(v.1.64.0), biomaRt(v.2.42.0), blockmodeling(v.0.3.4), vctrs(v.0.2.2), remotes(v.2.1.0), withr(v.2.1.2), ggforce(v.0.3.1), triebeard(v.0.3.0), robustbase(v.0.93-5), checkmate(v.2.0.0), GenomicAlignments(v.1.22.1), prettyunits(v.1.1.1), cluster(v.2.1.0), DOSE(v.3.12.0), lazyeval(v.0.2.2), crayon(v.1.3.4), genefilter(v.1.68.0), pkgconfig(v.2.0.3), labeling(v.0.3), tweenr(v.1.0.1), GenomeInfoDb(v.1.22.0), nlme(v.3.1-144), pkgload(v.1.0.2), nnet(v.7.3-12), devtools(v.2.2.1), rlang(v.0.4.4), lifecycle(v.0.1.0), BiocFileCache(v.1.10.2), directlabels(v.2020.1.31), cellranger(v.1.1.0), rprojroot(v.1.3-2), polyclip(v.1.10-0), matrixStats(v.0.55.0), rngtools(v.1.5), graph(v.1.64.0), Matrix(v.1.2-18), urltools(v.1.7.3), boot(v.1.3-24), base64enc(v.0.1-3), ggridges(v.0.5.2), processx(v.3.4.2), png(v.0.1-7), viridisLite(v.0.3.0), bitops(v.1.0-6), R.oo(v.1.23.0), KernSmooth(v.2.23-16), pander(v.0.6.3), EBSeq(v.1.26.0), Biostrings(v.2.54.0), blob(v.1.2.1), doRNG(v.1.8.2), stringr(v.1.4.0), qvalue(v.2.18.0), gProfileR(v.0.7.0), readr(v.1.3.1), jpeg(v.0.1-8.1), gridGraphics(v.0.4-1), S4Vectors(v.0.24.3), scales(v.1.1.0), memoise(v.1.1.0), magrittr(v.1.5), plyr(v.1.8.5), gplots(v.3.0.1.2), gdata(v.2.18.0), zlibbioc(v.1.32.0), compiler(v.3.6.1), RColorBrewer(v.1.1-2), lme4(v.1.1-21), DESeq2(v.1.26.0), Rsamtools(v.2.2.1), cli(v.2.0.1), XVector(v.0.26.0), ps(v.1.3.0), htmlTable(v.1.13.3), Formula(v.1.2-3), MASS(v.7.3-51.5), mgcv(v.1.8-31), tidyselect(v.1.0.0), stringi(v.1.4.5), yaml(v.2.2.1), GOSemSim(v.2.12.0), askpass(v.1.1), locfit(v.1.5-9.1), latticeExtra(v.0.6-29), ggrepel(v.0.8.1), grid(v.3.6.1), fastmatch(v.1.1-0), tools(v.3.6.1), rstudioapi(v.0.11), foreach(v.1.4.8), foreign(v.0.8-75), gridExtra(v.2.3), farver(v.2.0.3), Rtsne(v.0.15), ggraph(v.2.0.1), digest(v.0.6.23), rvcheck(v.0.1.7), BiocManager(v.1.30.10), shiny(v.1.4.0), quadprog(v.1.5-8), Rcpp(v.1.0.3), GenomicRanges(v.1.38.0), later(v.1.0.0), httr(v.1.4.1), AnnotationDbi(v.1.48.0), colorspace(v.1.4-1), XML(v.3.99-0.3), fs(v.1.3.1), IRanges(v.2.20.2), splines(v.3.6.1), RBGL(v.1.62.1), graphlayouts(v.0.5.0), ggplotify(v.0.0.4), plotly(v.4.9.1), sessioninfo(v.1.1.1), xtable(v.1.8-4), jsonlite(v.1.6.1), nloptr(v.1.2.1), tidygraph(v.1.1.2), corpcor(v.1.6.9), testthat(v.2.3.1), Vennerable(v.3.1.0.9000), R6(v.2.4.1), Hmisc(v.4.3-1), pillar(v.1.4.3), htmltools(v.0.4.0), mime(v.0.9), glue(v.1.3.1), fastmap(v.1.0.1), minqa(v.1.2.4), clusterProfiler(v.3.14.3), BiocParallel(v.1.20.1), codetools(v.0.2-16), fgsea(v.1.12.0), pkgbuild(v.1.0.6), lattice(v.0.20-38), tibble(v.2.1.3), sva(v.3.34.0), pbkrtest(v.0.4-7), curl(v.4.3), colorRamps(v.2.3), gtools(v.3.8.1), zip(v.2.0.4), GO.db(v.3.10.0), openxlsx(v.4.1.4), openssl(v.1.4.1), survival(v.3.1-8), limma(v.3.42.2), rmarkdown(v.2.1), desc(v.1.2.0), munsell(v.0.5.0), DO.db(v.2.9), fastcluster(v.1.1.25), GenomeInfoDbData(v.1.2.2), iterators(v.1.0.12), variancePartition(v.1.16.1), reshape2(v.1.4.3) and gtable(v.0.3.0)
## If you wish to reproduce this exact build of hpgltools, invoke the following:
## > git clone http://github.com/abelew/hpgltools.git
## > git reset c484b08e559f3bc0530c1fc765f6b05ba584bb80
## This is hpgltools commit: Wed Feb 12 15:10:29 2020 -0500: c484b08e559f3bc0530c1fc765f6b05ba584bb80