Difference between revisions of "Nanguiano Week 8"

From LMU BioDB 2015
Jump to: navigation, search
(Conclusion: changed header)
(Files as of 10/22/15: added GO and MAPPFinder files)
 
(31 intermediate revisions by the same user not shown)
Line 1: Line 1:
== Overview of Microarray Data Analysis ==
 
 
This is a list of steps required to analyze DNA microarray data.
 
 
# Quantitate the fluorescence signal in each spot in the microarray image.
 
#* Typically performed by the scanner software, although third party software packages do exist.
 
#* The image of the microarray slide and this quantitation are considered the "raw-est" form of the data.
 
#* Ideally, this type of raw data would be made publicly available upon publication. 
 
#* In practice, the image data is usually not made available because the raw image file of one slide could be up to 100 MB in size.
 
#* Also, some journals do not require data deposition as a requirement for publication, so often published data are not actually available anywhere for download.
 
#* Microarray data is not centrally located on the web.  Some major sources are:
 
#** [http://www.ncbi.nlm.nih.gov/geo/ NCBI GEO]
 
#** [http://www.ebi.ac.uk/microarray-as/ae/ EBI ArrayExpress]
 
#** [http://smd.princeton.edu/ Stanford Microarray Database (now hosted by Princeton)]
 
#** [http://puma.princeton.edu/ PUMAdb (Princeton Microarray Database)]
 
#** In addition, microarray data can sometimes be found as supplementary information with a journal article or on an investigator's own web site.
 
# Calculate the ratio of red/green fluorescence
 
# Log(base 2) transform the ratios
 
# Normalize the log ratios on each microarray slide
 
# Normalize the log ratios for a set of slides in an experiment
 
# Perform statistical analysis on the log ratios
 
# Compare individual genes with known data
 
# Look for patterns (expression profiles) in the data (many programs are available to do this; we are going to skip this step)
 
# Perform Gene Ontology term enrichment analysis (we will use MAPPFinder for this)
 
# Map onto biological pathways (we will use GenMAPP for this)
 
 
In this week's exercise, we will do steps 5-7 (part 1, using Microsoft Excel) and 9-10 (part 2, using GenMAPP & MAPPFinder).
 
 
 
== Statistical Analysis of ''Vibrio cholerae'' Microarray Data (Part 1) ==
 
== Statistical Analysis of ''Vibrio cholerae'' Microarray Data (Part 1) ==
  
