20220719_CMFI_Seminar_Data_Cleanup_Prelim_Stats.ipynb

{"cells":[{"cell_type":"markdown","id":"55766981","metadata":{"id":"55766981"},"source":["## Data Clean up\n","Authors: Abzer Kelminal (abzer.shah@uni-tuebingen.de), Madeleine Ernst(MAET@ssi.dk) <br>\n","Edited by: Daniel Petras (daniel.petras@uni-tuebingen.de)  <br>\n","Input file format: .csv files or .txt files <br>\n","Outputs: .csv files  <br>\n","Dependencies: ggplot2, dplyr, ecodist, vegan, svglite"]},{"cell_type":"code","execution_count":null,"id":"BvfG-ihq2PZP","metadata":{"id":"BvfG-ihq2PZP"},"outputs":[],"source":["#installing and calling the necessary packages:\n","install.packages(\"ggplot2\")\n","install.packages(\"dplyr\")\n","install.packages(\"ecodist\") #for PCoA using Bray Curtis distance\n","install.packages(\"vegan\") #for PermANOVA\n","install.packages(\"svglite\") # for saving ggplots as svg files"]},{"cell_type":"code","source":["require(\"ggplot2\")\n","require(\"dplyr\")\n","require(\"ecodist\")\n","require(\"vegan\")\n","require(\"svglite\")"],"metadata":{"id":"-Rym0mXNMZ7B"},"id":"-Rym0mXNMZ7B","execution_count":null,"outputs":[]},{"cell_type":"markdown","id":"863a44d4","metadata":{"id":"863a44d4"},"source":["## Setting a local working directory and creating an automatic result directory:\n","Works well with Jupyter Notebook. \n","For Google Collab, we can upload the necessary files into a new folder using the 'Files' icon on the left and set the folder as working directory."]},{"cell_type":"code","execution_count":null,"id":"4734f4b6","metadata":{"id":"4734f4b6"},"outputs":[],"source":["# setting the current directory as the working directory\n","Directory <- normalizePath(readline(\"Enter the path of the folder with input files: \"),\"/\",mustWork=FALSE)\n","setwd(Directory)"]},{"cell_type":"code","source":["getwd()"],"metadata":{"id":"og4ryyATcSju"},"id":"og4ryyATcSju","execution_count":null,"outputs":[]},{"cell_type":"code","execution_count":null,"id":"dd010938","metadata":{"id":"dd010938"},"outputs":[],"source":["# Getting all the files in the folder\n","dirs <- dir(path=paste(getwd(), sep=\"\"), full.names=TRUE, recursive=TRUE)\n","folders <- unique(dirname(dirs))\n","files <- list.files(folders, full.names=TRUE)\n","files_1 <- basename((files))\n","files_2 <- dirname((files))\n","# Creating a Result folder\n","dir.create(path=paste(files_2[[1]], \"_Results\", sep=\"\"), showWarnings = TRUE)\n","fName <-paste(files_2[[1]], \"_Results\", sep=\"\")\n","\n","print(files_1)"]},{"cell_type":"markdown","id":"55ab1133","metadata":{"id":"55ab1133"},"source":["**<font color='red'> In the following line, enter the required file ID numbers separated by commas. For example as: 1,2,3 </font>**"]},{"cell_type":"code","execution_count":null,"id":"45d07e82","metadata":{"id":"45d07e82"},"outputs":[],"source":["input <- as.double(unlist(strsplit(readline(\"Specify the file index of gapfilled & non-gapfilled feature-file, metadata:\"), split=\",\")))\n","\n","#Gets the extension of each file. Ex:csv\n","pattern <- c()\n","for (i in files_1){\n","  sep_file <- substr(i, nchar(i)-2,nchar(i))\n","  pattern <- rbind(pattern,sep_file)\n","}\n","#pattern\n","\n","ft <- read.csv(files_1[input[1]],sep = ifelse(pattern[input[1]]!=\"csv\",\"\\t\",\",\"), header=TRUE,check.names = FALSE) # By applying 'row.names = 1', the 1st column 'ID' becomes the row names\n","nft<- read.csv(files_1[input[2]],sep=ifelse(pattern[input[2]]!=\"csv\",\"\\t\",\",\"), header = TRUE,check.names = FALSE)\n","md <-read.csv(files_1[input[3]], sep = ifelse(pattern[input[3]]!=\"csv\",\"\\t\",\",\"), header=TRUE,check.names = FALSE)"]},{"cell_type":"markdown","id":"3a31eafc","metadata":{"id":"3a31eafc"},"source":["## Reading the input data using URL (from Github):"]},{"cell_type":"code","execution_count":null,"id":"7988ffd6","metadata":{"id":"7988ffd6"},"outputs":[],"source":["## Non-gap filled\n","nft_url <- 'https://raw.githubusercontent.com/Functional-Metabolomics-Lab/CMFI_Seminar_Multivariate_Statistics/main/Test_Data/20220716_Xenobiotic_metabolism_non_gapfilled_quant_Bsub_quant.csv'\n","## Gap filled\n","ft_url <- 'https://raw.githubusercontent.com/Functional-Metabolomics-Lab/CMFI_Seminar_Multivariate_Statistics/main/Test_Data/20220716_Xenobiotic_metabolism_gapfilled_quant_Bsub.csv'\n","md_url <- 'https://raw.githubusercontent.com/Functional-Metabolomics-Lab/CMFI_Seminar_Multivariate_Statistics/main/Test_Data/20220716_Xenobiotic_Metabolism_metadata_Bsub.txt'"]},{"cell_type":"code","execution_count":null,"id":"3ad9e910","metadata":{"id":"3ad9e910"},"outputs":[],"source":["nft <- read.csv(nft_url, header = T, check.names = F)\n","ft <- read.csv(ft_url, header = T, check.names = F)\n","md <- read.csv(md_url, header = T, check.names = F, sep = '\\t')"]},{"cell_type":"markdown","id":"4f12e24c","metadata":{"id":"4f12e24c"},"source":["Lets check if the data has been read correclty!!"]},{"cell_type":"code","execution_count":null,"id":"429ff705","metadata":{"id":"429ff705"},"outputs":[],"source":["head(ft)\n","dim(ft)"]},{"cell_type":"code","execution_count":null,"id":"864b7231","metadata":{"id":"864b7231"},"outputs":[],"source":["head(nft)\n","dim(nft)"]},{"cell_type":"code","execution_count":null,"id":"a7bf0865","metadata":{"id":"a7bf0865"},"outputs":[],"source":["head(md)\n","dim(md)"]},{"cell_type":"markdown","id":"f19073e0","metadata":{"id":"f19073e0"},"source":["Trying to bring the feature table and metadata in the correct format:"]},{"cell_type":"code","execution_count":null,"id":"25ffd93c","metadata":{"id":"25ffd93c"},"outputs":[],"source":["#Removing Peak area extensions\n","colnames(ft) <- gsub(' Peak area','',colnames(ft))\n","colnames(nft) <- gsub(' Peak area','',colnames(nft))\n","md$filename<- gsub(' Peak area','',md$filename)\n","\n","#Removing if any NA columns present in the md file\n","ft <- ft[,colSums(is.na(ft))<nrow(ft)]\n","nft <- nft[,colSums(is.na(nft))<nrow(nft)]\n","md <- md[,colSums(is.na(md))<nrow(md)]\n","\n","#Changing the row names of the files\n","rownames(md) <- md$filename\n","md <- md[,-1]\n","rownames(ft) <- paste(ft$'row ID',round(ft$'row m/z',digits = 3),round(ft$'row retention time',digits = 3), sep = '_')\n","rownames(nft) <- paste(nft$'row ID',round(nft$'row m/z',digits = 3),round(nft$'row retention time',digits = 3), sep = '_')\n","\n","#Picking only the files with column names containing 'mzML'\n","ft <- ft[,grep('mzML',colnames(ft))]\n","nft <- nft[,grep('mzML',colnames(nft))]\n","\n","# Converting replicate attributes into factors (categorical data)\n","md$ATTRIBUTE_replicates <- as.factor(md$ATTRIBUTE_replicates)"]},{"cell_type":"markdown","id":"e927cf45","metadata":{"id":"e927cf45"},"source":["Lets check the files once again!!"]},{"cell_type":"code","execution_count":null,"id":"e76fb0e3","metadata":{"id":"e76fb0e3"},"outputs":[],"source":["head(nft)\n","dim(nft)"]},{"cell_type":"code","execution_count":null,"id":"499d8a48","metadata":{"id":"499d8a48"},"outputs":[],"source":["head(ft)\n","dim(ft)"]},{"cell_type":"code","execution_count":null,"id":"55093982","metadata":{"id":"55093982"},"outputs":[],"source":["head(md)\n","dim(md)"]},{"cell_type":"markdown","id":"96275453","metadata":{"id":"96275453"},"source":["### Creating a function named FrequencyPlot:  \n","The below function takes in the two input datatables: for example, gapfilled and non-gapfilled, calculates the frequency distribution of the data in the order of 10 and produces a grouped barplot showing the distribution as output. The frequency plot shows where the features are present in higher number."]},{"cell_type":"code","execution_count":null,"id":"e6d2b3b4","metadata":{"id":"e6d2b3b4"},"outputs":[],"source":["#'Global' settings for plot size in the output cell\n","options(repr.plot.width=10, repr.plot.height=8,res=600) #For google collab\n","#options(repr.plot.width=5, repr.plot.height=3) #For Jupyter Notebook"]},{"cell_type":"code","execution_count":null,"id":"4755a764","metadata":{"id":"4755a764"},"outputs":[],"source":["FrequencyPlot <- function(x1,x2){\n","  \n","   #creating bins from -1 to 10^10 using sequence function seq()\n","    bins <- c(-1,0,(1 * 10^(seq(0,10,1)))) \n","    \n","    #cut function cuts the give table into its appropriate bins\n","    scores_x1 <- cut(as.matrix(x1),bins,labels = c('0','1','10','1E2','1E3','1E4','1E5','1E6','1E7','1E8','1E9','1E10')) \n","    \n","    #transform function convert the tables into a column format: easy for visualization \n","    Table_x1<-transform(table(scores_x1)) #contains 2 columns: \"scores_x1\", \"Freq\"\n","    \n","    #Repeating the same steps for x2\n","    scores_x2 <- cut(as.matrix(x2),bins,labels = c('0','1','10','1E2','1E3','1E4','1E5','1E6','1E7','1E8','1E9','1E10'))\n","    Table_x2<-transform(table(scores_x2))\n","  \n","    #Getting the names of x1 and x2\n","    arg1 <- deparse(substitute(x1))\n","    arg2 <-deparse(substitute(x2))\n","    \n","    #Creating a data frame for plotting\n","    data_plot <- as.data.frame(c(Table_x1$Freq,Table_x2$Freq)) #Concatenating the frequency info of both tables rowwise\n","    colnames(data_plot) <- \"Freq\" #naming the 1st column as 'Freq'\n","    data_plot$Condition <- c(rep(arg1,12),rep(arg2,12)) #adding a 2nd column 'Condition', which just repeats the name of x1 and x2 accordingly\n","    data_plot$Range_bins <- rep(Table_x1$scores_x1,2) #Adding 3rd column 'Range Bins'\n","    data_plot$Log_Freq <- log(data_plot$Freq+1) #Log scaling the frequency values\n","    \n","    ## GGPLOT2\n","    BarPlot <- ggplot(data_plot, aes(Range_bins, Log_Freq, fill = Condition)) + \n","    geom_bar(stat=\"identity\", position = \"dodge\", width=0.4) + \n","    scale_fill_brewer(palette = \"Set1\") +\n","    ggtitle(label=\"Frequency plot\") +\n","    xlab(\"Range\") + ylab(\"(Log)Frequency\") + labs(fill = \"Data Type\") + \n","    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +   # setting the angle for the x label\n","    theme(axis.text.y = element_text(angle = 45, vjust = 0.5, hjust=1)) +   # setting the angle for the y label\n","    theme(plot.title = element_text(hjust = 0.5)) # centering the plot title\n","  \n","    print(BarPlot)\n","}  "]},{"cell_type":"markdown","id":"ea08fc3d","metadata":{"id":"ea08fc3d"},"source":["**About the experiment:**\n","- Bacteria (B.subtilis) was treated with a pool of antibiotics (Sulfamethoxazole, sulfadimethoxine, cyproconazole) including a herbicide Asulam, taken at a concentration lower than their MIC (minimum inhibitory concentration).\n","- The samples were collected at different timepoints, the compounds were extracted (with 50% EtOAc) and measured using LC-MS/MS.\n","- The goal of the experiment was to look for any potential biotransformation. eg: Drug or xenobiotic metabolism"]},{"cell_type":"markdown","id":"fa549984","metadata":{"id":"fa549984"},"source":["## Splitting the data into Control and Samples using Metadata:"]},{"cell_type":"code","source":["input_data <- ft"],"metadata":{"id":"F7LHJ3dOf1EO"},"id":"F7LHJ3dOf1EO","execution_count":null,"outputs":[]},{"cell_type":"code","execution_count":null,"id":"f2dc8ca3","metadata":{"id":"f2dc8ca3"},"outputs":[],"source":["head(md)\n","print(matrix(data=colnames(md),nrow=length(colnames(md))))\n","\n","#These lines are not needed in R console, but in Jupyter Notebook to get the previous print statement working\n","flush.console()  \n","Sys.sleep(0.2)\n","\n","Condition <- as.double(unlist(readline(\"Enter the index of the attribute to split sample and control:\")))\n","\n","flush.console()  \n","Sys.sleep(0.2)\n","\n","Levels_Cdtn <- levels(as.factor(md[,Condition[1]]))\n","print(matrix(Levels_Cdtn,length(Levels_Cdtn)))\n","\n","flush.console()  \n","Sys.sleep(0.2)\n","    \n","#Among the shown levels of an attribute, select the ones to keep\n","Ctrl_id <- as.double(unlist(readline(\"Enter the index of your BLANK:\")))\n","paste0('You chosen blank is:',Levels_Cdtn[Ctrl_id])\n","\n","#Splitting the data into control and samples based on the metadata\n","md_Ctrl <- md[(md[,Condition] == Levels_Cdtn[Ctrl_id]),]\n","Ctrl <- input_data[,which(colnames(input_data)%in%rownames(md_Ctrl))] \n","md_Samples <- md[(md[,Condition] != Levels_Cdtn[Ctrl_id]),]\n","Samples <- input_data[,which(colnames(input_data)%in%rownames(md_Samples))] "]},{"cell_type":"code","execution_count":null,"id":"3c2c31f4","metadata":{"id":"3c2c31f4"},"outputs":[],"source":["head(Ctrl)\n","dim(Ctrl)"]},{"cell_type":"code","execution_count":null,"id":"85457a6b","metadata":{"id":"85457a6b"},"outputs":[],"source":["head(Samples)\n","dim(Samples)"]},{"cell_type":"code","execution_count":null,"id":"8b376c71","metadata":{"id":"8b376c71"},"outputs":[],"source":["FrequencyPlot(Samples,Ctrl)"]},{"cell_type":"markdown","id":"43099574","metadata":{"id":"43099574"},"source":["## Blank Removal:\n","\n","(Note: In LC-MS/MS, we use solvents also called as Blanks which are usually injected time-to-time to prevent carryover of the sample) </br>\n","\n","For the Blank removal step, we need to split the data as control blanks and samples. </br>\n","\n","**The blanks we are referring to here, is the control blanks in the experiment and not the LC-MS/MS blanks.**\n","- The control blanks here is the sample without treatment. \n","- Samples are biological replicates with treatment and we have two sets of data: B.sub and E.coli. </br>\n","\n","In general, having multiple control blanks helps us to compare any variation in the data. Comparing control to the sample helps us to identify the background features that contribute to any technical variation. A common filtering method is to use a cutoff to remove features that are not present sufficient enough in our biological samples.\n","\n","1. We find an average for all the feature intensities in your control set and sample set.\n","Therefore, for n no.of features in a control or sample set, we get n no.of averaged features.\n","2. Next, we get a ratio of this average_control vs average_sample. This ratio Control/sample tells us how much of that particular feature of a sample gets its contribution from control. If it is more than 30% (or Cutoff as 0.3), we consider the feature as noise.\n","3. The resultant information (if ratio > Cutoff or not) is stored in a bin such as **1 == Noise or background signal, 0 == Feature Signal**\n","4. We count the no.of features in the bin that satisfies the condition ratio > cutoff, and consider those features as 'noise or background features' and remove them.\n","\n","For a dataset containing several batches, the filtering steps are performed batch-wise."]},{"cell_type":"code","execution_count":null,"id":"af1aa72f","metadata":{"id":"af1aa72f"},"outputs":[],"source":["if(readline('Do you want to perform Blank Removal- Y/N:')=='Y'){\n","    \n","    #When cutoff is low, more noise (or background) detected; With higher cutoff, less background detected, thus more features observed\n","    Cutoff <- as.numeric(readline('Enter Cutoff value between 0.1 & 1:')) # (i.e. 10% - 100%). Ideal cutoff range: 0.1-0.3\n","    \n","    #Getting mean for every feature in Ctrl and Samples\n","    Avg_ctrl <- rowMeans(Ctrl, na.rm= FALSE, dims = 1) # set na.rm = FALSE to check if there are NA values. When set as TRUE, NA values are changed to 0\n","    Avg_samples <- rowMeans(Samples, na.rm= FALSE, dims = 1)\n","    \n","    #Getting the ratio of Ctrl vs Sample\n","    Ratio_Ctrl_Sample <- (Avg_ctrl+1)/(Avg_samples+1)\n","    \n","    # Creating a bin with 1s when the ratio>Cutoff, else put 0s\n","    Bg_bin <- ifelse(Ratio_Ctrl_Sample > Cutoff, 1, 0 )\n","    Blank_removal <- cbind(Samples,Bg_bin)\n","\n","    # Checking if there are any NA values present. Having NA values in the 4 variables will affect the final dataset to be created\n","    temp_NA_Count <-cbind(Avg_ctrl ,Avg_samples,Ratio_Ctrl_Sample,Bg_bin)\n","    \n","    print('No of NA values in the following columns:')\n","    print(colSums(is.na(temp_NA_Count)))\n","\n","     #Calculating the number of background features and features present\n","    print(paste(\"No.of Background or noise features:\",sum(Bg_bin ==1,na.rm = TRUE)))\n","    print(paste(\"No.of features after excluding noise:\",(nrow(Samples) - sum(Bg_bin ==1,na.rm = TRUE)))) \n","\n","    Blank_removal <- Blank_removal %>% filter(Bg_bin == 0) # Taking only the feature signals\n","    Blank_removal <- as.matrix(Blank_removal[,-ncol(Blank_removal)]) # removing the last column Bg_bin \n","}"]},{"cell_type":"code","source":["dim(Blank_removal)"],"metadata":{"id":"eh8Gd6QD_Dt3"},"id":"eh8Gd6QD_Dt3","execution_count":null,"outputs":[]},{"cell_type":"markdown","id":"b4210007","metadata":{"id":"b4210007"},"source":["## Imputation: \n","\n","For several reasons, real world datasets might have some missing values in it, in the form of NA, NANs or 0s. Eventhough the gapfilling step of MZmine fills the missing values, we still end up with some missing values or 0s in our feature table. This could be problematic for statistical analysis. \n","In order to have a better dataset, we cannot simply discard those rows or columns with missing values as we will lose a chunk of our valuable data.\n","Instead we can try imputing those missing values. Imputation involves replacing the missing values in the data with a meaningful, reasonable guess. There are several methods, such as:  \n","1) Mean imputation (replacing the missing values in a column with the mean or average of the column)  \n","2) Replacing it with the most frequent value  \n","3) Several other machine learning imputation methods such as k-nearest neighbors algorithm(k-NN), Hidden Markov Model(HMM)\n","\n","One such method, we are going to use is: **to replace the zeros from the gapfilled quant table with the non-gap filled table** we get from MZmine. In order to do that, we can visualize our data distribution using the frequenct plot"]},{"cell_type":"code","execution_count":null,"id":"9d83ae53","metadata":{"id":"9d83ae53"},"outputs":[],"source":["GapFilled <-Blank_removal\n","NotGapFilled <- nft"]},{"cell_type":"code","execution_count":null,"id":"e527abdb","metadata":{"id":"e527abdb"},"outputs":[],"source":["if(readline('Do you want to perform Imputation with minimum value of NonGapFilled table? - Y/N:')=='Y'){\n","    \n","    plot<- FrequencyPlot(GapFilled,NotGapFilled)\n","    \n","    Arg1 = plot$data$Condition[1]\n","    Arg2 = plot$data$Condition[13]\n","    \n","    # accessing the datatable of plot and subsetting with the condition: Eliminating the Range (or bin) 0 and Ranges with zero frequencies \n","    plotData_New <- subset(plot$data,plot$data$Freq!=0 & plot$data$Range_bins !=0) \n","    \n","    #getting the first appearing value of this new plot datatable\n","    First_val_temp <- aggregate(plotData_New$Freq, by=list(plotData_New$Condition), FUN=first) \n","    \n","    # Subsetting the rows in the plotData_New that has the first appearing values\n","    First_val <- plotData_New[plotData_New$Freq %in% c(First_val_temp$x[1],First_val_temp$x[2]),]\n","  \n","    # getting the 2nd minimum value of non-gap filled data. (The first minimum value in the data table is usually zero)\n","    RawLOD <- round(min(NotGapFilled[NotGapFilled!=min(NotGapFilled)]))\n","    print(paste0(\"The minimum value greater than 0 for \",Arg1,\":\", round(min(GapFilled[GapFilled!=min(GapFilled)]))))\n","    print(paste0(\"The minimum value greater than 0 for \",Arg2,\":\", RawLOD))\n","    \n","    Imputed <- GapFilled\n","    Imputed[Imputed<RawLOD] <- RawLOD # Replacing values<RawLOD with RawLOD\n","} else return(GapFilled)"]},{"cell_type":"code","execution_count":null,"id":"2698000e","metadata":{"id":"2698000e"},"outputs":[],"source":["write.csv(Imputed, file=paste0('Quant_Table_filled_with_MinValue_',RawLOD,'.csv'),row.names =TRUE) "]},{"cell_type":"code","execution_count":null,"id":"a25d5a9c","metadata":{"id":"a25d5a9c"},"outputs":[],"source":["head(Imputed)\n","dim(Imputed)"]},{"cell_type":"markdown","id":"f06a2dd7","metadata":{"id":"f06a2dd7"},"source":["## Normalization:\n","The following code performs sample-centric (column-wise) normalisation:"]},{"cell_type":"code","execution_count":null,"id":"3a008a7a","metadata":{"id":"3a008a7a"},"outputs":[],"source":["if (readline(\"Do you want to perform Normalization: Y/N:\") == 'Y'){\n","    \n","    #Getting column-wise sums of the input-data\n","    sample_sum <- colSums(Imputed, na.rm= TRUE, dims = 1)\n","    \n","    #Dividing each element of a particular column with its column sum\n","    Normalized_data <- c()\n","    for (i in 1:ncol(Imputed)){\n","        x <- Imputed[,i] / sample_sum[i]\n","        Normalized_data <- cbind(Normalized_data, x)\n","    }\n","    colnames(Normalized_data) <- names(sample_sum)\n","    \n","} else return(Imputed)\n","  \n","print(paste('No.of NA values in Normalized data:',sum(is.na(Normalized_data)== TRUE)))"]},{"cell_type":"code","execution_count":null,"id":"58fc407e","metadata":{"id":"58fc407e"},"outputs":[],"source":["write.csv(Normalized_data,file='Normalised_Quant_table.csv',row.names =TRUE) "]},{"cell_type":"markdown","id":"0d1d4bbc","metadata":{"id":"0d1d4bbc"},"source":["## Principal Coordinate analysis:"]},{"cell_type":"markdown","source":["Principal coordinates analysis (PCoA) is a metric multidimensional scaling (MDS) method that attempts to represent sample dissimilarities in a low-dimensional space. It converts a distance matrix consisting of pair-wise distances (dissimilarities) across samples into a 2- or 3-D graph ([Gower, 2005](https://doi.org/10.1002/0470011815.b2a13070)). Different distance metrics can be used to calculate dissimilarities among samples (e.g. Euclidean, Canberra, Minkowski). Performing a principal coordinates analysis using the Euclidean distance metric is the same as performing a principal components analysis (PCA). Selecting the best distance metric for a given dataset is part of the 'art' of data science.\n","\n","Within the metabolomics field the Euclidean, Bray-Curtis, Jaccard or Canberra distances are most commonly used. The Jaccard distance is an unweighted metric (presence/absence) whereas Euclidean, Bray-Curtis and Canberra distances take into account relative abundances (weighted). "],"metadata":{"id":"hpRW_hXmiPjs"},"id":"hpRW_hXmiPjs"},{"cell_type":"code","execution_count":null,"id":"420c9541","metadata":{"id":"420c9541"},"outputs":[],"source":["#Making sure the metadata rownames are identical to that of filenames in our featuretable in order to perform multivariate statistics\n","md_Stats <- md[which(rownames(md)%in%colnames(Normalized_data)),]\n","md_Stats <- md_Stats[match(colnames(Normalized_data),rownames(md_Stats)),]\n","identical(colnames(Normalized_data),rownames(md_Stats))\n","\n","#Checking the data sparsity (amount of zeros in our data matrix):\n","sum(Normalized_data == 0)/(dim(Normalized_data)[1]*dim(Normalized_data)[2])"]},{"cell_type":"code","execution_count":null,"id":"d0e4828b","metadata":{"id":"d0e4828b"},"outputs":[],"source":["#Metadata subsetting based on condition:\n","print(matrix(data=colnames(md_Stats),nrow=length(colnames(md_Stats))))\n","Condition <- as.double(unlist(strsplit(readline(\"Enter the IDs of interested attributes separated by commas to subset the data:\"),split=\",\")))\n","for(i in 1:length(Condition)){\n","  #Shows the different levels within each selected condition:\n","  Levels_Cdtn <- levels(as.factor(md_Stats[,Condition[i]]))\n","  print(matrix(Levels_Cdtn,length(Levels_Cdtn)))\n","  \n","  #These lines are not needed in R console, but in Jupyter Notebook to get the previous print statement working\n","  flush.console()  \n","  Sys.sleep(0.2)\n","  \n","  #Among the shown levels of an attribute, select the ones to keep\n","  Cdtn <- as.double(unlist(strsplit(readline(\"Enter the IDs of condition(s) you want to KEEP (separated by commas):\"), split=',')))\n","  Levels_Cdtn[Cdtn]\n","  \n","  #Selecting only rows in meta_filtered that match the condition\n","  md_Stats <- md_Stats[(md_Stats[,Condition[i]] == Levels_Cdtn[Cdtn]),]\n","}"]},{"cell_type":"code","execution_count":null,"id":"ff2bd0b0","metadata":{"id":"ff2bd0b0"},"outputs":[],"source":["dist_matrix <- bcdist(t(Normalized_data)) # transposed in order to compute the distance between the columns of a data matrix\n","pcoa<- cmdscale(dist_matrix, eig = TRUE, x.ret=TRUE)\n","pcoa.var.per <-round(pcoa$eig/sum(pcoa$eig)*100,1)\n","pcoa.values <- pcoa$points"]},{"cell_type":"markdown","source":["Calculating the pairwise distances across all samples using the Bray-Curtis distance metric:"],"metadata":{"id":"AURSgrlrjhcm"},"id":"AURSgrlrjhcm"},{"cell_type":"code","execution_count":null,"id":"0c3bce6e","metadata":{"id":"0c3bce6e"},"outputs":[],"source":["#PCoA calculation\n","md_data <- Normalized_data[,which(colnames(Normalized_data)%in%rownames(md_Stats))] # the corresponding column files for the filtered metadata is picked from the normalized data\n","\n","dist_matrix <- as.matrix(bcdist(t(md_data))) # transposed in order to compute the distance between the columns of a data matrix\n","pcoa<- cmdscale(dist_matrix, eig = TRUE, x.ret=TRUE)\n","pcoa.var.per <-round(pcoa$eig/sum(pcoa$eig)*100,1)\n","pcoa.values <- pcoa$points"]},{"cell_type":"code","execution_count":null,"id":"1cbcddff","metadata":{"id":"1cbcddff"},"outputs":[],"source":["at_int <- as.double(readline('Enter the index of your interested attribute for PCoA visualisation:'))"]},{"cell_type":"code","execution_count":null,"id":"fbb2b56f","metadata":{"id":"fbb2b56f"},"outputs":[],"source":["#PCoA plot:\n","pcoa.data <- data.frame(md_Stats[,at_int],\n","                        X=pcoa.values[,1],\n","                        Y=pcoa.values[,2])\n","\n","PCoA_plot <- ggplot(pcoa.data, aes(x=X, y=Y, col= as.factor(md_Stats[,at_int]))) + \n","  geom_point(size=4,alpha=0.8)  +\n","  ggtitle(label=\"MDS plot using Bray-Cutis Distance\") +\n","  xlab(paste0(\"MDS1 : \",pcoa.var.per[1],\"%\",sep=\"\")) + \n","  ylab(paste0(\"MDS2 : \",pcoa.var.per[2],\"%\",sep=\"\")) + \n","  labs(color = 'Timepoint') + \n","  theme(plot.title = element_text(hjust = 0.5)) \n","\n","PCoA_plot"]},{"cell_type":"markdown","id":"9557a7c1","metadata":{"id":"9557a7c1"},"source":["## Permutational multivariate analysis of variance (PERMANOVA):"]},{"cell_type":"markdown","source":["PERMANOVA is a non-parametric method for multivariate analysis of variance, where P-values are obtained using permutations. The metric was originally developed within the field of ecology ([Anderson, 2008](https://onlinelibrary.wiley.com/doi/10.1111/j.1442-9993.2001.01070.pp.x)) but is today widely used in other fields, including the microbiome and metabolomics field. PERMANOVA is used to compare groups of samples and tests whether the centroid and/or the spread of the samples is different between the groups. "],"metadata":{"id":"mzJvolcGj-6l"},"id":"mzJvolcGj-6l"},{"cell_type":"markdown","source":["The adonis2() function in the vegan package can be used to perform a PERMANOVA. The input is any dissimilarity matrix and the test-statistic retrieved is a multivariate analogue to Fisher's F-ratio as well as an R2 value (Adonis R2)."],"metadata":{"id":"mjaTwY3Tkde9"},"id":"mjaTwY3Tkde9"},{"cell_type":"code","execution_count":null,"id":"b3f215de","metadata":{"id":"b3f215de"},"outputs":[],"source":["adonres <- adonis2(dist_matrix  ~ md_Stats[,at_int],permutations = 999, distance='bray')\n","rownames(adonres)[1] <- colnames(md_Stats)[at_int]\n","adonres"]},{"cell_type":"markdown","source":["The PERMANOVA test result tells us that, for P < 0.05, there is a significant difference amng the different conditions within the given attribute (significant variation expressed as the percentage of  Adonis R2)."],"metadata":{"id":"dtA0FEaRkyuA"},"id":"dtA0FEaRkyuA"},{"cell_type":"code","execution_count":null,"id":"4f52d72d","metadata":{"id":"4f52d72d"},"outputs":[],"source":["PCoA_plot + labs(subtitle = paste0(\"p=\",round(adonres$'Pr(>F)'[1],4),', ' ,\"adonis-R2=\",round(adonres$'R2'[1],4)))"]},{"cell_type":"code","execution_count":null,"id":"4de93365","metadata":{"id":"4de93365"},"outputs":[],"source":["ggsave(PCoA_plot,filename=\"MDS_plot.svg\", width = 10, height = 8)"]}],"metadata":{"colab":{"collapsed_sections":[],"name":"20220719_CMFI_Seminar_Data_Cleanup_Prelim_Stats.ipynb","provenance":[{"file_id":"https://github.com/Functional-Metabolomics-Lab/CMFI_Seminar_Multivariate_Statistics/blob/main/20220719_CMFI_Seminar_Data_Cleanup_Prelim_Stats.ipynb","timestamp":1658162793428},{"file_id":"https://github.com/Functional-Metabolomics-Lab/CMFI_Seminar_Multivariate_Statistics/blob/main/20220719_CMFI_Seminar_Data_Cleanup_Prelim_Stats.ipynb","timestamp":1658154516086},{"file_id":"https://github.com/Functional-Metabolomics-Lab/CMFI_Seminar_Multivariate_Statistics/blob/main/0220716_CMFI_Seminar_Data_Cleanup.ipynb","timestamp":1658152498637},{"file_id":"https://github.com/abzer005/Summer-School_Functional-Metabolomics/blob/main/20220715_Data_CleanUp_SummerSchool_Metabolomics.ipynb","timestamp":1657987503674}]},"kernelspec":{"display_name":"R","language":"R","name":"ir"},"language_info":{"codemirror_mode":"r","file_extension":".r","mimetype":"text/x-r-source","name":"R","pygments_lexer":"r","version":"3.6.1"}},"nbformat":4,"nbformat_minor":5}