gse <- getGEO("GSE5220", destdir = "./")
## Found 1 file(s)
## GSE5220_series_matrix.txt.gz
## Using locally cached version: .//GSE5220_series_matrix.txt.gz
## Using locally cached version of GPL96 found here:
## .//GPL96.soft.gz
gpl <- getGEO("GPL96", destdir = "./")
## Using locally cached version of GPL96 found here:
## .//GPL96.soft.gz
 eset <- gse[[1]]

Extract Data

gene_counts <- exprs(eset)
gene_identifiers <- fData(eset)
phenotype_data <- pData(eset)
Ensure Data is formatted correctly
gene.entrez_df <- gene_identifiers[, c("Gene Symbol", "ENTREZ_GENE_ID")]
gene.entrez_df$ENTREZ_GENE_ID <- sapply(strsplit(as.character(gene.entrez_df$ENTREZ_GENE_ID), " /// "), `[`, 1)
gene.entrez_df$'Gene Symbol' <- sapply(strsplit(as.character(gene.entrez_df$'Gene Symbol'), " /// "), `[`, 1)

phenotype_data$group <- ifelse(grepl("aviremic", phenotype_data$source_name_ch1, ignore.case = TRUE), "aviremic", "viremic")

# Convert this new column to a factor
phenotype_data$group <- as.factor(phenotype_data$group)
Filter and Align Data
keep <- rowSums(cpm(gene_counts) > 1) >= 2  # Or another appropriate threshold
gene_counts <- gene_counts[keep,]
gene.entrez_df <- gene.entrez_df[keep,]

# Make sure all data frames/objects are aligned
if (!all(rownames(gene_counts) == rownames(gene.entrez_df))) {
    stop("Gene identifiers between count data and gene annotations do not match.")
}

# Create a DGEList object for edgeR analysis
DGE <- DGEList(counts=gene_counts, samples = phenotype_data, genes = gene.entrez_df)
Normalise Data
# Normalise data with Z-scores
z_scores <- t(scale(t(gene_counts)))  # Transpose twice to scale rows (genes)

# Update the gene_counts object
gene_counts <- z_scores

Differential Expression Analysis

# Assuming you have a factor in your phenotype data that defines the groups
design <- model.matrix(~ group, data = phenotype_data)
# Perform the robust analysis using voom with quality weights and then lmFit
v <- voom(DGE, design, plot=TRUE)

fit <- lmFit(v, design)
efit <- eBayes(fit, robust=TRUE)  # Ensure robustness in t-statistics

  • X-axis (Log2 Count Size + 0.5): Shows the average expression level of genes across samples.

  • Y-axis (Standard Deviation): Standard deviation of the log-transformed counts. It measures the variability of gene expression across samples.

  • The plot shows lower variability in gene expression at medium expression levels, as shown by the curve dipping in the middle. This suggests that genes with moderate expression levels tend to have more consistent expression across different samples.

  • There is an increase in variability at very low and very high expression levels. At low levels, this could be due to noise and measurement errors (more pronounced when dealing with small counts). At high levels, the increased variability might be due to greater impact of transcriptional regulation in highly expressed genes.

Differentially expressed genes
# Get results
results <- decideTests(fit, adjust.method="BH", p.value=0.05)

# Differentially expressed genes 
dt <- topTable(efit, adjust.method="BH", sort.by="P", number=Inf)
## Removing intercept from test coefficients
dt <- na.omit(dt)
sum(is.na(dt))
## [1] 0
Volcano plot for differentially expressed genes
# Rename for ggplot
dt$Log2FoldChange <- dt$logFC  
dt$NegLog10PValue <- -log10(dt$P.Value)

# Highlight significant points
significant <- dt[dt$P.Value < 0.05, ]  # Adjust the threshold as needed
dt$significant <- ifelse(dt$P.Value < 0.05, "Significant", "Not Significant")

ggplot(dt, aes(x = Log2FoldChange, y = NegLog10PValue, colour = significant)) +
  geom_point(alpha = 0.4) +  
  scale_color_manual(values = c("black", "red"), labels = c("Not Significant, P.value > 0.05", "Significant, P.value < 0.05")) +
  theme_minimal() +
  labs(title = "Volcano Plot",
       x = "Log2 Fold Change",
       y = "-Log10 P-value", 
       colour = "Significance Level ") +
  geom_hline(yintercept = -log10(0.05), col = "red", linetype = "dashed")

