R for Bioinformatics Assessment
Student ID: Pavankumar_W9646246 2024-05-18 1.
Answer: R, a versatile programming language, is highly regarded for its statistical analysis and data visualization capabilities. Researchers and data scientists frequently utilize it for tasks like data analysis, statistical modeling, and creating informative graphics.
Strengths of R:
Extensive Package Ecosystem: R boasts a vast library of packages available through CRAN, Bioconductor, and GitHub, expanding its functionality into various domains such as finance and genomics. Powerful Statistical Analysis: R excels at statistical analysis, providing a wide array of techniques for tasks including modeling, hypothesis testing, and time series analysis. High-Quality Graphics: R is renowned for its flexible and visually appealing graphical capabilities, enabling users to create insightful data visualizations. Efficient Data Handling: R offers a comprehensive suite of tools for importing data from various sources, cleaning messy datasets, and transforming data for analysis. Versatile Programming: R supports both functional and object-oriented programming paradigms, offering users the flexibility to choose the approach that best aligns with their needs. Limitations of R:
Memory-Intensive Operations: R tends to load data into memory, potentially causing bottlenecks when dealing with extremely large datasets. Security Concerns: R lacks built-in security features, making it less suitable for web applications where security is paramount. Steep Learning Curve: R’s syntax and paradigms may pose challenges for beginners, particularly those without prior programming or statistics experience. Core Data Structures in R:
Vectors: The fundamental building blocks in R, holding ordered sequences of elements of a single data type. Lists: Flexible containers capable of storing elements of various types, including other lists. Data Frames: The go-to structure for tabular data, with named columns that can accommodate different data types. Matrices: Two-dimensional arrays containing elements of the same data type, commonly used for mathematical operations and statistical modeling. Arrays: Multi-dimensional extensions of matrices for organizing data with multiple attributes. Factors: Designed for categorical data, storing values alongside distinct levels (categories). Essential Data Types in R:
Numeric: Represents a wide range of numerical values, including real numbers (with decimals) and integers (whole numbers). Integer: Specifically for whole numbers. Complex: Represents complex numbers with both real and imaginary components. Character (Strings): Stores text data. Logical (Boolean): Represents TRUE or FALSE values, essential for logical operations and conditional statements.
##4 7 2 ##5 9 2 ##6 11 2 ##7 1 4 ##8 3 4 ##9 5 4 ##10 7 4 ##11 9 4 ##12 11 4 ##13 1 6 ##14 3 6 ##15 5 6 ##16 7 6 ##17 9 6 ##18 11 6 ##19 1 8 ##20 3 8 ##21 5 8 ##22 7 8 ##23 9 8 ##24 11 8 ##25 1 10 ##26 3 10 ##27 5 10 ##28 7 10 ##29 9 10 ## 30 11 10 ##31 1 12 ##32 3 12 ##33 5 12 ##34 7 12 ##35 9 12 ## 36 11 12 b. uses your own function (not its internal function) to calculate the average value of all the sum of combinations. Answer: calculate_average_sum <- function(odd_numbers, even_numbers) { all_combinations <- expand.grid(odd_numbers, even_numbers) sum_combinations <- rowSums(all_combinations) average_sum <- mean(sum_combinations) return(average_sum) } average <- calculate_average_sum(odd_numbers, even_numbers) print(paste(“The average of all the combinationsn is”, {average})) ## [1] “The average of all the combinationsn is 13” 3
Make a function to create 20 moving-sum numbers: 1, 2, 3, 5, 8, 13, 21, 34, . . . How many even numbers? What is the total sum of these even numbers? Answer: #Definifn a moving sum function moving_sum <- function(n, numbers = c(1, 1)) { if (length(numbers) == n) { return(list(numbers = numbers, even_count = sum(numbers %% 2 == 0), even_sum = sum(numbers[numbers %% 2 == 0]))) } else { next_number <- numbers[length(numbers)] + numbers[length(numbers) - 1] return(moving_sum(n, c(numbers, next_number))) } } #Printing Results for 21 instead of 20 to accomodate two ones at the start result2 <- moving_sum(21) print(result2) ## $numbers ## [1] 1 1 2 3 5 8 13 21 34 55 89 144 ## [13] 233 377 610 987 1597 2584 4181 6765 10946 ## ## $even_count ## [1] 7 ## ## $even_sum ## [1] 14328 6
Tables GSE9476_series_matrix.txt (tab-delimited) and GSE9476_meta_info.txt are the microarray gene expression dataset and sample information table about patients with Acute Myeloid Leukemia (AML) and normal subjects as Controls. The expression data have been log2 transformed. Please complete the following tasks.
How many patients are with AML and how many are the control cases? How many genes(probesets) are we studying? (For the purpose of the assignment, each probeset is treated as a gene). Answer: # Data Loading expression_data <- read.table(“GSE9476_series_matrix.txt”, header = TRUE, row.names = 1, c meta_info <- read.table(“GSE9476_meta_info.txt”, header = TRUE, sep = “, fill= TRUE) #Counting fpr AML & Control aml_count <- sum(meta_info\(Disease_status == "AML") control_count <- sum(meta_info\)Disease_status ==”Control”) genes_count <- nrow(expression_data) #Printing the results cat(“Number of AML patients:”, aml_count, “”) ## Number of AML patients: 26 cat(“Number of control cases:”, control_count, “”) ## Number of control cases: 38 cat(“Number of genes/probesets:”, genes_count) ## Number of genes/probesets: 22283
Create your R function to calculate the average age of AML patients and Control subjects separately. Check your results by using R mean function. Answer: omment.char = #Definig calculate_average_age to find the average age of groups calculate_average_age <- function(data) { # Filtering data for AML and Control groups aml_data <- data[data\(Disease_status == "AML" & !is.na(data\)Age), ] control_data <- data[data\(Disease_status == "Control" & !is.na(data\)Age), ] # Calculating average age for each group 7
avg_age_aml <- if (nrow(aml_data) > 0) { mean(aml_data\(Age) } else { "No AML patients or missing age data" } avg_age_control <- if (nrow(control_data) > 0) { mean(control_data\)Age) } else { “No control subjects or missing age data” } # Returning results return(list( avg_age_AML = avg_age_aml, avg_age_Control = avg_age_control )) } # Calculating average calculate_average_age(meta_info) ## $avg_age_AML ## [1] 56.11538 ## ## \(avg_age_Control ## [1] 33.33333 ## [1] 56.11538 mean(meta_info\)Age[meta_info\(Disease_status == "Control"], na.rm = TRUE) # For Control subjects ## [1] 33.33333 The results above demonstrate that both the custom-written function and the built-in mean() function produced identical outcomes. 3. Is there any significant difference between sample ages of AMLs and Controls? Assume the age follows a normal distribution. Answer: # Verification using mean() function mean(meta_info\)Age[meta_info\(Disease_status == "AML"], na.rm = TRUE) # For AML patients #Applying t-test for significance t_test_result <- t.test(meta_info\)Age[meta_info\(Disease_status == "AML"], meta_info\)Age[meta_info$Disease_status == “Control”], print(t_test_result) na.action = na.exclude) 8
The exceptionally low p-value of 9.743e-08, significantly below the 0.05 threshold, provides compelling ev- idence to reject the null hypothesis. This indicates a statistically significant difference in age between the AML sample and the control group. 4. Make a boxplot of gene expression for all the samples, marking data points in red for AML samples and blue for Control samples. Answer: 9 #Loading library ggplot2 for boxplot library(ggplot2) # Creating a data frame for plotting boxplot plot_data <- data.frame( SampleID = colnames(expression_data), Gene_Expression = unlist(expression_data[2, ]), Group = factor(meta_info$Disease_status) ) # Creating the boxplot ggplot(plot_data, aes(x = Group, y = Gene_Expression, color = Group)) + geom_boxplot(outlier.shape = NA) + geom_jitter(width = 0.2, size = 2) + scale_color_manual(values = c(“AML” = “red”, “Control” = “blue”)) + labs(title = “Gene Expression in AML and Control Samples”, x = “Group”, y = “Gene Expression (Log2)”) + theme_classic()
8.0 7.5 7.0 6.5 6.0 Gene Expression in AML and Control Samples AML Control Group Group AML Control Gene Expression (Log2) 5. Create a table to include mean_AML, sd_AML, min_AML, max_AML, mean_Control, sd_Control, min_Control, max_Control, and Fold_change (i.e, mean_AML – mean_Control) for each gene. It is fine to use built-in functions. Hints: split the dataset into AML data and normal data sets, and then for each gene/probeset, calculate its mean, standard derivation, min, and max expression values across samples separately (using a built-in function), and further merge these statistical values for each gene into one table. Apply data.frame() and give the created table the same row names as the gene expression data table. Answer: # Spliting the data into disticnt groups aml_data <- expression_data[, meta_info\(Disease_status == "AML"] control_data <- expression_data[, meta_info\)Disease_status == “Control”] # Calculating all statistics for AML mean_AML <- apply(aml_data, 1, mean, na.rm = TRUE) sd_AML <- apply(aml_data, 1, sd, na.rm = TRUE) min_AML <- apply(aml_data, 1, min, na.rm = TRUE) max_AML <- apply(aml_data, 1, max, na.rm = TRUE) # Calculating all statistics for Control mean_Control <- apply(control_data, 1, mean, na.rm = TRUE) sd_Control <- apply(control_data, 1, sd, na.rm = TRUE) min_Control <- apply(control_data, 1, min, na.rm = TRUE) 10
max_Control <- apply(control_data, 1, max, na.rm = TRUE) # Calculate the fold change fold_change <- mean_AML - mean_Control # Combining all results into a data frame summary_table <- data.frame( mean_AML = mean_AML, sd_AML = sd_AML, min_AML = min_AML, max_AML = max_AML, mean_Control = mean_Control, sd_Control = sd_Control, min_Control = min_Control, max_Control = max_Control, fold_change = fold_change ) # Setting row names rownames(summary_table) <- rownames(expression_data) # Displaing the first 6 rows of the summary table head(summary_table, 6) ## mean_AML sd_AML min_AML max_AML mean_Control sd_Control ## 1007_s_at 3.190831 0.28688658 2.738629 4.189943 3.109452 0.13263540 6.671247 0.54927460 6.475112 2.63384079 7.176393 0.26198848 2.239653 0.07173247 7.715736 0.52599571 ## 1294_at ## ## 1007_s_at ## 1053_at ## 117_at ## 121_at ## 1255_g_at ## 1294_at 7.870432 0.63481967 6.757480 9.531410 min_Control max_Control fold_change 6.497534 0.43710532 5.707238 7.603169 5.352038 1.54329303 3.961892 9.644181 7.449728 0.35013512 6.665644 8.221074 ## 1053_at ## 117_at ## 121_at ## 1255_g_at 2.265448 0.07559105 2.068680 2.398069 2.843774 5.932384 3.923634 6.771577 2.099221 6.706907 3.448303 0.08137906 8.143075 -0.17371329 12.358038 -1.12307397 8.019825 0.27333450 2.512393 0.02579531 8.929810 0.15469537 6. With the table created above, make a scatterplot of the mean expression of AMLvs. Control samples. Are mean_AML and mean_Control significantly different? Answer: #Loading library ggplot2 for scastterplot library(ggplot2) # Scatter plot of mean AML and mean Control expression ggplot(summary_table, aes(x = mean_Control, y = mean_AML)) + geom_point() + 11
geom_abline(intercept = 0, slope = 1, color = “red”, linetype = “dashed”) + labs(title = “Scatterplot of Mean Gene Expression: AML vs. Control”, x = “Mean Control Expression”, y = “Mean AML Expression”) + theme_classic() Scatterplot of Mean Gene Expression: AML vs. Control 15 10 5 5 10 15 Mean Control Expression Mean AML Expression # Paired t-test to compare mean_AML and mean_Control t_test <- t.test(summary_table\(mean_AML, summary_table\)mean_Control, paired = TRUE) # Print t-test results print(t_test) ## ## Paired t-test ## ## data: summary_table\(mean_AML and summary_table\)mean_Control ## t = -6.8546, df = 22282, p-value = 7.336e-12 ## alternative hypothesis: true mean difference is not equal to 0 ## 95 percent confidence interval: ## -0.02303778 -0.01279226 ## sample estimates: ## mean difference ## -0.01791502 The combined evidence from the plot and t-test results strongly suggests a significant difference in mean gene expression between the AML and control groups. 12
# Counting the number of filtered genes num_filtered_genes <- nrow(filtered_genes) # Printing the result cat(“The number of genes/probsets left after filtering:”, num_filtered_genes) ## The number of genes/probsets left after filtering: 16026 b.) max_AML < 14 or max_Control < 14 Answer: # Filtering genes where max_AML < 14 OR max_Control < 14 filtered_genes <- summary_table[summary_table\(max_AML < 14 | summary_table\)max_Control < 1 # Counting the number of filtered genes num_filtered_genes <- nrow(filtered_genes) # Printing the result cat(“The number of genes/probsets left after filtering:”, num_filtered_genes) ## The number of genes/probsets left after filtering: 22171 c)Apply a) and b) filters together which is more stringent than applying a) or b) alone Answer: 4, ] # Combining filtering: mean_AML > 3 AND mean_Control > 3 AND (max_AML < 14 OR max_Control filtered_genes <- summary_table[ (summary_table\(mean_AML > 3 & summary_table\)mean_Control > 3) & (summary_table\(max_AML < 14 | summary_table\)max_Control < 14),] # Counting the number of filtered genes num_filtered_genes <- nrow(filtered_genes) # Printing the result cat(“The number of genes/probsets left after combined filtering:”, num_filtered_genes) ## The number of genes/probsets left after combined filtering: 15914 9. Following 8, if we further filter gene/probeset by applying p-value < 0.05 and the absolute value of log2 Fold_change bigger than 1, then how many genes will be left which we call them significantly expressed genes (DEGs)? Among them, how many are up-regulated (Fold_change > 0) and how many are down- regulated (Fold_change < 0)? List top 10 most upregulated DEGs and 10 most down-regulated DEGs including their p-value, log2FC, and mean values in each sample group. Answer: 14 < 14)
# Filtering for significant DEGs significant_DEGs <- summary_table[summary_table\(p_value_adj < 0.05 & abs(summary_table\)fol # Counting DEGs num_DEGs <- nrow(significant_DEGs) cat(“Number of differentially expressed genes (DEGs):”, num_DEGs, “”) ## Number of differentially expressed genes (DEGs): 576 ## Number of up-regulated DEGs: 117 ## Number of down-regulated DEGs: 459 ## [1] “Top 10 Up-regulated DEGs:” ## mean_AML mean_Control fold_change p_value_adj 1:10, ] ## 215489_x_at 6.237537 3.431159 6.876417 4.736884 4.096730 2.935262 4.054346 9.891301 4.204067 7.931564 6.624619 2.806378 1.463890e-05 2.793421 4.588981e-07 2.704195 1.038839e-06 2.667401 1.526208e-04 2.415503 1.047242e-03 2.337609 2.249909e-04 2.236181 6.339470e-06 2.150490 1.568099e-02 2.146922 1.557072e-02 2.107120 8.747943e-03 ## 206674_at ## 210365_at ## 204647_at ## 205899_at ## 206067_s_at 6.391956 ## 201105_at 12.127482 ## 205131_x_at 6.354558 ## 205382_s_at 10.078486 ## 211709_s_at 8.731739 9.669838 7.441080 6.764131 5.350765 ## [1] “Top 10 Down-regulated DEGs:” 15 d_change) > 1, # Counting up-regulated and down-regulated DEGs num_up_DEGs <- sum(significant_DEGs\(fold_change > 0) num_down_DEGs <- sum(significant_DEGs\)fold_change < 0) cat(“Number of up-regulated DEGs:”, num_up_DEGs, “”) cat(“Number of down-regulated DEGs:”, num_down_DEGs, “”) # Sorting and list top 10 up- and down-regulated DEGs top_up_DEGs <- significant_DEGs[order(significant_DEGs\(fold_change, decreasing = TRUE), ][ top_down_DEGs <- significant_DEGs[order(significant_DEGs\)fold_change),][1:10, ] print(“Top 10 Up-regulated DEGs:”) print(top_up_DEGs[, c(“mean_AML”, “mean_Control”, “fold_change”, “p_value_adj”)]) print(“Top 10 Down-regulated DEGs:”) print(top_down_DEGs[, c(“mean_AML”, “mean_Control”, “fold_change”, “p_value_adj”)])
12.489345 11.578075 8.710812 12.213326 11.419415 11.479616 11.758847 8.090973 11.802806 11.726737 -3.950883 1.538226e-04 -3.837153 1.593589e-03 -3.699978 8.645927e-07 -3.545056 6.262285e-05 -3.535036 1.624928e-03 -3.516957 1.593589e-03 -3.458872 1.695677e-03 -3.452333 3.077606e-04 -3.414876 1.562872e-03 -3.356328 2.784251e-03 10. If you apply the limma package rather than t-test, then how many DEGs will be discovered by limma? (p-value < 0.05 & abs(log2FC) > 1). How many over- lapped DEGs by these two different approaches? please draw a Venn diagram for demonstration. Answer: ## Loading required package: grid ## Loading required package: futile.logger plot.new() #Loading limma and venndiagram library library(limma) library(VennDiagram) # Preparing design matrix for limma group <- factor(meta_info\(Disease_status) design <- model.matrix(~0 + group) colnames(design) <- levels(group) # Fitting linear model fit <- lmFit(expression_data, design) # Seting up contrast matrix between both groups contrast_matrix <- makeContrasts(AML - Control, levels = design) # Fitting contrasts fit2 <- contrasts.fit(fit, contrast_matrix) # Computing moderated t-statistics and log-odds of differential expression fit2 <- eBayes(fit2) # Getting results limma_results <- topTable(fit2, adjust.method = "BH", n = Inf) # Filtering DEGs based on p-value < 0.05 and |log2FC| > 1 limma_DEGs <- limma_results[limma_results\)adj.P.Val < 0.05 & abs(limma_results$logFC) > 1, 16 ]
# Finding overlapped DEGs t_test_DEGs <- rownames(summary_table)[summary_table\(p_value_adj < 0.05 & abs(summary_tabl overlapped_DEGs <- intersect(rownames(limma_DEGs), t_test_DEGs) # Numbering of DEGs num_limma_DEGs <- nrow(limma_DEGs) cat("Number of DEGs found by limma:", num_limma_DEGs, "\n") ## Number of DEGs found by limma: 572 ## Number of overlapping DEGs: 544 ## [1] 1 17 e\)fold_change) num_overlapped <- length(overlapped_DEGs) cat(“Number of overlapping DEGs:”, num_overlapped, “”) # Venn Diagram venn.diagram( x = list( ‘T-test’ = t_test_DEGs, ‘Limma’ = rownames(limma_DEGs) ), filename = ‘DEG_overlap.png’, imagetype = ‘png’, fill = c(‘red’, ‘blue’), alpha = 0.5, cat.col = c(‘orange’, ‘blue’), cat.cex = 1.5, margin = 0.1 ) # Displaying the Venn diagram using grid.raster library(grid) img <- png::readPNG(“DEG_overlap.png”) grid.raster(img)
The Venn diagram reveals an overlap of approximately 544 differentially expressed genes (DEGs) identified by both approaches. 18