* We will begin this analysis in class on Thursday, October 15.
+
=== Files as of 10/22/15 ===
* The detailed instructions for the microarray data analysis we will carry out can be found on the [http://www.openwetware.org/wiki/BIOL398-01/S10:Sample_Microarray_Analysis_Vibrio_cholerae Sample Microarray Analysis for ''Vibrio cholerae'' page] hosted by [http://www.openwetware.org OpenWetWare.org].  
+
* [[Media:Merrell_Compiled_Raw_Data_Vibrio_NA_20151015.xls | Current Spreadsheet]]
 +
<!-- [[File:Merrell_Compiled_Raw_Data_Vibrio_NA_20151015.xls]] -->
 +
* [[Media:Merrell_Compiled_Raw_Data_Vibrio_NA_20151015_GenMAPP.txt | Text file for GenMAPP]]
 +
<!-- [[File:Merrell_Compiled_Raw_Data_Vibrio_NA_20151015_GenMAPP.txt]] -->
 +
* [[Media:Merrell_Compiled_Raw_Data_Vibrio_NA_20151015_GenMAPP_Exceptions.EX.txt | GenMAPP Errors File]]
 +
<!-- [[File:Merrell_Compiled_Raw_Data_Vibrio_NA_20151015_GenMAPP_Exceptions.EX.txt]] -->
 +
* [[Media:Merrell_Compiled_Raw_Data_Vibrio_NA_20151015_Expression_ColorSet.gex | Expression Color Set]]
 +
<!-- [[File:Merrell_Compiled_Raw_Data_Vibrio_NA_20151015_Expression_ColorSet.gex]] -->
 +
* [[Media:GenMAPP-Results-Criterion1-GO_NA_Decreased.txt | MAPPFinder Gene Ontology Results]]
 +
<!-- [[File:GenMAPP-Results-Criterion1-GO_NA_Decreased.txt]] -->
 +
* [[Media:Phosphoribosylformylglycinamidine_synthase_activity_NA.mapp | Phosphoribosylformylglycinamidine Synthase Activity MAPP]]
 +
<!-- [[File:Phosphoribosylformylglycinamidine_synthase_activity_NA.mapp]] -->
  
* '''The below instructions are from the [http://www.openwetware.org/wiki/BIOL398-01/S10:Sample_Microarray_Analysis_Vibrio_cholerae Sample Microarray Analysis for ''Vibrio cholerae'' page] hosted by [http://www.openwetware.org OpenWetWare.org].'''
 
 
=== Normalize the log ratios for the set of slides in the experiment ===
 
=== Normalize the log ratios for the set of slides in the experiment ===
  
To scale and center the data (between chip normalization) perform the following operations:
+
* '''The below instructions are from the [http://www.openwetware.org/wiki/BIOL398-01/S10:Sample_Microarray_Analysis_Vibrio_cholerae Sample Microarray Analysis for ''Vibrio cholerae'' page] hosted by [http://www.openwetware.org OpenWetWare.org].'''
  
* Insert a new Worksheet into your Excel file, and name it "scaled_centered".
+
* First, I downloaded the data file from [http://www.openwetware.org/wiki/Media:Merrell_Compiled_Raw_Data_Vibrio.xls Open Wet Ware].
* Go back to the "compiled_raw_data" worksheet, Select All and CopyGo to your new "scaled_centered" worksheet, click on the upper, left-hand cell (cell A1) and Paste.
+
* To scale and center the data (between chip normalization), I performed the following operations:
* Insert two rows in between the top row of headers and the first data row.
+
** I inserted a new Worksheet into the Excel file, and name it "scaled_centered".
* In cell A2, type "Average" and in cell A3, type "StdDev".
+
** In the "compiled_raw_data" worksheet, I pressed <code>Command + A</code>to select everything, and then copied it allI pasted the information in the "scaled_centered" sheet.
* You will now compute the Average log ratio for each chip (each column of data). In cell B2, type the following equation:
+
** I Inserted two rows in between the column headers and the first row of data.
=AVERAGE(B4:B5224)
+
** I typed "Average" in cell A2 and "StdDev" in cell A3.  
: and press "Enter".  Excel is computing the average value of the cells specified in the range given inside the parentheses.  Instead of typing the cell designations, you can click on the beginning cell, scroll down to the bottom of the worksheet, and shift-click on the ending cell.
+
** Now I needed to compute the Average Log Ration for each column. In cell B2, I typed the following equation: <pre> =AVERAGE(B4:B5224) </pre>
* You will now compute the Standard Deviation of the log ratios on each chip (each column of data).  In cell B3, type the following equation:
+
** To compute the Standard Deviation, I typed the following equation in cell B3: <pre> =STDEV(B4:B5224) </pre>
=STDEV(B4:B5224)
+
** To apply the formula to all of the columns, I dragged the small box on the bottom right-hand corner across each column to apply it to all columns.  
: and press "Enter". 
+
**Next, I copied the column headings and pasted them to the right of C4to create a second set of headers.  I edited the names of the columns so they read:  A1_scaled_centered, A2_scaled_centered, etc.
* Excel will now do some work for you.  Copy these two equations (cells B2 and B3) and paste them into the empty cells in the rest of the columns.  Excel will automatically change the equation to match the cell designations for those columns.
+
** In cell N4 under the column titled "A1_scaled_centered", I typed the following equation: <pre> =(B4-B$2)/B$3 </pre>
* You have now computed the average and standard deviation of the log ratios for each chip.  Now we will actually do the scaling and centering based on these values.
+
** I clicked cell N4, then double-clicked the small box on the bottom right when the cursor turned into a plus so that it applied to every row in the column.  
* Copy the column headings for all of your data columns and then paste them to the right of the last data column so that you have a second set of headers above blank colums of cellsEdit the names of the columns so that they now read:  A1_scaled_centered, A2_scaled_centered, etc.
+
** I dragged the formula from N4 across the remaining columns to apply it to all of them, then applied the formula to all of the rows of the columns.
* In cell N4, type the following equation:
+
=(B4-B$2)/B$3
+
: In this case, we want the data in cell B4 to have the average subtracted from it (cell B2) and be divided by the standard deviation (cell B3).  We use the dollar sign symbols in front of the "2" and "3" to tell Excel to always reference that row in the equation, even though we will paste it for the entire column of 5221 genes. '''''Why is this important?'''''
+
* Copy and paste this equation into the entire column.  One easy way to do this is to click on the original cell with your equation and position your cursor at the bottom right corner. You should see your cursor change to a thin black plus sign (not a chubby white one). When it does, double click, and the formula will magically be copied to the entire column of genes.
+
* Copy and paste the scaling and centering equation for each of the columns of data with the "_scaled_centered" column header.  Be sure that your equation is correct for the column you are calculating.
+
  
 
=== Perform statistical analysis on the ratios ===
 
=== Perform statistical analysis on the ratios ===
  
We are going to perform this step on the scaled and centered data you produced in the previous step.
+
* To perform the statistical analysis, I performed the following steps:
 
+
** First, I inserted a new worksheet and named it "statistics".
* Insert a new worksheet and name it "statistics".
+
** I copied the first column from the "scaled_centered" worksheet (the column titled ID), then pasted the column into the new sheet.
* Go back to the "scaling_centering" worksheet and copy the first column ("ID").
+
** Back in the "scaled_centered" sheet, I copied all of the columns with "_scaled_centered" in the titles and pasted them using "Edit -> Paste Special -> Values" in the "statistics" sheet, starting at cell B1.
* Paste the data into the first column of your new "statistics" worksheet.
+
** Then, I deleted rows 2 and 3 from "statistics" (Average/StdDev).
* Go back to the "scaling_centering" worksheet and copy the columns that are designated "_scaled_centered".
+
** To the right of my "_scaled_centered" columns, I created three columns, titled "Avg_LogFC_A", "Avg_LogFC_B", and "Avg_LogFC_C".
* Go to your new worksheet and click on the B1 cell.  Select "Paste Special" from the Edit menu.  A window will open: click on the radio button for "Values" and click OK.  This will paste the numerical result into your new worksheet instead of the equation which must make calculations on the fly.
+
** To compute the average log fold change for the replicates, I typed in cell N2 under "Avg_LogFC_A" the equation: <pre> =AVERAGE(B2:E2) </pre>
* Delete Rows 2 and 3 where it says "Average" and "StDev" so that your data rows with gene IDs are immediately below the header row 1.
+
** I applied the formula to the entire column, as well as to the next two columns.
* Go to a new column on the right of your worksheet.  Type the header "Avg_LogFC_A", "Avg_LogFC_B", and "Avg_LogFC_C" into the top cell of the next three columns.
+
** Now I needed to computer the average of the averages.  In the column after "Avg_LogFC_C", I created a new column titled "Avg_LogFC_all". Under the column header, I typed the equation: <pre> =AVERAGE(N2:P2) </pre>
* Compute the average log fold change for the replicates for each patient by typing the equation:
+
** After applying the formula to the whole column, I created a new column called "Tstat" next to the "Avg_LogFC_all". Under this column header, I entered the following equation: <pre> =AVERAGE(N2:P2)/(STDEV(N2:P2)/SQRT(3)) </pre>
=AVERAGE(B2:E2)
+
** I then applied it to all rows in the column.
: into cell N2.  Copy this equation and paste it into the rest of the column
+
** Next to the "Tstat" column, I labeled the next column "Pvalue".  In the cell below the label, I entered the equation:  <pre> =TDIST(ABS(R2),2,2) </pre> <!-- when sorting by p<0.05, should be 948 values-->
* Create the equation for patients B and C and paste it into their respective columns.
+
** Then, I applied the formula to all of the rows in the column.
* Now you will compute the average of the averages.  Type the header "Avg_LogFC_all" into the first cell in the next empty column.  Create the equation that will compute the average of the three previous averages you calculated and paste it into this entire column.
+
* Insert a new column next to the "Avg_LogFC_all" column that you computed in the previous step. Label the column "Tstat".  This will compute a T statistic that tells us whether the scaled and centered average log ratio is significantly different than 0 (no change).  Enter the equation:
+
=AVERAGE(N2:P2)/(STDEV(N2:P2)/SQRT(number of replicates))
+
: (NOTE: in this case the number of replicates is 3.  Be careful that you are using the correct number of parentheses.)  Copy the equation and paste it into all rows in that column.
+
* Label the top cell in the next column "Pvalue".  In the cell below the label, enter the equation:   
+
=TDIST(ABS(R2),degrees of freedom,2)
+
The number of degrees of freedom is the number of replicates minus one, so in our case there are 2 degrees of freedom.  Copy the equation and paste it into all rows in that column.
+
  
 
=== Calculate the Bonferroni p value Correction ===
 
=== Calculate the Bonferroni p value Correction ===
  
* Now we will perform adjustments to the p value to correct for the [https://xkcd.com/882/ multiple testing problem].  Label the next two columns to the right with the same label, Bonferroni_Pvalue.
+
* I labeled the next two columns to the right of the "Pvalue" column with the same label, "Bonferroni_Pvalue".
* Type the equation <code>=<(Pvalue>*5221</code>, Upon completion of this single computation, use the trick to copy the formula throughout the column.
+
* Under the first "Bonferroni_Pvalue" column, I typed the equation: <pre> =S2*5221 </pre>  
* Replace any corrected p value that is greater than 1 by the number 1 by typing the following formula into the first cell below the second Bonferroni_Pvalue header: <code>=IF(s2>1,1,s2)</code>.  Use the trick to copy the formula throughout the column.
+
* Then, I applied the formula to the whole column. In order to replace any corrected P value that was greater than 1 with the number 1, I typed the following formula under the header of the second "Bonferonni_Pvalue" column: <pre>=IF(T2>1,1,T2)</pre>.   
 +
* I applied the formula to the whole column.
  
 
=== Calculate the Benjamini & Hochberg p value Correction ===
 
=== Calculate the Benjamini & Hochberg p value Correction ===
  
* Insert a new worksheet named "B-H_Pvalue".
+
* To calculate the B&H p value correction, I began by creating a new worksheet and naming it  "B-H_Pvalue".
* Copy and paste the "ID" column from your previous worksheet into the first column of the new worksheet.  
+
* I coped and pasted the "ID" column from the "statistics" worksheet into the second column.  
* Insert a new column on the very left and name it "MasterIndex". We will create a numerical index of genes so that we can always sort them back into the same order.
+
* I named the first column "MasterIndex". I typed a 1 in cell A2 and a 2 in cell A3, selected both cells, then double-clicked the box in the right-hand corner to fill the entire column.
** Type a "1" in cell A2 and a "2" in cell A3.
+
* Using "Paste Special -> Values", I copied the Pvalue column from "statistics" into column C.
** Select both cells. Hover your mouse over the bottom-right corner of the selection until it makes a thin black + sign. Double-click on the + sign to fill the entire column with a series of numbers from 1 to 5221 (the number of genes on the microarray).  
+
* I selected columns A, B, and C, and sorted by ascending order on Column C. To do this, I selected Data -> Sort. I sorted by column C, smallest to largest.
* For the following, use Paste special > Paste values.  Copy your unadjusted p values from your previous worksheet and paste it into Column C.
+
* In cell D1, I typed the header "Rank". This column will contain the p value rank, from smallest to largest.  I repeated the same process that I used for the "MasterIndex" column in order to fill the column with numbers from 1 to 5,221.  
* Select all of columns A, B, and C. Sort by ascending values on Column C. Click the sort button from A to Z on the toolbar, in the window that appears, sort by column C, smallest to largest.
+
* In column E1, I typed "B-H_Pvalue". In cell E2, I typed the formula: <pre>=(C2*5221)/D2</pre>  
* Type the header "Rank" in cell D1. We will create a series of numbers in ascending order from 1 to 5221 in this column.  This is the p value rank, smallest to largest.  Type "1" into cell D2 and "2" into cell D3. Select both cells D2 and D3. Double-click on the plus sign on the lower right-hand corner of your selection to fill the column with a series of numbers from 1 to 5221.
+
* I applied the formula to the whole column.  
* Now you can calculate the Benjamini and Hochberg p value correction. Type B-H_Pvalue in cell E1. Type the following formula in cell E2: <code>=(C2*5221)/D2</code> and press enter. Copy that equation to the entire column.
+
* In cell F1, I typed "B-H_Pvalue".
* Type "B-H_Pvalue" into cell F1.  
+
* In cell F2, I typed the formula: <pre>=IF(E2>1,1,E2)</pre>
* Type the following formula into cell F2: <code>=IF(E2>1,1,E2)</code> and press enter. Copy that equation to the entire column.  
+
* I copied that equation into the entire column.
* Select columns A through F.  Now sort them by your MasterIndex in Column A in ascending order.
+
* I sorted the columns in ascending order by the MasterIndex column, then copied column F (the adjusted B-H_Pvalues) and used Paste -> Paste Special to paste it in the statistics sheet to the right of the Bonferroni_Pvalue.
* Copy column F and use Paste special > Paste values to paste it into the next column on the right of your "statistics" sheet.
+
  
 
=== Prepare file for GenMAPP ===
 
=== Prepare file for GenMAPP ===
  
* Insert a new worksheet and name it "forGenMAPP".
+
* I inserted a new worksheet and name it "forGenMAPP".
* Go back to the "statistics" worksheet and Select All and Copy.
+
* I copied everything from the statistics worksheet, and pasted it into the new worksheet using Paste -> Paste Special -> Values.
* Go to your new sheet and click on cell A1 and select Paste Special, click on the Values radio button, and click OK.  We will now format this worksheet for import into GenMAPP.
+
* Next, I began formatting the spreadsheet for GenMAPP.
* Select Columns B through Q (all the fold changes).  Select the menu item Format > Cells. Under the number tab, select 2 decimal places.  Click OK.
+
** I selected Columns B through Q (all the fold changes), then selected Format -> Cells. Under the number tab, I selected 2 decimal places.   
* Select all the columns containing p values. Select the menu item Format > Cells.  Under the number tab, select 4 decimal places. Click OK.
+
** Next, I selected Columns S through U (p values). I selected Format -> Cells, and under the number tab, selected 4 decimal places.
* Delete the left-most Bonferroni p value column, preserving the one that shows the result of your "if" statement.
+
** I deleted column T, which was the left-most Bonferroni P value column, preserving the one that had the results of the If statement.
* Insert a column to the right of the "ID" column. Type the header "SystemCode" into the top cell of this column.  Fill the entire column (each cell) with the letter "N".
+
** To the right of the "ID" column, I inserted a column. I typed the header "SystemCode" into the top cell, and replaced the entire column with the letter "N".
* Select the menu item File > Save As, and choose "Text (Tab-delimited) (*.txt)" from the file type drop-down menu.  Excel will make you click through a couple of warnings because it doesn't like you going all independent and choosing a different file type than the native .xls.  This is OK.  Your new *.txt file is now ready for import into GenMAPP.  But before we do that, we want to know a few things about our data as shown in the next section.
+
** Lastly, I saved the file as a "Text (Tab-delimited) (*.txt)" file.
** Upload both the .xls and .txt files that you have just created to your journal page in the class wiki.  Make sure that your file name is distinct from your other classmates so that nobody overwrites anyone else's file.
+
  
 
=== Sanity Check: Number of genes significantly changed ===
 
=== Sanity Check: Number of genes significantly changed ===
  
Before we move on to the GenMAPP/MAPPFinder analysis, we want to perform a sanity check to make sure that we performed our data analysis correctly.  We are going to find out the number of genes that are significantly changed at various p value cut-offs and also compare our data analysis with the published results of Merrell et al. (2002).
+
Now, I wanted to perform a sanity check to make sure I performed the analysis correctly.
  
* Open your spreadsheet and go to the "forGenMAPP" tab.
+
* To begin, I opened the spreadsheet and went to the "forGenMAPP" tab.  
* Click on cell A1 and select the menu item Data > Filter > Autofilter.  Little drop-down arrows should appear at the top of each column.  This will enable us to filter the data according to criteria we set.
+
* I selected cell A1, then selected Data -> Filter to activate the filter for each tab.  
* Click on the drop-down arrow on your "Pvalue" column.  Select "Custom".  In the window that appears, set a criterion that will filter your data so that the Pvalue has to be less than 0.05.
+
* On the "Pvalue" column, I selected the dropdown arrow and applied the following filters, recording the number of values that remained after filtering, as well as what percentage that is of the genes lie within the range.
** '''''How many genes have p value < 0.05? and what is the percentage (out of 5221)?'''''
+
** p < 0.05: 948 (18.15%)
** '''''What about p < 0.01? and what is the percentage (out of 5221)?'''''
+
** p < 0.01: 235 (4.50%)
** '''''What about p < 0.001? and what is the percentage (out of 5221)?'''''
+
** p < 0.001: 24 (0.45%)
** '''''What about p < 0.0001? and what is the percentage (out of 5221)?'''''
+
** p < 0.0001: 2 (0.03%)
* When we use a p value cut-off of p < 0.05, what we are saying is that you would have seen a gene expression change that deviates this far from zero less than 5% of the time.
+
* After clearing the filter on the "Pvalue" column, I then moved to filter the Bonferroni and B-H P value columns. I performed the following filters, recording the same values as before (number and percentage).
* We have just performed 5221 T tests for significance.  Another way to state what we are seeing with p < 0.05 is that we would expect to see this magnitude of a gene expression change in about 5% of our T tests, or 261 times. (Test your understanding: [http://xkcd.com/882/ http://xkcd.com/882/].) Since we have more than 261 genes that pass this cut off, we know that some genes are significantly changed.  However, we don't know ''which'' ones.  To apply a more stringent criterion to our p values, we performed the Bonferroni and Benjamini and Hochberg corrections to these unadjusted p values.  The Bonferroni correction is very stringent.  The Benjamini-Hochberg correction is less stringent. To see this relationship, filter your data to determine the following:
+
** Bonferroni p < 0.05: 0 (0%)
** '''''How many genes are p < 0.05 for the Bonferroni-corrected p value? and what is the percentage (out of 5221)?'''''
+
** B&H p < 0.05: 0 (0%)
** '''''How many genes are p < 0.05 for the Benjamini and Hochberg-corrected p value? and what is the percentage (out of 5221)?'''''
+
* Next, I filtered the "Avg_LogFC_all" column and recorded the number and percent of each search.
* In summary, the p value cut-off should not be thought of as some magical number at which data becomes "significant".  Instead, it is a moveable confidence level.  If we want to be very confident of our data, use a small p value cut-off.  If we are OK with being less confident about a gene expression change and want to include more genes in our analysis, we can use a larger p value cut-off. 
+
** p < 0.05 && Avg_LogFC_all > 0: 352 (6.74%)
* The "Avg_LogFC_all" tells us the size of the gene expression change and in which direction.  Positive values are increases relative to the control; negative values are decreases relative to the control.
+
** p < 0.05 && Avg_LogFC_all < 0: 596 (11.41%)
** Keeping the (unadjusted) "Pvalue" filter at p < 0.05, filter the "Avg_LogFC_all" column to show all genes with an average log fold change greater than zero. '''''How many are there? (and %)'''''
+
** p < 0.05 && Avg_LogFC_all > 0.25: 339 (6.49%)
** Keeping the (unadjusted) "Pvalue" filter at p < 0.05, filter the "Avg_LogFC_all" column to show all genes with an average log fold change less than zero. '''''How many are there? (and %)'''''
+
** p < 0.05 && Avg_LogFC_all < -0.25: 579 (11.08%)
** '''''What about an average log fold change of > 0.25 and p < 0.05? (and %)'''''
+
* '''''What criteria did [http://www.nature.com/nature/journal/v417/n6889/full/nature00778.html Merrell et al. (2002)] use to determine a significant gene expression change?  How does it compare to our method?'''''
** '''''Or an average log fold change of < -0.25 and p < 0.05? (and %)'''''  (These are more realistic values for the fold change cut-offs because it represents about a 20% fold change which is about the level of detection of this technology.)
+
** Merrel et al. (2002) used a two-class SAM (statistical analysis for microarrays) analysis to determine a significant expression change. Genes with a twofold change in expression or more in all patents were considered significantly changed. This seems to indicate that they used the change in expression to help determine statistical significance, whereas we used p values to determine significance.
* In summary, the p value cut-off should not be thought of as some magical number at which data becomes "significant".  Instead, it is a moveable confidence level.  If we want to be very confident of our data, use a small p value cut-off.  If we are OK with being less confident about a gene expression change and want to include more genes in our analysis, we can use a larger p value cut-off.  For the GenMAPP analysis below, we will use the fold change cut-off of greater than 0.25 or less than -0.25 and the unadjusted p value cut off of p < 0.05 for our analysis because we want to include several hundred genes in our analysis.
+
* '''''What criteria did Merrell et al. (2002) use to determine a significant gene expression change?  How does it compare to our method?'''''
+
  
 
=== Sanity Check:  Compare individual genes with known data ===
 
=== Sanity Check:  Compare individual genes with known data ===
Line 143: Line 110:
 
== MAPPFinder Analysis of ''Vibrio cholerae'' Microarray Data (Part 2) ==
 
== MAPPFinder Analysis of ''Vibrio cholerae'' Microarray Data (Part 2) ==
  
* We will begin this analysis in class on Tuesday, October 20.
+
=== Map Onto Biological Pathways (GenMAPP & MAPPFinder) ===
* The detailed instructions can be found on the [http://www.openwetware.org/wiki/BIOL367/F10:GenMAPP_and_MAPPFinder_Protocols GenMAPP and MAPPFinder Protocols page] hosted by [http://www.openwetware.org OpenWetWare.org].
+
 
+
 
* '''The below instructions are from the [http://www.openwetware.org/wiki/BIOL367/F10:GenMAPP_and_MAPPFinder_Protocols GenMAPP and MAPPFinder Protocols page] hosted by [http://www.openwetware.org OpenWetWare.org].'''
 
* '''The below instructions are from the [http://www.openwetware.org/wiki/BIOL367/F10:GenMAPP_and_MAPPFinder_Protocols GenMAPP and MAPPFinder Protocols page] hosted by [http://www.openwetware.org OpenWetWare.org].'''
=== Map Onto Biological Pathways (GenMAPP & MAPPFinder) ===
 
  
Each time you launch GenMAPP, you need to make sure that the correct Gene Database (.gdb) is loaded.
+
Each time I launched GenMAPP, I needed to make sure I had the correct Gene Database (.gdb) loaded. The database I was using was [http://sourceforge.net/projects/xmlpipedb/files/V.%20cholerae%20Gene%20Database/V.%20cholerae%2020101022/Vc-Std_External_20101022.zip/download Vc-Std_External_20101022.gdb], which I extracted prior to continuing.
* Look in the lower left-hand corner of the window to see which Gene Database has been selected.
+
* In the lower left hand corner of GenMAPP, it indicated that there was no database selected. I selected Data -> Choose Gene Database, then selected the .gdb file.
* If you need to change the Gene Database, select Data > Choose Gene Database.  Navigate to the directory C:\GenMAPP 2 Data\Gene Databases and choose the correct one for your species.
+
* For the exercise today, you will need to download the appropriate ''Vibrio cholerae'' Gene Database.
+
** Half of the class will use the Vc-Std_External_20090622.gdb Gene Database that was initially created by the Fall 2008 Biological Databases class.
+
*** To download this Gene Database, [http://sourceforge.net/projects/xmlpipedb/files/V.%20cholerae%20Gene%20Database/V.%20cholerae%2020090622/Vc-Std_External_20090622.zip/download follow '''''this link''''' to the XMLPipeDB SourceForge Download page].
+
** Half of the class will use a new Vc-Std_External_20101022.gdb Gene Database that was created by Drs. Dahlquist and Dionisio a year later.
+
*** To download this Gene Database, [http://sourceforge.net/projects/xmlpipedb/files/V.%20cholerae%20Gene%20Database/V.%20cholerae%2020101022/Vc-Std_External_20101022.zip/download follow '''''this link''''' to the XMLPipeDB SourceForge Download page].
+
** The members of a pair should each choose a different gene database.
+
* Click on the link for the Gene Database to which you have been assigned, download the file, and save it into the folder C:\GenMAPP 2 Data\Gene Databases, and extract it.
+
  
 
==== GenMAPP Expression Dataset Manager Procedure ====
 
==== GenMAPP Expression Dataset Manager Procedure ====
  
* Launch the GenMAPP Program.  Check to make sure the correct Gene Database is loaded.
+
* I selected Data -> Expression Dataset Manager, then chose Expression Dataset -> New Dataset. I selected the text file that I'd created (linked above) in the dialog box.  
** Look in the lower, left-hand corner of the main GenMAPP Drafting Board window to see the name of the Gene Database that is loaded.  If this is not the correct Gene Database or it says "No Gene Database", then go to the Data > Choose Gene Database menu item to select the Gene Database you need to perform the analysis.
+
** In the window that appeared, I just hit "ok" to continue.
** '''Remember, you and your partner are going to use ''different versions'' of the ''Vibrio cholerae'' Gene Database for this exercise.'''
+
** The file reported '''121 errors'''.
* Select the Data menu from the main Drafting Board window and choose Expression Dataset Manager from the drop-down list. The Expression Dataset Manager window will open.
+
*** In the EX.txt file, I opened it in excel and filtered by the "~Errors~" column. All 121 errors had the error message "Gene not found in OrderedLocusNames or any related system."
* Select New Dataset from the Expression Datasets menu. Select the tab-delimited text file that you formatted for GenMAPP (.txt) in the procedure above from the file dialog box that appears.
+
*** My partner, Emily Simso, had 772 errors. This is likely due to the names that were used for the genes. Since the file I used was newer, it's possible that the gene names were changed. Additionally, since the file I'm using was created by Dr. Dahlquist and Dr. Dionisio, it's possible that they knew the proper names and protocol to use to ensure the most successful files.
** You may need to download your .txt file from the wiki onto your Desktop if you have not already done so.
+
* Now, I needed to customize the color set for GenMAPP. Under the "Color Sets" "Name" header, I named it "Significant Genes" and selected "Avg_LogFC_All" under the Gene Value dropdown. I created two criterion that I was going to color:
* The Data Type Specification window will appear.  GenMAPP is expecting that you are providing numerical data.  If any of your columns has text (character) data, you would check the box next to the field (column) name.
+
** "Increased": [Avg_LogFC_all] > 0.25 AND [Pvalue] < 0.05
** ''The Vibrio data we have been working with does not have any text (character) data in it.''
+
** "Decreased": [Avg_LogFC_all] < -0.25 AND [Pvalue] < 0.05
* Allow the Expression Dataset Manager to convert your data.
+
** To create the criteria, I selected "New" under Criteria Builder. I entered the name, then selected a color. For "Increased", I used a magenta pink color, and for "Decreased", I used a cyan blue color. Then, I wrote in the criteria, selecting the first part of the expression under the "Columns" section, selecting the operation under the "Ops" column, and then typing in the number. I used the operation "AND" between the two pieces of the expression.
** This may take a few minutes depending on the size of the dataset and the computer’s memory and processor speed. When the process is complete, the converted dataset will be active in the Expression Dataset Manager window and the file will be saved in the same folder the raw data file was in, named the same except with a .gex extension; for example, MyExperiment.gex.
+
* Lastly, I saved the expression dataset by selecting Expression Datasets -> Save, then exited.
** A message may appear saying that the Expression Dataset Manager could not convert one or more lines of data. Lines that generate an error during the conversion of a raw data file are not added to the Expression Dataset. Instead, an exception file is created. The exception file is given the same name as your raw data file with .EX before the extension (e.g., MyExperiment.EX.txt). The exception file will contain all of your raw data, with the addition of a column named ~Error~. This column contains either error messages or, if the program finds no errors, a single space character.
+
*** '''Record the number of errors. For your journal assignment, open the .EX.txt file and use the Data > Filter > Autofilter function to determine what the errors were for the rows that were not converted.  Record this information in your individual journal page.'''
+
*** '''It is likely that you will have a different number of errors than your partner who is using a different version of the ''Vibrio cholerae'' Gene Database. Which of you has more errors?  Why do you think that is?  Record your answers in your journal page.'''
+
*** '''Upload your exceptions file:  <code>EX.txt</code> to your wiki page.
+
* Customize the new Expression Dataset by creating new Color Sets which contain the instructions to GenMAPP for displaying data on MAPPs.
+
** Color Sets contain the instructions to GenMAPP for displaying data from an Expression Dataset on MAPPs. Create a Color Set by filling in the following different fields in the Color Set area of the Expression Dataset Manager:  a name for the Color Set, the gene value, and the criteria that determine how a gene object is colored on the MAPP. Enter a name in the Color Set Name field that is 20 characters or fewer.
+
** The Gene Value is the data displayed next to the gene box on a MAPP. Select the column of data to be used as the Gene Value from the drop down list or select [none]. We will use "Avg_LogFC_all" for the Vibrio dataset you just created.
+
** Activate the Criteria Builder by clicking the New button.
+
** Enter a name for the criterion in the Label in Legend field.
+
** Choose a color for the criterion by left-clicking on the Color box. Choose a color from the Color window that appears and click OK.
+
** State the criterion for color-coding a gene in the Criterion field.
+
*** A criterion is stated with relationships such as "this column greater than this value" or "that column less than or equal to that value". Individual relationships can be combined using as many ANDs and ORs as needed. A typical relationship is
+
[ColumnName] RelationalOperator Value
+
::with the column name always enclosed in brackets and character values enclosed in single quotes. For example:
+
[Fold Change] >= 2
+
[p value] < 0.05
+
[Quality] = 'high'
+
::This is the equivalent to queries that you performed on the command line when working with the PostgreSQL movie database.  GenMAPP is using a graphical user interface (GUI) to help the user format the queries correctly.  The easiest and safest way to create criteria is by choosing items from the Columns and Ops (operators) lists shown in the Criteria Builder. The Columns list contains all of the column headings from your Expression Dataset. To choose a column from the list, click on the column heading. It will appear at the location of the cursor in the Criterion box. The Criteria Builder surrounds the column names with brackets.
+
  
::The Ops (operators) list contains the relational operators that may be used in the criteria: equals ( = )  greater than ( > ), less than ( < ), greater than or equal to ( >= ), less than or equal to ( <= ), is not equal to ( <> ). To choose an operator from the list, click on the symbol. It will appear at the location of the insertion bar (cursor) in the Criterion box. The Criteria Builder automatically surrounds the operators with spaces.
+
==== MAPPFinder Procedure ====
::The Ops list also contains the conjunctions AND and OR, which may be used to make compound criteria. For example:
+
[Fold Change] > 1.2 AND [p value] <= 0.05
+
::Parentheses control the order of evaluation. Anything in parentheses is evaluated first. Parentheses may be nested. For example:
+
[Control Average] = 100 AND ([Exp1 Average] > 100 OR [Exp2 Average] > 100)
+
::Column names may be used anywhere a value can, for example:
+
[Control Average] < [Experiment Average]
+
  
* After completing a new criterion, add the criterion entry (label, criterion, and color) to the Criteria List by clicking the Add button.
+
* Next, I needed to use MAPPFinder to look at '''decreased expression'''. I selected Tools -> MAPPFinder to launch the program.
** For the Vibrio dataset, you will create two criterion.  "Increased" will be [Avg_LogFC_all] > 0.25 AND [Pvalue] < 0.05 and "Decreased will be [Avg_LogFC_all] < -0.25 AND [Pvalue] < 0.05.
+
** I selected "Calculate New Results", then selected "Find File" and selected the ".gex" file that had been produced in the last step.  
** You may continue to add criteria to the Color Set by using the previous steps.
+
** I close the Color Set "Significant Genes", and filtered by "Decreased".
*** The buttons to the right of the list represent actions that can be performed on individual criteria. To modify a criterion label, color, or the criterion itself, first select the criterion in the list by left-clicking on it, and then click the Edit button. This puts the selected criterion into the Criteria Builder to be modified. Click the Save button to save changes to the modified criterion; click the Add button to add it to the list as a separate criterion. To remove a criterion from the list, left-click on the criterion to select it, and then click on the Delete button. The order of Criteria in the list has significance to GenMAPP. When applying an Expression Dataset and Color Set to a MAPP, GenMAPP examines the expression data for a particular gene object and applies the color for the first criterion in the list that is true. Therefore, it is imperative that when criteria overlap the user put the most important or least inclusive criteria in the list first. To change the order of the criteria in the list, left-click on the criterion to select it and then click the Move Up or Move Down buttons. No criteria met and Not found are always the last two positions in the list.
+
** I checked the boxes "Gene Ontology" and "Click here to calculate p values", then hit "Browse", and selected where I would save the new file, as well as gave it a name.
* Save the entire Expression Dataset by selecting Save from the Expression Dataset menu. Changes made to a Color Set are not saved until you do this.
+
** Then, I hit "Run MAPPFinder".
* Exit the Expression Dataset Manager to view the Color Sets on a MAPP. Choose Exit from the Expression Dataset menu or click the close box in the upper right hand corner of the window.
+
* I selected "Show Ranked List" after the Gene Ontology Browser opened. The top 10 GO terms were listed below:
* '''Upload your .gex file to your journal entry page for later retrieval.'''
+
*# glucose catabolic process
 +
*# hexose catabolic process
 +
*# glycolysis
 +
*# monosaccharide catabolic process
 +
*# cytoplasm
 +
*# alcohol catabolic process
 +
*# cellular carbohydrate catabolic process
 +
*# glucose metabolic process
 +
*# protein folding
 +
*# hexose metabolic process
 +
** '''Compare your list with your partner who used a different version of the Gene Database. Are your terms the same or different?  Why do you think that is?  Record your answer in your individual journal entry.'''
 +
 
 +
* In the MAPPFinder Browser window, I selected "Collapse Tree". To search for a gene, I typed the gene ID into the "Search for a specific Gene ID" input box, then selected "OrderedLocusNames" in the "Select Gene ID type" dropdown. I searched for the following terms, and recorded the corresponding GO terms:
 +
** VC0028:
 +
*** GO Terms: branched chain family amino acid biosynthetic process, cellular amino acid biosynthetic process, metabolic process, metal ion binding, iron-sulfur cluster binding, 4 iron, 4 sulfur cluster binding, catalytic activity, lyase activity, dihydroxy-acid dehydratase activity
 +
** VC0941:
 +
*** GO Terms: glycine metabolic process, L-serine metabolic process, one-carbon metabolic process, cytoplasm, pyridoxal phosphate binding, catalytic activity, transferase activity, glycine hydroxymethyltransferase activity,
 +
** VC0869:
 +
*** GO Terms: glutamine metabolic process, purine nucleotide biosynthetic process, 'de novo' IMP biosynthetic process, cytoplasm, nucleotide binding, ATP binding, catalytic activity, ligase activity, phosphoribosylformyglycinamidine synthase activity
 +
** VC0051:
 +
*** GO Terms: purine nucleotide biosynthetic process, 'de novo' IMP biosynthetic process, nucleotide binding, ATP binding, catalytic activity, lyase activity, carboxy-lyase activity,  phosphoribosylaminoimidazole carboxylase activity
 +
** VC0647:
 +
*** GO Terms: mRNA catabolic process, RNA processing, cytoplasm, mitochondion, RNA binding, 3'-5'-exoribonuclease activity, transferase activity, nucleotidyltransferase activity, polyribonucleotide nuclotidyltransferase activity
 +
** VC0468:
 +
*** GO Terms: glutathione biosynthetic process, metal ion binding, nucloetide binding, ATP binding, catalytic activity, ligase activity, glutathione synthase activity
 +
** VC2350:
 +
*** GO Terms: deoxyribonucleotide catabolic process, metabolic process, cytoplasm, catalytic activity, lyase activity, deoxyribose-phosphase aldolase activity
 +
** VCA0583:
 +
*** GO Terms: transport, outer membrane-bounded periplasmic space, transporter activity
 +
* The GO Terms I found were not the same as the ones my partner found. <!-- why or why not? -->
 +
* Searching for VC0869 again, I double-clicked the GO term "phosphoribosylformyglycinamidine synthase activity". This opened a window showing what genes went under that category. The expression did change significantly, in that the expression was increased.  
 +
** I double-clicked the gene, named PUR4_VIBCH.  
 +
*** The function of this gene, according to [http://www.ncbi.nlm.nih.gov/gene?cmd=Retrieve&dopt=Graphics&list_uids=2614536 NCBI], the function of the gene is to "catalyze the formation of 2-(formamido)-N1-(5-phospho-D-ribosyl)acetamidine from N2-formyl-N1-(5-phospho-D-ribosyl)glycinamide and L-glutamine in purine biosynthesis"
  
==== MAPPFinder Procedure ====
 
  
'''''Note: You and your partner will both do the same criterion, either "Increased" or "Decreased", but your group does not need to do both "Increased" and "Decreased"  Sign up for the criterion you want on the group list ([[BIOL367/F10#Groups_3 | Fall 2010]] or [https://xmlpipedb.cs.lmu.edu/biodb/fall2013/index.php/Week_8#Groups Fall 2013]) so that we can make sure that as a class we are covering both criteria.'''''
 
  
* Launch the MAPPFinder program (or from within GenMAPP, select Tools > MAPPFinder).
 
* Make sure that the Gene Database for the correct species is loaded.  The name of the Gene Database appears at the bottom of the window.  If this is not the right one, go to File > Choose Gene Database and choose the correct one.  (The Gene Databases are stored in the folder C:\GenMAPP 2 Data\Gene Databases\.)
 
* Click on the button "Calculate New Results".
 
* Click on "Find File" and choose the your Expression Dataset file, for example, "MyDataset.gex", and click OK.
 
** MAPPFinder may have found it for you already if you already had it open in GenMAPP, in which case, you just need to click OK.
 
* Choose the Color Set and Criteria with which to filter the data.  Click on either the "Increased" and "Decreased" criteria in the right-hand box, depending on which one your group is doing.  (You could select both by holding down the Control key while clicking).
 
* Check the boxes next to "Gene Ontology" and "p value".
 
* Click the "Browse" button and create a meaningful filename for your results.
 
* Click "Run MAPPFinder".  The analysis will take several minutes.  It may look like the computer is stalled; be patient, it will eventually start running.
 
* When the results have been calculated, a Gene Ontology browser will open showing your results.  All of the Gene Ontology terms that have at least 3 genes measured and a p value of less than 0.05 will be highlighted yellow.  A term with a p value less than 0.05 is considered a "significant" result.  Browse through the tree to see your results.
 
* To see a list of the most significant Gene Ontology terms, click on the menu item "Show Ranked List". 
 
** '''List the top 10 Gene Ontology terms in your individual journal entry.'''
 
** '''Compare your list with your partner who used a different version of the Gene Database.  Are your terms the same or different?  Why do you think that is?  Record your answer in your individual journal entry.'''
 
* One of the things you can do in MAPPFinder is to find the Gene Ontology term(s) with which a particular gene is associated. First, in the main MAPPFinder Browser window, click on the button "Collapse the Tree". Then, you can search for the genes that were mentioned by Merrell et al. (2002), VC0028, VC0941, VC0869, VC0051, VC0647, VC0468, VC2350, and VCA0583. Type the identifier for one of these genes into the MAPPFinder browser gene ID search field. Choose "OrderedLocusNames" from the drop-down menu to the right of the search field.  Click on the GeneID Search button.  The GO term(s) that are associated with that gene will be highlighted in blue. '''List the GO terms associated with each of those genes in your individual journal.  (Note: they might not all be found.)  Are they the same as your partner who is using a different Gene Database?  Why or why not?'''
 
* Click on one of the GO terms that are associated with one of the genes you looked up in the previous step. A MAPP will open listing all of the genes (as boxes) associated with that GO term. The genes named within the map are based on the UniProt identification system. To match the gene of interest to its identification go to the [http://www.uniprot.org/ UniProt site] and type in your gene ID into the search bar. Moreover, the genes on the MAPP will be color-coded with the gene expression data from the microarray experiment.  '''List in your journal entry the name of the GO term you clicked on and whether the expression of the gene you were looking for changed significantly in the experiment.'''
 
** Double-click on the gene box.  This will open a Internet Explorer window called the "Backpage" for this gene.  This page has links to pages for this gene in the public databases.  '''Click on the links to find out the function of this gene and record your answer in your individual journal page.'''
 
** The MAPP that has just been created is stored in the directory, C:\GenMAPP 2 Data\MAPPs\VC GO.  '''Upload this file and link to it in your journal.'''
 
* In Windows, make a copy of your results (XXX-CriterionX-GO.txt) file. 
 
** "XXX" refers to the name you gave to your results file.
 
** "CriterionX" refers to either "Criterion0" or "Criterion1".  Since computers start counting at zero, "Criterion0" is the first criterion in the list you clicked on ("Increased" if you followed the directions) and "Criterion1" is the second criterion in the list you clicked on ("Decreased" if you followed the directions).
 
** '''Upload your results file to your journal page.'''
 
 
* Launch Microsoft Excel.  Open the copies of the .txt files in Excel (you will need to "Show all files" and click "Finish" to the wizard that will open your file).  This will show you the same data that you saw in the MAPPFinder Browser, but in tabular form.
 
* Launch Microsoft Excel.  Open the copies of the .txt files in Excel (you will need to "Show all files" and click "Finish" to the wizard that will open your file).  This will show you the same data that you saw in the MAPPFinder Browser, but in tabular form.
 
* Look at the top of the spreadsheet.  There are rows of information that give you the background information on how MAPPFinder made the calculations.  '''Compare this information with your partner who used a different version of the Vibrio Gene Database.  Which numbers are different?  Why are they different?  Record this information in your individual journal entry.'''
 
* Look at the top of the spreadsheet.  There are rows of information that give you the background information on how MAPPFinder made the calculations.  '''Compare this information with your partner who used a different version of the Vibrio Gene Database.  Which numbers are different?  Why are they different?  Record this information in your individual journal entry.'''
Line 251: Line 192:
  
 
* Write a paragraph that briefly summarizes and gives a scientific conclusion for the work that you did this week.
 
* Write a paragraph that briefly summarizes and gives a scientific conclusion for the work that you did this week.
 
=== Optional: downloading and installing the GenMAPP and MAPPFinder Software ===
 
 
* We will be using GenMAPP and MAPPFinder version 2.1 (http://genmapp.org).  This software is Windows-only and is already installed on the machines in the Seaver 120 computer lab.
 
** This version is now called "GenMAPP Classic" and can be downloaded [https://github.com/GenMAPPCS/genmapp/blob/master/GenMAPPv2Setup.exe?raw=true from this page].
 
** Follow the instructions in the installer.
 
** During installation, the installer will open a window called the GenMAPP Data Acquisition Tool.  It will not function because it cannot connect to the server.  This is OK, you will download your ''Vibrio cholerae'' Gene Database from the XMLPipeDB project at SourceForge.org.
 
*** Half of the class will use the Vc-Std_External_20090622.gdb Gene Database that was created by the Fall 2008 Biological Databases class.
 
**** To download this Gene Database, [http://sourceforge.net/projects/xmlpipedb/files/V.%20cholerae%20Gene%20Database/V.%20cholerae%2020090622/Vc-Std_External_20090622.zip/download follow '''''this link''''' to the XMLPipeDB SourceForge Download page].
 
*** Half of the class will use a more recent Vc-Std_External_20101022.gdb Gene Database that was created by Drs. Dahlquist and Dionisio in 2010.
 
**** To download this Gene Database, [http://sourceforge.net/projects/xmlpipedb/files/V.%20cholerae%20Gene%20Database/V.%20cholerae%2020101022/Vc-Std_External_20101022.zip/download follow '''''this link''''' to the XMLPipeDB SourceForge Download page].
 
*** The members of a pair should each choose a different gene database.
 
* Click on the link for the Gene Database to which you have been assigned, download the file, and save it into the folder C:\GenMAPP 2 Data\Gene Databases (if you accepted the default folders during the installation), and extract it.
 
  
 
== Links ==
 
== Links ==
 
{{Template:Nanguiano}}
 
{{Template:Nanguiano}}

Latest revision as of 22:58, 22 October 2015

Statistical Analysis of Vibrio cholerae Microarray Data (Part 1)

Files as of 10/22/15

Normalize the log ratios for the set of slides in the experiment

  • First, I downloaded the data file from Open Wet Ware.
  • To scale and center the data (between chip normalization), I performed the following operations:
    • I inserted a new Worksheet into the Excel file, and name it "scaled_centered".
    • In the "compiled_raw_data" worksheet, I pressed Command + Ato select everything, and then copied it all. I pasted the information in the "scaled_centered" sheet.
    • I Inserted two rows in between the column headers and the first row of data.
    • I typed "Average" in cell A2 and "StdDev" in cell A3.
    • Now I needed to compute the Average Log Ration for each column. In cell B2, I typed the following equation:
       =AVERAGE(B4:B5224) 
    • To compute the Standard Deviation, I typed the following equation in cell B3:
       =STDEV(B4:B5224) 
    • To apply the formula to all of the columns, I dragged the small box on the bottom right-hand corner across each column to apply it to all columns.
    • Next, I copied the column headings and pasted them to the right of C4to create a second set of headers. I edited the names of the columns so they read: A1_scaled_centered, A2_scaled_centered, etc.
    • In cell N4 under the column titled "A1_scaled_centered", I typed the following equation:
       =(B4-B$2)/B$3 
    • I clicked cell N4, then double-clicked the small box on the bottom right when the cursor turned into a plus so that it applied to every row in the column.
    • I dragged the formula from N4 across the remaining columns to apply it to all of them, then applied the formula to all of the rows of the columns.

Perform statistical analysis on the ratios

  • To perform the statistical analysis, I performed the following steps:
    • First, I inserted a new worksheet and named it "statistics".
    • I copied the first column from the "scaled_centered" worksheet (the column titled ID), then pasted the column into the new sheet.
    • Back in the "scaled_centered" sheet, I copied all of the columns with "_scaled_centered" in the titles and pasted them using "Edit -> Paste Special -> Values" in the "statistics" sheet, starting at cell B1.
    • Then, I deleted rows 2 and 3 from "statistics" (Average/StdDev).
    • To the right of my "_scaled_centered" columns, I created three columns, titled "Avg_LogFC_A", "Avg_LogFC_B", and "Avg_LogFC_C".
    • To compute the average log fold change for the replicates, I typed in cell N2 under "Avg_LogFC_A" the equation:
       =AVERAGE(B2:E2) 
    • I applied the formula to the entire column, as well as to the next two columns.
    • Now I needed to computer the average of the averages. In the column after "Avg_LogFC_C", I created a new column titled "Avg_LogFC_all". Under the column header, I typed the equation:
       =AVERAGE(N2:P2) 
    • After applying the formula to the whole column, I created a new column called "Tstat" next to the "Avg_LogFC_all". Under this column header, I entered the following equation:
       =AVERAGE(N2:P2)/(STDEV(N2:P2)/SQRT(3)) 
    • I then applied it to all rows in the column.
    • Next to the "Tstat" column, I labeled the next column "Pvalue". In the cell below the label, I entered the equation:
       =TDIST(ABS(R2),2,2) 
    • Then, I applied the formula to all of the rows in the column.

Calculate the Bonferroni p value Correction

  • I labeled the next two columns to the right of the "Pvalue" column with the same label, "Bonferroni_Pvalue".
  • Under the first "Bonferroni_Pvalue" column, I typed the equation:
     =S2*5221 
  • Then, I applied the formula to the whole column. In order to replace any corrected P value that was greater than 1 with the number 1, I typed the following formula under the header of the second "Bonferonni_Pvalue" column:
    =IF(T2>1,1,T2)
    .
  • I applied the formula to the whole column.

Calculate the Benjamini & Hochberg p value Correction

  • To calculate the B&H p value correction, I began by creating a new worksheet and naming it "B-H_Pvalue".
  • I coped and pasted the "ID" column from the "statistics" worksheet into the second column.
  • I named the first column "MasterIndex". I typed a 1 in cell A2 and a 2 in cell A3, selected both cells, then double-clicked the box in the right-hand corner to fill the entire column.
  • Using "Paste Special -> Values", I copied the Pvalue column from "statistics" into column C.
  • I selected columns A, B, and C, and sorted by ascending order on Column C. To do this, I selected Data -> Sort. I sorted by column C, smallest to largest.
  • In cell D1, I typed the header "Rank". This column will contain the p value rank, from smallest to largest. I repeated the same process that I used for the "MasterIndex" column in order to fill the column with numbers from 1 to 5,221.
  • In column E1, I typed "B-H_Pvalue". In cell E2, I typed the formula:
    =(C2*5221)/D2
  • I applied the formula to the whole column.
  • In cell F1, I typed "B-H_Pvalue".
  • In cell F2, I typed the formula:
    =IF(E2>1,1,E2)
  • I copied that equation into the entire column.
  • I sorted the columns in ascending order by the MasterIndex column, then copied column F (the adjusted B-H_Pvalues) and used Paste -> Paste Special to paste it in the statistics sheet to the right of the Bonferroni_Pvalue.

Prepare file for GenMAPP

  • I inserted a new worksheet and name it "forGenMAPP".
  • I copied everything from the statistics worksheet, and pasted it into the new worksheet using Paste -> Paste Special -> Values.
  • Next, I began formatting the spreadsheet for GenMAPP.
    • I selected Columns B through Q (all the fold changes), then selected Format -> Cells. Under the number tab, I selected 2 decimal places.
    • Next, I selected Columns S through U (p values). I selected Format -> Cells, and under the number tab, selected 4 decimal places.
    • I deleted column T, which was the left-most Bonferroni P value column, preserving the one that had the results of the If statement.
    • To the right of the "ID" column, I inserted a column. I typed the header "SystemCode" into the top cell, and replaced the entire column with the letter "N".
    • Lastly, I saved the file as a "Text (Tab-delimited) (*.txt)" file.

Sanity Check: Number of genes significantly changed

Now, I wanted to perform a sanity check to make sure I performed the analysis correctly.

  • To begin, I opened the spreadsheet and went to the "forGenMAPP" tab.
  • I selected cell A1, then selected Data -> Filter to activate the filter for each tab.
  • On the "Pvalue" column, I selected the dropdown arrow and applied the following filters, recording the number of values that remained after filtering, as well as what percentage that is of the genes lie within the range.
    • p < 0.05: 948 (18.15%)
    • p < 0.01: 235 (4.50%)
    • p < 0.001: 24 (0.45%)
    • p < 0.0001: 2 (0.03%)
  • After clearing the filter on the "Pvalue" column, I then moved to filter the Bonferroni and B-H P value columns. I performed the following filters, recording the same values as before (number and percentage).
    • Bonferroni p < 0.05: 0 (0%)
    • B&H p < 0.05: 0 (0%)
  • Next, I filtered the "Avg_LogFC_all" column and recorded the number and percent of each search.
    • p < 0.05 && Avg_LogFC_all > 0: 352 (6.74%)
    • p < 0.05 && Avg_LogFC_all < 0: 596 (11.41%)
    • p < 0.05 && Avg_LogFC_all > 0.25: 339 (6.49%)
    • p < 0.05 && Avg_LogFC_all < -0.25: 579 (11.08%)
  • What criteria did Merrell et al. (2002) use to determine a significant gene expression change? How does it compare to our method?
    • Merrel et al. (2002) used a two-class SAM (statistical analysis for microarrays) analysis to determine a significant expression change. Genes with a twofold change in expression or more in all patents were considered significantly changed. This seems to indicate that they used the change in expression to help determine statistical significance, whereas we used p values to determine significance.

Sanity Check: Compare individual genes with known data

  • Merrell et al. (2002) report that genes with IDs: VC0028, VC0941, VC0869, VC0051, VC0647, VC0468, VC2350, and VCA0583 were all significantly changed in their data. Look these genes up in your spreadsheet. What are their fold changes and p values? Are they significantly changed in our analysis?

MAPPFinder Analysis of Vibrio cholerae Microarray Data (Part 2)

Map Onto Biological Pathways (GenMAPP & MAPPFinder)

Each time I launched GenMAPP, I needed to make sure I had the correct Gene Database (.gdb) loaded. The database I was using was Vc-Std_External_20101022.gdb, which I extracted prior to continuing.

  • In the lower left hand corner of GenMAPP, it indicated that there was no database selected. I selected Data -> Choose Gene Database, then selected the .gdb file.

GenMAPP Expression Dataset Manager Procedure

  • I selected Data -> Expression Dataset Manager, then chose Expression Dataset -> New Dataset. I selected the text file that I'd created (linked above) in the dialog box.
    • In the window that appeared, I just hit "ok" to continue.
    • The file reported 121 errors.
      • In the EX.txt file, I opened it in excel and filtered by the "~Errors~" column. All 121 errors had the error message "Gene not found in OrderedLocusNames or any related system."
      • My partner, Emily Simso, had 772 errors. This is likely due to the names that were used for the genes. Since the file I used was newer, it's possible that the gene names were changed. Additionally, since the file I'm using was created by Dr. Dahlquist and Dr. Dionisio, it's possible that they knew the proper names and protocol to use to ensure the most successful files.
  • Now, I needed to customize the color set for GenMAPP. Under the "Color Sets" "Name" header, I named it "Significant Genes" and selected "Avg_LogFC_All" under the Gene Value dropdown. I created two criterion that I was going to color:
    • "Increased": [Avg_LogFC_all] > 0.25 AND [Pvalue] < 0.05
    • "Decreased": [Avg_LogFC_all] < -0.25 AND [Pvalue] < 0.05
    • To create the criteria, I selected "New" under Criteria Builder. I entered the name, then selected a color. For "Increased", I used a magenta pink color, and for "Decreased", I used a cyan blue color. Then, I wrote in the criteria, selecting the first part of the expression under the "Columns" section, selecting the operation under the "Ops" column, and then typing in the number. I used the operation "AND" between the two pieces of the expression.
  • Lastly, I saved the expression dataset by selecting Expression Datasets -> Save, then exited.

MAPPFinder Procedure

  • Next, I needed to use MAPPFinder to look at decreased expression. I selected Tools -> MAPPFinder to launch the program.
    • I selected "Calculate New Results", then selected "Find File" and selected the ".gex" file that had been produced in the last step.
    • I close the Color Set "Significant Genes", and filtered by "Decreased".
    • I checked the boxes "Gene Ontology" and "Click here to calculate p values", then hit "Browse", and selected where I would save the new file, as well as gave it a name.
    • Then, I hit "Run MAPPFinder".
  • I selected "Show Ranked List" after the Gene Ontology Browser opened. The top 10 GO terms were listed below:
    1. glucose catabolic process
    2. hexose catabolic process
    3. glycolysis
    4. monosaccharide catabolic process
    5. cytoplasm
    6. alcohol catabolic process
    7. cellular carbohydrate catabolic process
    8. glucose metabolic process
    9. protein folding
    10. hexose metabolic process
    • Compare your list with your partner who used a different version of the Gene Database. Are your terms the same or different? Why do you think that is? Record your answer in your individual journal entry.
  • In the MAPPFinder Browser window, I selected "Collapse Tree". To search for a gene, I typed the gene ID into the "Search for a specific Gene ID" input box, then selected "OrderedLocusNames" in the "Select Gene ID type" dropdown. I searched for the following terms, and recorded the corresponding GO terms:
    • VC0028:
      • GO Terms: branched chain family amino acid biosynthetic process, cellular amino acid biosynthetic process, metabolic process, metal ion binding, iron-sulfur cluster binding, 4 iron, 4 sulfur cluster binding, catalytic activity, lyase activity, dihydroxy-acid dehydratase activity
    • VC0941:
      • GO Terms: glycine metabolic process, L-serine metabolic process, one-carbon metabolic process, cytoplasm, pyridoxal phosphate binding, catalytic activity, transferase activity, glycine hydroxymethyltransferase activity,
    • VC0869:
      • GO Terms: glutamine metabolic process, purine nucleotide biosynthetic process, 'de novo' IMP biosynthetic process, cytoplasm, nucleotide binding, ATP binding, catalytic activity, ligase activity, phosphoribosylformyglycinamidine synthase activity
    • VC0051:
      • GO Terms: purine nucleotide biosynthetic process, 'de novo' IMP biosynthetic process, nucleotide binding, ATP binding, catalytic activity, lyase activity, carboxy-lyase activity, phosphoribosylaminoimidazole carboxylase activity
    • VC0647:
      • GO Terms: mRNA catabolic process, RNA processing, cytoplasm, mitochondion, RNA binding, 3'-5'-exoribonuclease activity, transferase activity, nucleotidyltransferase activity, polyribonucleotide nuclotidyltransferase activity
    • VC0468:
      • GO Terms: glutathione biosynthetic process, metal ion binding, nucloetide binding, ATP binding, catalytic activity, ligase activity, glutathione synthase activity
    • VC2350:
      • GO Terms: deoxyribonucleotide catabolic process, metabolic process, cytoplasm, catalytic activity, lyase activity, deoxyribose-phosphase aldolase activity
    • VCA0583:
      • GO Terms: transport, outer membrane-bounded periplasmic space, transporter activity
  • The GO Terms I found were not the same as the ones my partner found.
  • Searching for VC0869 again, I double-clicked the GO term "phosphoribosylformyglycinamidine synthase activity". This opened a window showing what genes went under that category. The expression did change significantly, in that the expression was increased.
    • I double-clicked the gene, named PUR4_VIBCH.
      • The function of this gene, according to NCBI, the function of the gene is to "catalyze the formation of 2-(formamido)-N1-(5-phospho-D-ribosyl)acetamidine from N2-formyl-N1-(5-phospho-D-ribosyl)glycinamide and L-glutamine in purine biosynthesis"


  • Launch Microsoft Excel. Open the copies of the .txt files in Excel (you will need to "Show all files" and click "Finish" to the wizard that will open your file). This will show you the same data that you saw in the MAPPFinder Browser, but in tabular form.
  • Look at the top of the spreadsheet. There are rows of information that give you the background information on how MAPPFinder made the calculations. Compare this information with your partner who used a different version of the Vibrio Gene Database. Which numbers are different? Why are they different? Record this information in your individual journal entry.
  • You will filter this list to show the top GO terms represented in your data for both the "Increased" and "Decreased" criteria. You will need to filter your list down to about 20 terms. Click on a cell in the row of headers for the data. Then go to the Data menu and click "Filter > Autofilter". Drop-down arrows will appear in the row of headers. You can now choose to filter the data. Click on the drop-down arrow for the column you wish to filter and choose "(Custom…)". A window will open giving you choices on how you want to filter. You must set these two filters:
Z Score (in column N) greater than 2
PermuteP (in column O) less than 0.05
You will use these two filters depending on the number of terms you have:
Number Changed (in column I) greater than or equal to 4 or 5 AND less than 100
Percent Changed (in column L) greater than or equal to 25-50%
  • Save your changes to an Excel spreadsheet. Select File > Save As and select Excel workbook (.xls) from the drop-down menu. Your filter settings won’t be saved in a .txt file.
  • Are any of your filtered GO terms closely related to one another, meaning are they a direct child or parent to another term in the list? You can judge this by comparing your spreadsheet with the MAPPFinder browser. Highlight the terms that fit this relationship with the same color in your Excel spreadsheet. Upload your .xls file to your journal page.
  • Interpret your results. Look up the definitions for any GO terms that are unfamiliar to you. The "official" definitions for GO terms can be found at http://www.geneontology.org. You can use one of the online biological dictionaries as a supplement, if needed. Write a paragraph relating the results of this GO analysis to the experiment performed (comparing laboratory-grown and patient-derived Vibrio cholerae. You need to give a biological interpretation of what do each of these GO terms in your filtered list have to to with the pathogenecity of the bacterium? You may consult with your partner on this, but your explanation on your individual journal page needs to be in your own words. This is where the real "brain power" comes in with interpreting DNA microarray data. Even experienced scientists struggle with this part. Use your creativity as a scientist to stretch your brain in this question.
  • There is one other file you need to save to your journal page. It has a .gmf extension and should be in the same fold as the .gex file that you created with the GenMAPP Expression Dataset Manager. You will need this file to re-open your results in MAPPFinder.

Conclusion

  • Write a paragraph that briefly summarizes and gives a scientific conclusion for the work that you did this week.

Links

Nicole Anguiano
BIOL 367, Fall 2015

Assignment Links
Individual Journals
Shared Journals