Phenotype information

Phenotype information is central to all aspects of the analysis. In the broad sense, it encompasses sample-level information beyond the experimental design (e.g., files of read counts, quality control metrics) to coordinate statistical analyses and data visualisation throughout downstream analyses. Defined in the earliest steps of an analytical workflow, it ensures traceability of sample data from files of raw data to every output and result.

Let us import experimental information:

pheno <- xlsx::read.xlsx("expdata/samples.xlsx", sheetName = "pheno")
colnames(pheno)

## [1] "File"      "Infection" "Status"    "Time"      "Lane"      "Plate"    
## [7] "Well"

To facilitate later use, let us retype some fields, reorder some factor levels, and add a set of additional fields:

• Sample: contains a shortened unique identifier derived from the count file name for each sample,
• Treatment: contains a factor that combines Infection and Status phenotypes
• Group: contains a factor that combines all experimental phenotypes (i.e., Time, Infection and Status)

pheno$File <- as.character(pheno$File)
pheno$Infection <-
  factor(pheno$Infection, c("Mock","STM-LT2","STM-D23580","Blank"))
pheno$Status <-
  factor(pheno$Status, c("Uninfected","Exposed","Infected","Blank","Bulk"))
pheno$Lane <- as.factor(pheno$Lane)
pheno$Sample <- gsub("WTCHG_[[:digit:]]+_", "Cell_", pheno$File) # character
treatmentLevels <- c(
  "Mock_Uninfected",
  "STM-LT2_Exposed","STM-LT2_Infected",
  "STM-D23580_Exposed","STM-D23580_Infected",
  "Blank_Blank", "Mock_Bulk", "STM-LT2_Bulk", "STM-D23580_Bulk"
  )
pheno$Treatment <- with(pheno, factor(
  droplevels(interaction(Infection, Status, sep = "_")), treatmentLevels
))
# factor
pheno$Group <- with(pheno, factor(
  droplevels(interaction(pheno$Time, pheno$Treatment, sep = "_"))
))
# factor

QC metrics

Importantly, MultiQC was used to collate metrics for each sample at each step of preprocessing—from raw reads to assigned read counts—; those QC metrics provide key information to identify technical issues and sample outliers.

Let us import the QC metrics produced by MultiQC:

multiqc <- read.delim("MultiQC/final_pipeline_data/multiqc_general_stats.txt")
colnames(multiqc)

## [1] "Sample"                          "featureCounts_Assigned"         
## [3] "FastQC_percent_gc"               "Bowtie.2_overall_alignment_rate"
## [5] "featureCounts_percent_assigned"  "FastQC_total_sequences"         
## [7] "FastQC_percent_duplicates"       "FastQC_percent_fails"           
## [9] "FastQC_avg_sequence_length"

Let us shorten the sample identifier used by MultiQC to a minimal unique identifier (as done above for experimental phenotype information):

multiqc$Sample <- gsub("WTCHG_[[:digit:]]+_", "Cell_", multiqc$Sample)
rownames(multiqc) <- multiqc$Sample
colnames(multiqc)

## [1] "Sample"                          "featureCounts_Assigned"         
## [3] "FastQC_percent_gc"               "Bowtie.2_overall_alignment_rate"
## [5] "featureCounts_percent_assigned"  "FastQC_total_sequences"         
## [7] "FastQC_percent_duplicates"       "FastQC_percent_fails"           
## [9] "FastQC_avg_sequence_length"

Let us append those QC metrics to the phenotype information, that will later be attached to the read counts in a SingleCellExperiment (SingleCellExperiment):

pheno <- merge(pheno, multiqc, by = "Sample", sort = FALSE)

Counts

To import count data that is stored in a separate file for each sample, let us use the edgeR readDGE method and the phenotype information prepared above:

RG <- with(
  pheno,
  readDGE(File, "expdata/counts", c(1, 3), Group, Sample)
)
dim(RG)

## [1] 58084   384