3 How to use TcGSA for one treatment group ?
3.1 Data preparation
This example comes from Hejblum et al., 2015 and DALIA-1 HIV therapeutic vaccine trial. The aim of this study is to evaluate the immune response to HIV vaccine. To conduct this study, 19 patients contaminated by the HIV have been followed for 48 weeks split into 2 phases of 24 each:
- 1st phase from week 0 to week 24 (before treatment interruption): the patients were under anti-retroviral treatment and received trial vaccine on week 0, 4, 8 and 12.
- 2nd phase from week 24 to week 48 (after treatment interruption): for the follow-up of patients, none of trial vaccine has been injected and anti-retroviral treatments have been interrupted (except for health problems).
Blood samples have been collected at each of the different measurement time points, for each subject, to study the dynamic of gene expression over time. For more details, check the article from Hejblum et al. here.
The data are publicly available on GEO website with GEO access number ‘GSE46734’. We will be using the GEOquery package to get the data files from GEO website (see appendix for more details on GEOquery)
3.1.1 Import of data files
In this example, we need import the supplementary files available on GEOwith the getGEOSuppFiles function (we only need the three following files, hence the regular expression filter: gene expression pre-ATI, gene expression post-ATI, experimental design):
GEOquery::getGEOSuppFiles('GSE46734', filter_regex="(*NonParamCombat*)|(*DESIGN*)")3.1.2 Design data matrix
The design data matrix (called design_preATI) is extracted from one of the GEO supplementary files. It contains the needed experimental variables, plus some additional information regarding this study. Some data processing is performed made according to the source of this paper.
design <- read.delim(gzfile("GSE46734/GSE46734_DALIA1longitudinalTranscriptome_DESIGN_anonym.txt.gz"))
design_preATI <- design[-which(design$TimePoint<0 | design$TimePoint==16 | design$TimePoint>22), ]
head(design_preATI,5)
#> Sample_name Chip_ID HYB_Chamber HYB_Day Patient_ID TimePoint
#> 1 X6116733008_L C4 CH1 D1 P001 4
#> 3 X6116733022_E C13 CH4 D1 P001 22
#> 5 X6116733029_A C18 CH1 D2 P001 8
#> 6 X6116733068_L C25 CH3 D2 P001 12
#> 8 X6116733073_K C29 CH4 D2 P001 0This data frame contains 90 samples and 6 experimental variables :
Sample_namefor the name of samplesPatient_IDfor the identification of patientsTimePointfor the time measurementsChip_IDHYB_ChamberHYB_Dayare the variables not required for TcGSA commands
3.1.3 Gene expression matrix
The gene expression matrix (called expr_preATI) is extracted from one of the GEO supplementary files (namely the “GSE46734_DALIA1longitudinalTranscriptome_PALO01_PreATI_NEQC_NonParamCombat.txt.gz” file).
NB: The data is already normalized.
expr_preATI <- read.delim(gzfile("GSE46734/GSE46734_DALIA1longitudinalTranscriptome_PALO01_PreATI_NEQC_NonParamCombat.txt.gz"))
rownames(expr_preATI) <- expr_preATI$PROBE_ID
expr_preATI <- expr_preATI[,as.character(design_preATI$Sample_name)]
expr_preATI[1:4,1:4]
#> X6116733008_L X6116733022_E X6116733029_A X6116733068_L
#> ILMN_1762337 4.967101 4.398911 4.384093 4.927489
#> ILMN_2055271 5.215882 4.744463 5.038972 5.209726
#> ILMN_2383229 4.565933 4.578123 4.696965 4.306206
#> ILMN_1806310 4.402265 4.473504 4.598726 4.589389We have:
- in row
ILMN_xxxxxxxfor each probe identifier - in column
Xxxxxxxxxxx_Xfor the name of each sample - each cell contains the normalized gene expression
The entire matrix contains \(32,978\) genes and \(90\) samples (number of samples should be the same as in design data matrix)
identical(ncol(expr_preATI), nrow(design_preATI))
#> [1] TRUE3.2 Likelihood ratios test
This function provides the result of likelihood ratio test using the linear mixed model for each gene set. For this example, we use gene sets data from Chaussabel’s modules (Chaussabel et al., 2008). TcGSA.LR function requires:
expr: name of the gene expression matrixgmt: name of the gmt gene set objectdesign: name of the design data matrixsubject_name: name of the identification of patients in the design data matrixtime_name: name of the time measurements in the design data matrix
tcgsa_result <- TcGSA::TcGSA.LR(expr = expr_preATI,
gmt = gmt_modulesV2,
design = design_preATI,
subject_name = "Patient_ID",
time_name = "TimePoint")#> A TcGSA object
#> Form of the time trend:
#> linear
#> Number of treatment groups:
#> 1
#> Number of gene sets tested for significant time trend:
#> 260Now tcgsa_result is a tcgsa object containing, in addition to the likelihood ratio test results:
- the form of tested time trends (default option is
linear) - the number of treatment groups and the number of gene sets tested (depending on which gene set base is defined)
- a few other information
To get the number of significant gene sets, one can use summary function on a tcgsa object:
summary(tcgsa_result)
#> A TcGSA object
#> Form of the time trend:
#> linear
#> Number of treatment groups:
#> 1
#> Number of gene sets tested for significant time trend:
#> 260
#>
#> Number of significant gene sets at a 5% threshold (BY procedure):
#> 17 out of 260 gene setsTo get more details on the significant gene sets, use the signifLRT.TcGSA() function. It returns information such as the significant gene sets among all the gene sets tested, along their p-values with adjustment for multiple testing (default option is BY for Benjamini-Yekutieli correction Benjamini et Yekutieli, 2001 and 5% threshold). Below is an example of five significant gene sets:
head(TcGSA::signifLRT.TcGSA(tcgsa_result)$mixedLRTadjRes)
#> GeneSet AdjPval desc
#> 1 M1.1 8.296095e-09 Platelets
#> 2 M3.2 4.096412e-08 Inflammation
#> 3 M3.5 3.952334e-02 Cell Cycle
#> 6 M4.6 1.337921e-06 Inflammation
#> 7 M4.9 3.820230e-04 Undetermined
#> 4 M4.11 7.411800e-07 Plasma CellsYou can also use the multtest.TcGSA function to provide the likelihood ratios, the raw and adjusted p-values for the whole gene sets with 5% threshold. Below is an example displaying only for five results:
head(TcGSA::multtest.TcGSA(tcgsa_result))
#> LR CVG_H0 CVG_H1 raw_pval adj_pval
#> 1 48.5951385 0 0 1.559071e-11 8.296095e-09
#> 2 0.4574752 0 0 6.471718e-01 1.000000e+00
#> 3 3.1626395 0 0 1.405224e-01 1.000000e+00
#> 4 1.6726757 0 0 3.145975e-01 1.000000e+00
#> 5 0.1412185 0 0 8.194493e-01 1.000000e+00
#> 6 0.1036336 0 0 8.485074e-01 1.000000e+00CVG_H0 and CVG_H1 are the convergence of the model under null and alternative hypotheses. 0 indicates a good convergence of the model (based on lme4 output).
3.3 Graphical outputs
plot1GS() plots the different representations of gene expression in a specific gene set of interest (specified by the geneset.name argument). This function requires the following:
expr: either the name of the gene expression matrix or the estimations of linear mixed model (in this example, we used the raw data from the gene expression matrix)gmt: the name of the gmt gene set objectSubject_ID: the name of the identification of patients in the design data matrixTimePoint: the name of the time measurements of the design data matrixgeneset.name: the name of gene set (significant ones can be found withsignifLRT.TcGSA(tcgsa_result)$mixedLRTadjRes)time_unit: string to be displayed before to the values ofTimePointon the x-axis (such as ‘D’ for ‘days’ for instance - optional)
TcGSA::plot1GS(expr = expr_preATI,
#plot1GS(expr = tcgsa_result$Estimations,
gmt = gmt_modulesV2,
Subject_ID = design_preATI$Patient_ID,
TimePoint = design_preATI$TimePoint,
clustering = FALSE,
time_unit = "W",
geneset.name = "M3.2",
title="",
margins=0.4,
lab.cex=0.37,
axis.cex=0.37,
line.size=0.45,
gg.add=list(ggplot2::theme(legend.position="none"),
ggplot2::ylim(-1.26,1.26)
))
Dotted line shows the median gene expression across subjects, in the gene set over time.
