Part 1: Data Exploration and Visualization

Q1. Load the data and examine the structure of Y and patient_data. Provide a brief summary of what each object represents. Identify key features of the data, such as dimensions, column names, or unique values. (2 Marks)

Answer

Provide your answer below. When working in R Markdown, you can include code chunks to run your R code and display the results directly in the document. For example:

#load dataset
load("Coursework0525_Data.Rdata")

# Check which objects were loaded
ls()
## [1] "patient_data" "Y"
#object: Y
head(Y, 5) # see first 5 rows of data
##                 Sample001 Sample002 Sample003 Sample004 Sample005 Sample006
## ENSG00000157106      1555     29378         0     18618     62411     99033
## ENSG00000122565      3831     12505         0     15067      5170     23181
## ENSG00000124275       739      1643         0      2754       243     21602
## ENSG00000092036       583       784         0       505         0      1658
## ENSG00000160392        48       203         0       303       392      3527
##                 Sample007 Sample008 Sample009 Sample010 Sample011 Sample012
## ENSG00000157106      6605     69717      7174    119884     40832     17033
## ENSG00000122565     17583     13217     10208     73716     10901     15817
## ENSG00000124275      1413     14545       355      6099     11534       597
## ENSG00000092036       198      3513        50      5395         0       349
## ENSG00000160392        98      1140        19       674       681         0
##                 Sample013 Sample014 Sample015 Sample016 Sample017 Sample018
## ENSG00000157106     12519     18978         0    620483     22644     16391
## ENSG00000122565     26740     52647     34957     17486      1839      2815
## ENSG00000124275       905       524      4301         2     38044      7060
## ENSG00000092036       317        32      1203       268       181        44
## ENSG00000160392       289         0         0       416       815       230
##                 Sample019 Sample020 Sample021 Sample022 Sample023 Sample024
## ENSG00000157106     32627     35763      4542    138034      7703     13399
## ENSG00000122565      2676      1623      1350       896       599      1385
## ENSG00000124275         0     23599      4518     22728      2465      4968
## ENSG00000092036       145        88        41        14        20        53
## ENSG00000160392      1425      1391       442       705       523       269
##                 Sample025 Sample026 Sample027 Sample028 Sample029 Sample030
## ENSG00000157106     22736     22946      7874     77667      4893      4334
## ENSG00000122565      2032      5905      2227      2018      2245      1536
## ENSG00000124275      7694     21317      7347      6716     18238      2664
## ENSG00000092036        19       275        50       152        44        41
## ENSG00000160392       406         0       647       327       262         1
dim(Y)  # column names Sample001:030 with rows 1:4568 gene expression counts
## [1] 4568   30
class(Y) # matrix array 30 across and 4568 down.
## [1] "matrix" "array"
#object: patient_data
head(patient_data, 5) # see first 5 rows of data
##   tissue patient_id sample_id
## 1 Tumour  Patient.1 Sample001
## 2 Tumour  Patient.2 Sample002
## 3 Tumour  Patient.3 Sample003
## 4 Tumour  Patient.4 Sample004
## 5 Tumour  Patient.5 Sample005
dim(patient_data)  # 3 column named tissue, patient_id and sample_id with rows 30 patient data
## [1] 30  3
class(patient_data) # matrix array 3 across and 30 down.
## [1] "data.frame"
# Check unique values in key columns
unique(patient_data$tissue_type)
## NULL
unique(patient_data$patient_ID)
## NULL

You can run a code chunk by clicking the green “Run” button above it. This will execute the code and show the output in your R Markdown file. When you knit the document to HTML, the code chunks are run, and the output is included in the final HTML file.

Q2. Explain the concept of library size in the context of this dataset. Discuss potential reasons why library size may vary across samples. Plot library size across the sample size for this dataset. Using the provided code, normalize the data by library size. Mathematically describe the normalization process you applied. Finally, plot library size across the sample size for this dataset after normalization. Comment on how normalization affects library size distribution across samples. (2 Marks)

# Calculate library sizes (sum of counts per sample)
lib_sizes <- colSums(Y)

# Normalize counts by library size and scale to 1 million (CPM)
Y_cpm <- sweep(Y, 2, lib_sizes, FUN = "/") * 1e6

Answer