Upregulated genes

upregulated <- sum(significant$logFC > 1.5)
cat("Upregulated genes:", upregulated, "\n")
## Upregulated genes: 181
upregulated_genes <- significant[significant$logFC > 1.5, ]
Top 20 upregulated genes
# Sort by logFC in descending order to get the top upregulated genes
top_upregulated_genes <- significant[order(significant$logFC, decreasing = TRUE), ][1:20, ]
top_upregulated_genes$Gene.Symbol
##  [1] "IFI44L"  "SIGLEC1" "IFI27"   "NAB1"    "TOX"     "DNASE2B" "DBR1"   
##  [8] "GHR"     "IGK"     "DET1"    "HSPA4L"  "SULT1C2" "KAL1"    "ACTR5"  
## [15] "MMP8"    "CEP290"  "VGF"     "FAM182A" "RUFY3"   "SORCS3"

Downregulated genes

downregulated <- sum(significant$logFC < -1.5)
cat("Downregulated genes:", downregulated, "\n")
## Downregulated genes: 59
downregulated_genes <- significant[significant$logFC < -1.5, ]
Top 20 downregulated genes
# Sort by logFC in ascending order to get the top downregulated genes
top_downregulated_genes <- significant[order(significant$logFC, decreasing = FALSE), ][1:20, ]
top_downregulated_genes$Gene.Symbol
##  [1] "LOC100996756" "SYNJ1"        "CD1C"         "TRIM58"       "MINPP1"      
##  [6] "EPB41"        "ANO1"         "DPEP2"        "TSR1"         "FBXL5"       
## [11] "PATZ1"        "ZNF764"       "PEX7"         "KCNJ4"        "SLMO2"       
## [16] "BCR"          "ICT1"         "FAM3C"        "ZCCHC11"      "VDR"

Up-regulated genes network

g_up <- graph_from_data_frame(d = data.frame(source = "HIV", target = top_upregulated_genes$Gene.Symbol), directed = FALSE)
plot(g_up, vertex.size=45, vertex.label.cex=0.5, vertex.color="white", edge.arrow.size=.1)

Down-regulated genes network

g_down <- graph_from_data_frame(d = data.frame(source = "HIV", target = top_downregulated_genes$Gene.Symbol), directed = FALSE)
plot(g_down, vertex.size=45, vertex.label.cex=0.5, vertex.color="white", edge.arrow.size=0.1)

Hierarchical clustering of up-regulated genes

upregulated_expr_matrix <- gene_counts[rownames(gene_counts) %in% rownames(top_upregulated_genes), ]
rownames(upregulated_expr_matrix) <- top_upregulated_genes$Gene.Symbol

pheatmap(upregulated_expr_matrix, 
         scale = "row",  # Normalizes the data by rows (genes)
         clustering_distance_rows = "euclidean",  # Distance metric for rows
         clustering_distance_cols = "euclidean",  # Distance metric for columns
         clustering_method = "complete",  # Clustering method
         show_rownames = TRUE,  # Shows row names (gene symbols)
         show_colnames = TRUE,  # Shows column names (sample IDs)
         color = colorRampPalette(c("blue", "white", "red"))(255),  # Color scale
         fontsize_row = 6,
         cellheight = ,# Adjust font size for row labels if necessary
         fontsize_col = 6)  # Adjust font size for column labels if necessary

Gene symbol on the y axis and samples on the x axis.

Four main clusters of gene symbols:

  • Group 1 (IF127,MMP8, SORCS3,IFI44L, FAM182A).

  • Group 2 (SIGLEC1, RUFY3, GHR, HSPA4L, DNASE2B, ACTR5).

  • Group 3 (DET1, DBR1, TOX, CEP290, NAB1).

  • Group 4 (SULT1C2, KAL1, IGK, VGF).

Three main clusters of samples:

  • Group 1 (GSM118290, GSM118292, GSM118293, GSM118295, GSM118291, GSM118294, GSM118296, GSM118297).

  • Group 2 (GSM118301, GSM118305).

  • Group 3 (GSM118300, GSM118299, GSM118302, GSM118304, GSM118298, GSM118303).

