BIOL478/S14:Microarray Data Analysis

Laboratory 10: Microarray Data Analysis Tuesday, April 29, and Thursday, May 1

• Please complete the CURE post-course survey.

Background

This is a list of steps required to analyze DNA microarray data.

1. Quantitate the fluorescence signal in each spot
2. Calculate the ratio of red/green fluorescence
3. Log transform the ratios
4. Normalize the ratios on each microarray slide
   • Steps 1-4 are performed by the GenePix Pro software.
5. Normalize the ratios for a set of slides in an experiment
   • Step 5 was performed for you by Katrina using a script written in the R statistical programming language
   • You will perform the following steps:
6. Perform statistical analysis on the ratios
7. Compare individual genes with known data
   • Steps 6-7 are performed in Microsoft Excel
8. Pattern finding algorithms (clustering)
   • We will use software called STEM for clustering
9. Map onto biological pathways
   • We will use software called GenMAPP for mapping onto biological pathways

We will combine the data that you obtained in the Δswi4 cold shock and recovery experiment for flasks 3 and 4 with data obtained by the Spring 2013 class for flasks 1 and 2.

You will download an Excel spreadsheet that contains the normalized log fold changes from MyLMUConnect <insert link> for use in the following protocol.

Experimental Design

• On the spreadsheet, each row contains the data for one gene (one spot on the microarray). The first column (labeled "MasterIndex") numbers the rows in the spreadsheet so that we can match the data from different experiments together later. The second column (labeled "ID") contains the gene identifier from the Saccharomyces Genome Database. Each subsequent column contains the normalized log₂ ratio of the red/green fluorescence from each microarray hybridized in the experiment (steps 1-5 above having been done for you by the scanner software and R).
• Each of the column headings from the data begin with the experiment name ("dSWI4" for the Δswi4 data). "LogFC" stands for "Log₂ Fold Change" which is the Log₂ red/green ratio. The timepoints are designated as "t" followed by a number in minutes. Replicates are numbered as "-1", "-2", etc. after the timepoint, indicating which flask the cells were harvested from.
• The timepoints are t30, t60 (cold shock at 13°C) and t90 and t120 (cold shock at 13°C followed by 30 or 60 minutes of recovery at 30°C).
• There are four replicates for each timepoint.

Perform statistical analysis on the ratios

We are going to perform this step on the normalized data that you downloaded from MyLMUConnect.

• Insert a new worksheet into your Excel spreadsheet and name it "statistics".
• Go back to the "normalized" worksheet, Select All and Copy. Go to your new "statistics" worksheet, click on the upper, left-hand cell (cell A1) and 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.
  • There may be some non-numerical values in some of the cells in your worksheet. This is due to errors created when Excel tries to compute an equation on a cell that has no data. We need to go through and remove these error messages before going on to the next step.
  • Scan through your spreadsheet to find an example of the error message. Then go to the Edit menu and Select Replace. A window will open, type the text you are replacing in the "Find what:" field. In the "Replace with:" field, enter a single space character. Click on the button "Replace All" and record the number of replacements in your lab notebook.
• Go to the empty columns to the right on your worksheet. Create new column headings in the top cells to label the average log fold changes that you will compute. Name them with the pattern <dSWI4>_<AvgLogFC>_<tx> where you use the appropriate text within the <> and where x is the time. For example, "dSWI4_AvgLogFC_t30".
• Compute the average log fold change for the replicates for each timepoint by typing the equation:

=AVERAGE(range of cells in the row for that timepoint)

into the second cell below the column heading. For example, your equation might read

=AVERAGE(C2:F2)

Copy this equation and paste it into the rest of the column.

