This analysis examines microarray samples from two different breast cancer studies, and then compares the expression of differentiation between normal and malignant cells. All the samples used in this analysis were obtained from NCBI GO.

Out of the 10 microarray sample data, six (6) were obtained from the accession code “GSE19697” and the remaining four (4) which included the two samples used as control in this analysis from “GSE1045192”. A GSE19697_family.soft file was also obtained to collate probe and gene information from the study. ComBat, a function from the sva bioconductor package was used in the analysis of the two different study samples to eliminate experimental “noise” whilst still retaining all biologically relevant information.

After obtaining samples and and saving the .CEL files along with the .soft file into a working folder, analysis commenced.

For any more information, contact me on linked-in linkedin profile

WORKFLOW IN R

Routine Analytical Steps

  1. Set the working directory in R
setwd("C:/Users/ALFRED/Documents/Bioinformatics_with_R")
  1. Calling the important packages using the library function.
  • Oligo and limma perform almost the same function, but there are a few unique commands in each
library(affy)
library(oligo)
library(limma)
library(tidyverse)
library(sva)
  1. Creating a metadata table for the sample files consisting of filenames, sample and type
filename<- c("GSM1045191_Normal.CEL", "GSM1045192_Normal.CEL", "GSM1045274_Cancer.CEL", "GSM1045267_Cancer.CEL", "GSM1045234_Cancer.CEL",
             "GSM491610.CEL", "GSM491611.CEL", "GSM491612.CEL", "GSM491613.CEL", "GSM491614.CEL")
sample<- c("GSE1045192", "GSE1045192","GSE1045192", "GSE1045192", "GSE1045192",
           "GSE19697", "GSE19697", "GSE19697", "GSE19697", "GSE19697")
type<- c("CNT", "CNT", "cas", "cas", "cas", "cas", "cas", "cas", "cas", "cas")


metdata<- data.frame(filename, sample, type)
  1. Saving a backup of the generated dataframe which will be used later in the analysis
write.csv(metdata, file = "metdata.csv")
  1. Reading the .CEL files in the working directory using the created .csv file as a template
target<- readTargets("metdata.csv", sep = ",", row.names = "sample") 
affyc<- ReadAffy()
  • visualizing the affy data to show the variance between the probsets
hist(affyc)
  1. Normalize and correct the affy microarray with rma (Robust Multi-array Avarage)
  • This the standard program to remove background “noise” from the relevant biologically important prob signals. After normalization and correction, the expression information is read and saved as a .text file
eset<- affy::rma(affyc) 
expres_eset<- exprs(eset)
write.exprs(eset, file = "Expressed_eset.txt")
  • visualize the eset data to show the corrected variance between the prob sets
hist(eset)

Running Batch correction for the two data samples (GSE19697 & GSE1045192)

  • Running rma alone is not enough because data is being sourced from two different microarray flow-cells. There is higher experimental “noise” to be corrected.
  • The GSE1045192 data samples contain the normal cell sample used in this analysis but in order for them to be normalized alongside the cancer cells from GSE19697, two cancer cell samples were added from GSE1045192. This still meets the 4:1 test-to-reference reference ratio needed to avoid heavy reliance on one reference data sample.
batch<- as.factor(target$sample)
condition<- as.factor(target$type)
mod<- model.matrix(~factor(condition))
corr_data<- ComBat(dat = expres_eset, batch = batch, mod = mod, par.prior = TRUE)
  1. Annotating the prob IDs against their respective genes using the .soft file.

  • All description information in the .soft file was removed leaving behind the table with its relevant information.

    1. Reading and assigning the file to a variable
reference<- read.delim("GSE19697_family.soft", check.names = FALSE)
ref<- data.frame(reference)
    1. Match the prob ID to the gene name by combining “expres_eset” & “ref” by the id column in “ref” and identifier tab in “expres_eset”.
  • convert the “expres_eset” to a .csv and add the column name ‘ID’ to the prob_info list column and save
write.csv(expres_eset, "Expressed_eset.csv")
probe_info<- read.delim("Expressed_eset.csv", sep = ",", check.names = FALSE)
  1. Create a new object “refer” with only the relevant information to be used for annotating the genes by splitting the Gene & ENTREZ IDs of “ref” into main and alternatives before selecting important parts (main)
refer<- ref %>% 
  separate(Gene.Symbol,into = c("Gene.Symbol", "Gene_alt"), sep = "/") %>% 
  separate(ENTREZ_GENE_ID,into = c("ENTREZ_GENE_ID", "ENTREZ_alt"), sep = "/") %>% 
  select(ID, GB_ACC, Gene.Symbol, ENTREZ_GENE_ID)
  1. Merge “prob_info” and “refer” by ID to combine annotations acting as the go-to reference for all the genes.
annot<- left_join(probe_info, refer, by= 'ID', copy = TRUE)
write.csv(annot, "annotation.csv")

Differential Gene Expression Analysis

  1. Create a design/model for the analysis
  • The model.matrix creates a dummy variable consisting of 0s & 1s. Since the “condition” object is a factor composed of two levels (i.e: “CNT” & “cas”), the model.matrix assigns the dummy variables based on that
condition<- as.factor(target$type) 
design<- model.matrix(~0+condition) 
colnames(design)= levels(case)
  1. Running the differential analysis on our the corr_data using the “design” model

  • Using “corr_data” instead of “eset” in running the analysis because a re-correction and normalization was done to properly combine the two datasets (“GSE1045192” & “GSE19697”).

  • fdr (false discovery rate) tightens the p-value %, allowing for lower false positives.

    fit<- lmFit(corr_data, design = design)
    fit<- eBayes(fit) # eBayes generates the statistical values of importance
    options(digits = 2)
    top<- topTable(fit, coef = 2, n= Inf, adjust='fdr')
  1. Adding the analysis to the annotation
data<- merge(top, annot, by.x= 0, by.y = "ID", all= TRUE)

# Going to select the order in-which I want the variable to appear
data<- data %>% 
  mutate(Gene_ID= data$Row.names) %>% 
  select(Gene_ID, Gene.Symbol, GB_ACC, ENTREZ_GENE_ID, logFC, AveExpr, t, P.Value, adj.P.Val, B, everything() )
  1. Backing-up the data as a .CSV file
write.table(data, file = "DGE_Cancer_Final.csv", row.names = TRUE, col.names = TRUE, sep = ",")