I will use this sheet for exporting data to collaborators. Currently, that is only the folks interested in doing some analyses with the first PLOS NTD paper.
It is worth noting that I recently found some errors in my accounting for these and the plosntd2 samples (I had the pathogen strains reversed).
hs_annot <- load_biomart_annotations()$annotation## 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)Note that I need to overwrite the savefile of the annotations and not use the archive server, since 38v99 came out in 202001.
sample_sheet <- glue::glue("sample_sheets/UMD_leishmania_host_metasheet_{ver}.xlsx")
## As of 20200222, I have only performed hisat2 mapping for the plos ntd2 data.
prefix <- "excel/plos_ntd_host_hisat2-v"
savefile <- glue::glue("mbio_expressionset-v{ver}.rda")
all_expt <- create_expt(metadata=sample_sheet,
file_column="hg3891salmon",
gene_info=new_hs_annot,
tx_gene_map=hs_tx_gene,
savefile=savefile)## Reading the sample metadata.## Dropped 2 rows from the sample metadata because they were blank.## The sample definitions comprises: 441 rows(samples) and 68 columns(metadata fields).## Reading count tables.## Using the transcript to gene mapping.## Reading salmon data with tximport.## Finished reading count data.## Warning in create_expt(metadata = sample_sheet, file_column = "hg3891salmon", :
## Some samples were removed when cross referencing the samples against the count
## data.## Matched 16933 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'.## The final expressionset has 16933 rows and 267 columns.mbio_expt <- subset_expt(all_expt, subset="study=='mbio'")## Using a subset expression.## There were 267, now there are 82 samples.mbio_expt <- set_expt_conditions(mbio_expt, fact="infectionstatus")
mbio_expt <- set_expt_batches(mbio_expt, fact="donor")
mbio_norm <- normalize_expt(mbio_expt, filter=TRUE, convert="cpm", norm="quant", transform="log2")## This function will replace the expt$expressionset slot with:## log2(cpm(quant(cbcb(data))))## It will save copies of each step along the way
## in expt$normalized with the corresponding libsizes. Keep libsizes in mind
## when invoking limma. The appropriate libsize is non-log(cpm(normalized)).
## This is most likely kept at:
## 'new_expt$normalized$intermediate_counts$normalization$libsizes'
## A copy of this may also be found at:
## new_expt$best_libsize## Not correcting the count-data for batch effects. If batch is
## included in EdgerR/limma's model, then this is probably wise; but in extreme
## batch effects this is a good parameter to play with.## Step 1: performing count filter with option: cbcb## Removing 5857 low-count genes (11076 remaining).## Step 2: normalizing the data with quant.## Step 3: converting the data with cpm.## The method is: raw.## Step 4: not doing batch correction.## Step 4: transforming the data with log2.## transform_counts: Found 59018 values equal to 0, adding 1 to the matrix.tt <- plot_pca(mbio_norm, plot_labels="normal")pander::pander(sessionInfo())
message(paste0("This is hpgltools commit: ", get_git_commit()))
this_save <- paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
message(paste0("Saving to ", this_save))
tmp <- sm(saveme(filename=this_save))