Hierarchical clustering of down-regulated genes

downregulated_expr_matrix <- gene_counts[rownames(gene_counts) %in% rownames(top_downregulated_genes), ]
rownames(downregulated_expr_matrix) <- top_downregulated_genes$Gene.Symbol

pheatmap(downregulated_expr_matrix, 
         scale = "row",  # Normalizes the data by rows (genes)
         clustering_distance_rows = "euclidean",  # Distance metric for rows
         clustering_distance_cols = "euclidean",  # Distance metric for columns
         clustering_method = "complete",  # Clustering method
         show_rownames = TRUE,  # Shows row names (gene symbols)
         show_colnames = TRUE,  # Shows column names (sample IDs)
         color = colorRampPalette(c("blue", "white", "red"))(255),  # Color scale
         fontsize_row = 6,
         cellheight = ,# Adjust font size for row labels if necessary
         fontsize_col = 6)  # Adjust font size for column labels if necessary

Gene symbol on the y-axis and samples on the x-axis.

Four main clusters of gene symbols:

  • Group 1 (ZNF764, PATZ1, TRIM58, BCR, PEX7, LOC100996756).

  • Group 2 (ICT1, EPB41, SYNJ1, TSR1, ANO1).

  • Group 3 (ZCCHC11, FBXL5, DPEP2).

  • Group 4 (SLMO2, KCNJ4, FAM3C, CD1C, VDR, MIPP1).

Four main clusters of samples:

  • Group 1 (GSM118293, GSM118295).

  • Group 2 (GSM118290).

  • Group 3 (GSM118299, GSM118300, GSM118298, GSM118302, GSM118303, GSM118304).

  • Group 4 (GSM118292, GSM118305, GSM118301, GSM118297, GSM118296, GSM118291, GSM118294).

Protien-Protein Interactions with STRINGdb

library(STRINGdb)
## Warning: package 'STRINGdb' was built under R version 4.3.3
high_string_db <- STRINGdb$new(version = "11.0b", species = 9606, score_threshold = 800)
class(high_string_db)
## [1] "STRINGdb"
## attr(,"package")
## [1] "STRINGdb"
Protein-Protein Interactions (threshold 800) for Upregulated Genes
high_confidence_up <- high_string_db$map(upregulated_genes, "Gene.Symbol", removeUnmappedRows = TRUE)
## Warning:  we couldn't map to STRING 7% of your identifiers
hits3 <- high_confidence_up$STRING_id
high_string_db$plot_network(hits3)

Protein-Protein Interactions (threshold 800) for Down-regulated Genes

high_confidence_down <- high_string_db$map(downregulated_genes, "Gene.Symbol", removeUnmappedRows = TRUE)
## Warning:  we couldn't map to STRING 8% of your identifiers
hits4 <- high_confidence_down$STRING_id
high_string_db$plot_network(hits4)

Clusters of Up-regulated Genes
upclustersList <- high_string_db$get_clusters(high_confidence_up$STRING_id)
for(i in seq(4)){
   high_string_db$plot_network(upclustersList[[i]])
   }

  • Protein Cluster 1 (11 interactions): TRIM31, GBP1, IFI27, ISG15, IFI44L, IFI6

  • Protein Cluster 2 (5 interactions): CYFIP2, IGLL5, PRKCA, ITPR2, ANK2, SPTA1

  • Protein Cluster 3 (13 interactions): VGF, ALB, TEK, SERPING1, FLT1, HGF, MMRN1

  • Protein Cluster 4 (6 interactions): RPAP2, SNAPC5, SP1, CDK9

Clusters of Down-regulated Genes
downclustersList <- high_string_db$get_clusters(high_confidence_down$STRING_id)
for(i in seq(3)){
   high_string_db$plot_network(downclustersList[[i]])
   }

  • Protein Cluster 1 (3 interactions): ANO1, CFTR, SYNJI, MINPP1

  • Protein Cluster 2 (4 interactions): RNF111, MYLIP, FBXL5, CCT6A

  • Protein Cluster 3 (1 interactions): PRDX6, FAF2