• Create the equation for the rest of the timepoints and paste it into their respective columns. Note that you can save yourself some time by completing the first equation for all of the averages and then copy and paste all the columns at once.
• Go to the empty columns to the right on your worksheet. Create new column headings in the top cells to label the T statistic that you will compute. Name them with the pattern <dSWI4>_<Tstat>_<tx> where you use the appropriate text within the <> and where x is the time. For example, "dSWI4_Tstat_t30". You will now compute a T statistic that tells you whether the normalized average log fold change is significantly different than 0 (no change in expression). Enter the equation into the second cell below the column heading:

=AVERAGE(range of cells)/(STDEV(range of cells)/SQRT(number of replicates))

For example, your equation might read:

=AVERAGE(C2:F2)/(STDEV(C2:F2)/SQRT(4))

(NOTE: in this case the number of replicates is 4. Be careful that you are using the correct number of parentheses.) Copy the equation and paste it into all rows in that column. Create the equation for the rest of the timepoints and paste it into their respective columns. Note that you can save yourself some time by completing the first equation for all of the T statistics and then copy and paste all the columns at once.

• Go to the empty columns to the right on your worksheet. Create new column headings in the top cells to label the P value that you will compute. Name them with the pattern <dSWI4>_<Pval>_<tx> where you use the appropriate text within the <> and where x is the time. For example, "dSWI4_Pval_t30". In the cell below the label, enter the equation:

=TDIST(ABS(cell containing T statistic),degrees of freedom,2)

For example, your equation might read:

=TDIST(ABS(AE2),3,2)

The number of degrees of freedom is the number of replicates minus one, so in our case there are 3 degrees of freedom. Copy the equation and paste it into all rows in that column.

• Insert a new worksheet and name it "final".
• Go back to the "statistics" worksheet and Select All and Copy.
• Go to your new sheet and click on cell A1 and select Paste Special, click on the Values radio button, and click OK. This is your final worksheet from which we will perform biological analysis of the data.
• Select all of the columns containing Fold Changes. Select the menu item Format > Cells. Under the number tab, select 2 decimal places. Click OK.
• Select all of the columns containing T statistics or P values. Select the menu item Format > Cells. Under the number tab, select 4 decimal places. Click OK.
• Upload the .xls file that you have just created to LionShare. Give Dr. Dahlquist (username kdahlqui) permission to download your file. Send an e-mail to Dr. Dahlquist with the link to the file.

Sanity Check: Number of genes significantly changed

Before we move on to the biological analysis of the data, 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 have significant changes in gene expression at different p value cut-offs.

• Open your spreadsheet and go to the "final" worksheet.
• 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.
• Click on the drop-down arrow on one of your "Pval" columns. Select "Custom". In the window that appears, set a criterion that will filter your data so that the P value has to be less than 0.05.
  • How many genes have p value < 0.05?
  • What about p < 0.01?
  • What about p < 0.001?
  • What about p < 0.0001?
    • Answer these questions for each timepoint in your dataset.
• 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.
• We have just performed 6189 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 309 times. 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.
  • There is a simple correction that can be made to the p values to increase the stringency called the Bonferroni correction. To perform this correction, multiply the p value by the number of statistical tests performed (in our case 6189) and see whether any of the p values are still less than 0.05.
    • Perform this correction and determine whether and how many of the genes are still significantly changed at p < 0.05 after the Bonferroni correction.
• The "AvgLogFC" tells us the magnitude 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. For the timepoint that had the greatest number of genes significantly changed at p < 0.05, answer the following:
  • Keeping the "Pval" filter at p < 0.05, filter the "AvgLogFC" column to show all genes with an average log fold change greater than zero. How many meet these two criteria?
  • Keeping the "Pval" filter at p < 0.05, filter the "AvgLogFC" column to show all genes with an average log fold change less than zero. How many meet these two criteria?
  • Keeping the "Pval" filter at p < 0.05, How many have an average log fold change of > 0.25 and p < 0.05?
  • How many have an average log fold change of < -0.25 and p < 0.05? (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.)