Q3. For normal and cancer samples of patient 5, provide summary statistics for expression values and visualize the distributions. Use the normalized dataset. Comment on skewness and multimodality. Transform the gene expression matrix using \(log2(Y + 1)\). For normal and cancer samples of patient 5, visualize the distributions and compare with the distributin in the previous task. Comment on skewness and multimodality. Now create a boxplot to show the distribution of gene expression across all samples before and after log transformation. You will have 30 boxplots side by side. Identify any problematic sample. Provide evidence for your selection. Remove the patient associated with the problematic sample. (4 Marks)

Answer

Part 2: Differential Expression Analysis

Q5. Study the following code and add comments to describe what each part does. Do not change the code. Derive the mathematical equation that represents the relationship between gene expression and the predictor variables based on the given Negative Binomial regression model. Justify the use of NB over Poisson by checking mean vs. variance or dispersion parameter on the dataset after removing problematic samples. Analyze the model output to identify potential signs of overfitting and collinearity between variables. Justify your observations. (4 Marks)

library(tidyverse)
library(MASS)

load("~/Coursework0125_Data.Rdata")
idx <- 20
c_cl <- setdiff(1:30, c(3, 18))
x <- patient_data$tissue[c_cl]
z <- patient_data$patient[c_cl]
tmp <- data.frame(y = Y[idx, c_cl], x = x, z = z, lib_size = colSums(Y[, c_cl]))
out <- glm.nb(y ~ x + lib_size, data = tmp)
p_val <- summary(out)$coefficients[2, 4]

summary(out)

Answer

Q6. Employing the code above, perform a Negative Binomial regression analysis that includes the following covariates: Tissue type and library size. Use your raw data after removing problematic samples. Compute the p-value and Bonferroni-adjusted p-value for the tissue type covariate for each gene. Capture and record any warnings or fitting issues that arise during model fitting (you can use the function withCallingHandlers() which allows you to intercept warnings as they occur inside a loop, without stopping the entire process). Report how many genes generated warnings. For genes with warnings, extract their adjusted p-values. Briefly discuss possible reasons why some genes generated warnings. Should genes with warnings be filtered out for further analyses, or should they be retained with caution? Provide a brief justification for your decision. (5 Marks)

Answer

Q7. Plot \(-log10(p-value)\) vs gene index and \(-log10(adjusted p-value)\) vs gene index for all the genes. Compare two plots and explain any visible differences or similarities between the raw p-value and adjusted p-value distributions. Justify your observations. Using a threshold of 0.05 for adjusted p-values, identify genes significantly associated with the tissue type. You may add a line on your plots to visulaize this. Using a code, find how many such genes you have? For your plots, use a different color to highlight genes that generated warnings during model fitting. What do you observe? (5 Marks)

Answer

Q8. As an alternative modeling approach, consider using simple linear regression on log-transformed expression values. Fit two models for each gene: (i) \(log2(CPM+1)\sim tissue\) and (ii) \(log2(CPM+1)\sim tissue + patient ID\) , compute the Bonferroni-adjusted p-value for the tissue type covariate for each gene for the two models. Create a scatter plot comparing the \(-log10(adjustedPvalues)\) from the two models (Model 1 vs. Model 2). Add a reference line y=x. Comment on the differences observed. Is it appropriate to include patient ID in the model? Why? Besides, test the homoscedasticity assumptions of the linear regression models for Model 1, comment on your findings. Now think about what is the advantage of NB-GLM in our case. (5 Marks)

Answer

Q9 (Optional). Co-expression Analysis Using Cleaned log2(CPM+1) Data: Firstly, select the top 100 most variable genes across all samples. 1. Compute the Pearson correlation matrix separately for Normal and Tumour samples. Visualize each correlation matrix using a heatmap (e.g., pheatmap, you may use its default clustering options). Identify and report the top 20 gene pairs with the largest absolute difference in correlation between Normal and Tumour tissues. Discussion: What potential issues arise if we only focus on the pairs with the largest correlation differences? 2. For each patient, compute the expression difference for each gene between Tumour and Normal tissue. Using these differences, compute the Pearson correlation matrix across genes and visualize it with a heatmap. Interpretation: Compare the two approaches above. What biological questions does each method address? How do they differ in what they reveal about gene co-expression? If we compute the correlation between the normal samples and tumour samples for each gene, could you try to intepret this?

Answer