• 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.
• The expression of the gene NSR1 (ID: YGR159C)is known to be induced by cold shock. Find NSR1 in your dataset. Is it's expression significantly changed at any timepoint? Record the average fold change and p value for NSR1 for each timepoint in your dataset.
• Which gene has the smallest p value in your dataset (at any timepoint)? You can find this by sorting your data based on p value (but be careful that you don't cause a mismatch in the rows of your data!) Look up the function of this gene at the Saccharomyces Genome Database and describe it in your report. Why do you think the cell is changing this gene's expression upon cold shock?

Clustering and Gene Ontology Analysis with STEM

1. Begin by downloading and extracting the STEM software. Click here to go to the STEM web site.
   • Click on the download link, register, and download the stem.zip file to your Desktop.
   • Unzip the file. In Seaver 120, you can right click on the file icon and select the menu item 7-zip > Extract Here.
   • This will create a folder called stem. Inside the folder, double-click on the stem.cmd to launch the STEM program.
2. Prepare your microarray data file for loading into STEM.
   • Using the Excel spreadsheet from the analysis above, insert a new worksheet and name it "stem".
   • Select all of the data from your "final" worksheet and paste it into your "stem" worksheet.
     • Your leftmost column should have the column header "MasterIndex". Rename this column to "SPOT". Column B should be named "ID". Rename this column to "Gene Symbol".
     • Rename the data columns with just the time and units (for example, 30m, 60m, etc.).
     • We are going to delete all of the genes that do not have a significant change in gene expression at at least one timepoint. Do the following:
       1. Turn on filtering.
       2. For each timepoint, set a filter on the "Pval" column to > 0.05 (that's greater than). What this will do is select all of the genes that have a p value > 0.05 at each timepoint, which means conversely that the genes with p < 0.05 at at least on timepoint will not be selected.
       3. Record the number of genes that result from this filter in your lab notebook. Select all of the rows with genes meeting the filter and delete those rows. Then remove the filter.
     • Delete all of the data columns EXCEPT for the AvgLogFC columns for each timepoint.
     • Save your work. Then use Save As to save this spreadsheet as Text (Tab-delimited) (*.txt). Click OK to the warnings and close your file.
3. Running STEM
   1. In section 1 (Expression Data Info) of the the main STEM interface window, click on the Browse... button to navigate to and select your file.
      • Click on the radio button No normalization/add 0.
      • Check the box next to Spot IDs included in the data file.
   2. In section 2 (Gene Info) of the main STEM interface window, select Saccharomyces cerevisiae (SGD), from the drop-down menu for Gene Annotation Source. Select No cross references, from the Cross Reference Source drop-down menu. Select No Gene Locations from the Gene Location Source drop-down menu.
   3. In section 3 (Options) of the main STEM interface window, make sure that the Clustering Method says "STEM Clustering Method" and do not change the defaults for Maximum Number of Model Profiles or Maximum Unit Change in Model Profiles between Time Points.
   4. In section 4 (Execute) click on the yellow Execute button to run STEM.
4. Viewing and Saving STEM Results
   1. A new window will open called "All STEM Profiles (1)". Each box corresponds to a model expression profile. Colored profiles have a statistically significant number of genes assigned; they are arranged in order from most to least significant p value. Profiles with the same color belong to the same cluster of profiles. The number in each box is simply an ID number for the profile.
      • Click on the button that says "Interface Options...". At the bottom of the Interface Options window that appears below where it says "X-axis scale should be:", click on the radio button that says "Based on real time". Then close the Interface Options window.
      • Take a screenshot of this window (on a PC, simultaneously press the Alt and PrintScreen buttons to save the view in the active window to the clipboard) and paste it into a new PowerPoint presentation to save your figures.
   2. Click on each of the profiles to open a window showing a more detailed plot containing all of the genes in that profile.
      • Take a screenshot of each of the individual profile windows and save the images in your PowerPoint presentation.
      • At the bottom of each profile window, there are two yellow buttons "Profile Gene Table" and "Profile GO Table". For each of the profiles, click on the "Profile Gene Table" button to see the list of genes belonging to the profile. In the window that appears, click on the "Save Table" button and save the file to your desktop. Make your filename descriptive of the contents, e.g. "dSWI4_profile#_genelist.txt", where you replace the number symbol with the actual profile number.
        • Upload this file to LionShare and provide a link to Dr. Dahlquist.
      • For each of the profiles, click on the "Profile GO Table" to see the list of Gene Ontology terms belonging to the profile. In the window that appears, click on the "Save Table" button and save the file to your desktop. Make your filename descriptive of the contents, e.g. "dSWI4_profile#_GOlist.txt", where you replace the number symbol with the actual profile number. At this point you have saved all of the primary data from the STEM software and it's time to interpret the results!
        • Upload this file to LionShare and provide a link to Dr. Dahlquist.
5. Analyzing and Interpreting STEM Results
   1. Select one of the profiles you saved in the previous step for further intepretation of the data. We suggest that you choose one that has a pattern of up- or down-regulated genes at the early (first three) timepoints. Answer the following:
      • Why did you select this profile? In other words, why was it interesting to you?
      • How many genes belong to this profile?
      • How many genes were expected to belong to this profile?
      • What is the p value for the enrichment of genes in this profile? Bear in mind that in the protocol above, you computed p values to determine whether each individual gene had a significant change in gene expression at each time point. This p value determines whether the number of genes that show this particular expression profile across the time points is significantly more than expected.
      • Open the GO list file you saved for this profile in Excel. This list shows all of the Gene Ontology terms that are associated with genes that fit this profile. Select the third row and then choose from the menu Data > Filter > Autofilter. Filter on the "p-value" column to show only GO terms that have a p value of < 0.05. How many GO terms are associated with this profile at p < 0.05? The GO list also has a column called "Corrected p-value". This correction is needed because the software has performed thousands of significance tests. Filter on the "Corrected p-value" column to show only GO terms that have a corrected p value of < 0.05. How many GO terms are associated with this profile with a corrected p value < 0.05?
      • Select 10 Gene Ontology terms from your filtered list (either p < 0.05 or corrected p < 0.05). Look up the definitions for each of the terms at http://geneontology.org. Write a paragraph that describes the biological interpretation of these GO terms. In other words, why does the cell react to cold shock by changing the expression of genes associated with these GO terms?

Map onto Biological Pathways using GenMAPP

Each time you launch GenMAPP, you need to make sure that the correct Gene Database (.gdb) is loaded.

• Look in the lower left-hand corner of the window to see which Gene Database has been selected.
• 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, if the yeast Gene Database is not present on your computer, you will need to download it. Click this link to download the yeast Gene Database.
• Save the file, Sc-Std_20060526.gdb, to the folder C:\GenMAPP 2 Data\Gene Databases.

GenMAPP Expression Dataset Manager Procedure

• Launch the GenMAPP Program. Check to make sure the correct Gene Database is loaded.
• 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.
• 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.
• 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, check the box next to the field (column) name.
  • The yeast data we have been working with does not have any text (character) data in it.
• Allow the Expression Dataset Manager to convert your data.
  • 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.
  • 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.
• 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. You will have one Color Set per time point.
  • 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_" for the the appropriate time point.
  • 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. 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.
  • For the yeast dataset, you will create two criterion for each Color Set. "Increased" will be [Avg_LogFC_] > 0.25 AND [Pvalue] < 0.05 and "Decreased will be [Avg_LogFC_] < -0.25 AND [Pvalue] < 0.05. Make sure that the increased and decreased average log fold change values match the timepoint of the Color Set.
  • You may continue to add criteria to the Color Set by using the previous steps.
    • 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.
• 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.
• 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.
• Upload your .gex file to Lionshare and share it with Dr. Dahlquist. E-mail the link to the file to Dr. Dahlquist.