Rlegaspi Week 8

From LMU BioDB 2015
Jump to: navigation, search

Sample Microarray Analysis Vibrio cholerae

Statistical Analysis

  • The data from the Merrell et al. (2002) paper was accessed from this page at the Stanford Microarray Database.
  • The Log2 of R/G Normalized Ratio (Median) has been copied from the raw data files downloaded from the Stanford Microarray Database.
  • Downloaded the Merrell_Compiled_Raw_Data_Vibrio.xls file to my flashdrive.
    • Saved a copy of the file with a different filename that included my initials and the date of download - "Merrell_Compiled_Raw_Data_Vibrio_RARL_20151015.xls"

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

To scale and center the data (between chip normalization) I performed the following operations:

  • Inserted a new Worksheet into my Excel file, and named it "scaled_centered".
  • Went back to the "compiled_raw_data" worksheet, Select All and Copy. Went to my new "scaled_centered" worksheet, clicked on the upper, left-hand cell (cell A1) and Paste.
  • Inserted two rows in between the top row of headers and the first data row.
  • In cell A2, typed "Average" and in cell A3, typed "StdDev".
  • I computed the Average log ratio for each chip (each column of data). In cell B2, typed the following equation:
=AVERAGE(B4:B5224)
and pressed "Enter". Excel is computing the average value of the cells specified in the range given inside the parentheses. Instead of typing the cell designations, I clicked on the beginning cell, scrolled down to the bottom of the worksheet, and shift-click on the ending cell.
  • I then computed the Standard Deviation of the log ratios on each chip (each column of data). In cell B3, typed the following equation:
=STDEV(B4:B5224)
and pressed "Enter". I copied these two equations (cells B2 and B3) and pasted them into the empty cells in the rest of the columns. Excel automatically changed the equation to match the cell designations for those columns.
  • Since I computed the average and standard deviation of the log ratios for each chip, I began to actually do the scaling and centering based on these values.
  • I copied the column headings for all of my data columns and then pasted them to the right of the last data column so that I have a second set of headers above blank columns of cells. I edited the names of the columns so that they now read: A1_scaled_centered, A2_scaled_centered, etc.
  • In cell N4, I typed 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? - If we were to not include the dollar sign symbols in front of the "2" and "3" then equations that would be copied and pasted in the entire column would not use the values we needed which are found in cell B2 and B3. For instance, without the dollar sign symbols the equation for the next cell would become =(B5-B3/B4), the wrong values are being used for what we want.
  • Copied and pasted this equation into the entire column.
  • Copied and pasted the scaling and centering equation for each of the columns of data with the "_scaled_centered" column header.

Perform statistical analysis on the ratios

  • Inserted a new worksheet and named it "statistics".
  • Went back to the "scaling_centering" worksheet and copied the first column ("ID").
  • Pasted the data into the first column of my new "statistics" worksheet.
  • Went back to the "scaling_centering" worksheet and copied the columns that are designated "_scaled_centered".
  • Went to my 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.
  • Deleted Rows 2 and 3 where it says "Average" and "StDev" so that my data rows with gene IDs are immediately below the header row 1.
  • Went to a new column on the right of my worksheet. Typed the header "Avg_LogFC_A", "Avg_LogFC_B", and "Avg_LogFC_C" into the top cell of the next three columns.
  • Computed the average log fold change for the replicates for each patient by typing the equation:
=AVERAGE(B2:E2)
into cell N2. Copied this equation and pasted it into the rest of the column.
  • Created the equation for patients B and C and pasted it into their respective columns.
  • I computed the average of the averages. Typed the header "Avg_LogFC_all" into the first cell in the next empty column. Created the equation that will compute the average of the three previous averages I calculated and pasted it into this entire column.
  • Inserted a new column next to the "Avg_LogFC_all" column that I computed in the previous step. Labeled the column "Tstat". This computed a T statistic that tells us whether the scaled and centered average log ratio is significantly different than 0 (no change). Entered the equation: 
=AVERAGE(N2:P2)/(STDEV(N2:P2)/SQRT(number of replicates))
(NOTE: in this case the number of replicates is 3.) Copied the equation and pasted it into all rows in that column.
  • Labeled the top cell in the next column "Pvalue". In the cell below the label, entered 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. Copied the equation and pasted it into all rows in that column.

Calculate the Bonferroni p value Correction

  • Now I performed adjustments to the p value to correct for the multiple testing problem. Labeled the next two columns to the right with the same label, Bonferroni_Pvalue.
  • Typed the equation =S2*5221, Upon completion of this single computation, I used the trick to copy the formula throughout the column.
  • Replaced 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: =IF(T2>1,1,T2). Used the trick to copy the formula throughout the column.

Calculate the Benjamini & Hochberg p value Correction

  • Inserted a new worksheet named "B-H_Pvalue".
  • Copied and pasted the "ID" column from my previous worksheet into the first column of the new worksheet.
  • Inserted a new column on the very left and named it "MasterIndex".
    • Typed a "1" in cell A2 and a "2" in cell A3.
    • Selected both cells. Hovered my mouse over the bottom-right corner of the selection until it made a thin black + sign. Double-click'ed on the + sign to fill the entire column with a series of numbers from 1 to 5221 (the number of genes on the microarray).
  • For the following, use Paste special > Paste values. Copied my unadjusted p values from my previous worksheet and pasted it into Column C.
  • Selected all of columns A, B, and C. Sort by ascending values on Column C. Clicked the sort button from A to Z on the toolbar, in the window that appeared, sorted by column C, smallest to largest.
  • Typed the header "Rank" in cell D1. I created a series of numbers in ascending order from 1 to 5221 in this column. This is the p value rank, smallest to largest. Typed "1" into cell D2 and "2" into cell D3. Selected both cells D2 and D3. Double-click'ed on the plus sign on the lower right-hand corner of my selection to fill the column with a series of numbers from 1 to 5221.
  • I then calculated the Benjamini and Hochberg p value correction. Typed B-H_Pvalue in cell E1. Typed the following formula in cell E2: =(C2*5221)/D2 and pressed enter. Copied that equation to the entire column.
  • Typed "B-H_Pvalue" into cell F1.
  • Typed the following formula into cell F2: =IF(E2>1,1,E2) and pressed enter. Copied that equation to the entire column.
  • Selected columns A through F. Sorted them by my MasterIndex in Column A in ascending order.
  • Copied column F and used Paste special > Paste values to paste it into the next column on the right of my "statistics" sheet.

Prepare file for GenMAPP

  • Inserted a new worksheet and name it "forGenMAPP".
  • Went back to the "statistics" worksheet and Selected All and Copied.
  • Went to my new sheet and clicked on cell A1 and selected Paste Special, clicked on the Values radio button, and clicked OK. From here, I formatted this worksheet for import into GenMAPP.
  • Selected Columns B through Q (all the fold changes). Selected the menu item Format > Cells. Under the number tab, selected 2 decimal places. Clicked OK.
  • Selected all the columns containing p values. Selected the menu item Format > Cells. Under the number tab, selected 4 decimal places. Click OK.
  • Deleted the left-most Bonferroni p value column, preserving the one that showed the result of my "if" statement.
  • Inserted a column to the right of the "ID" column. Typed the header "SystemCode" into the top cell of this column. Filled the entire column (each cell) with the letter "N".
  • Selected the menu item File > Save As, and chose "Text (Tab-delimited) (*.txt)" from the file type drop-down menu. My new *.txt file was now ready for import into GenMAPP.
    • Uploaded both the .xls and .txt files that I created to the class wiki. My two files can be found here: .xls file and .txt file

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).

  • Open your spreadsheet and go 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.
  • 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.
    • How many genes have p value < 0.05? and what is the percentage (out of 5221)?
    • What about p < 0.01? and what is the percentage (out of 5221)?
    • What about p < 0.001? and what is the percentage (out of 5221)?
    • What about p < 0.0001? and what is the percentage (out of 5221)?
  • 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 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/.) 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:
    • How many genes are p < 0.05 for the Bonferroni-corrected p value? and what is the percentage (out of 5221)?
    • How many genes are p < 0.05 for the Benjamini and Hochberg-corrected p value? and what is the percentage (out of 5221)?
  • 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 "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.
    • 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 %)
    • 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 %)
    • What about an average log fold change of > 0.25 and p < 0.05? (and %)
    • 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.)
  • 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

  • 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

Preparation of GenMAPP

Launched GenMAPP, I made sure that the correct Gene Database (.gdb) is loaded.

  • Looked in the lower left-hand corner of the window to see which Gene Database has been selected.
  • I needed to change the Gene Database, so I selected Data > Choose Gene Database.
  • For the exercise in class on Thursday (10/22/2015), I needed to download the appropriate Vibrio cholerae Gene Database:
  • I clicked on the link for the Gene Database to which I had been assigned, downloaded the file, and saved it into the folder C:\GenMAPP 2 Data\Gene Databases, and extracted it.

GenMAPP Expression Dataset Manager Procedure

  • Launched the GenMAPP Program. Checked to make sure the correct Gene Database is loaded.
    • 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.
  • Selected the Data menu from the main Drafting Board window and chose Expression Dataset Manager from the drop-down list. The Expression Dataset Manager window opened.
  • Selected New Dataset from the Expression Datasets menu. Selected the tab-delimited text file that I formatted for GenMAPP .txt file created from Excel file in the procedure above from the file dialog box that appeared.
    • The Data Type Specification window appeared. 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; however, the Vibrio data we have been working with does not have any text (character) data in it.
  • Allowed the Expression Dataset Manager to convert my 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 - .gex file created after data conversion.
    • 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. 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: In my case, the total number of errors was 772 and all the error images indicated the following: Gene not found in OrderedLocusNames or any related system.
      • Error Comparison with HW Partner: It was likely that I had a different number of errors than my partner who was using the more update version of the Vibrio cholerae Gene Database. In comparison, my partner had 121 errors and I had 772 errors. Given the message Gene not found in OrderedLocusNames or any related system. Because the Gene Database was updated a year later, I would expect to be an update of the number of genes identified over time; therefore, Josh Kuroda's Gene Database allowed for the identification of more genes within his data meaning less errors.
      • Uploaded my exceptions file: EX.txt to my wiki page. File can be found here: Ron Legaspi's Exceptions File
  • Customized the new Expression Dataset by creating new Color Sets which contained 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]. I used "Avg_LogFC_all" for the Vibrio dataset I just created.
    • Activatde the Criteria Builder by clicking the New button.
    • Entered a name for the criterion in the Label in Legend field.
    • Chose a color for the criterion by left-clicking on the Color box. Chose a color from the Color window that appears and clicked OK.
    • Stated the criterion for color-coding a gene in the Criterion field.
  • After completing a new criterion, add the criterion entry (label, criterion, and color) to the Criteria List by clicking the Add button.
    • For the Vibrio dataset, I created 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.
  • Saved 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.
  • Exited 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.
  • Uploaded my .gex file. The .gex file can be found here: Ron Legaspi's .gex file

MAPPFinder Procedure

Note: My partner and I both did the same criterion, "Increased."

  • Launched 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\.)
Under construction
NEEDS EDITTING STARTING HERE
  • 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.
Top 10 Gene Ontology terms:
  1. Localization
  2. Transporter Activity
  3. Amino Acid Metabolic Process
  4. Cellular Nitrogen Compound Metabolic Process
  5. Cellular Amine Metabolic Process
  6. Cellular Amino Acid and Derivative Metabolic Process
  7. Biopolymer Biosynthetic Process
  8. Macromolecule Biosynthetic Process
  9. Cellular Macromolecule Metabolic Process
  10. Macromolecule 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.
  • 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 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.
  • 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 for part 1 and 2 this week.

List of Files to Upload

It may be easier to zip all of these files together and then upload them as a single zipped file, rather than zipping and uploading individually (for filetypes not allowed by OpenWetware).

  1. Your exceptions file when you imported your data into GenMAPP: .EX.txt
  2. Your Expression Dataset file: .gex
  3. Your GO results file: XXX-CriterionX-GO.txt
  4. Your GO results saved as an Excel spreadsheet with filters applied: .xls
  5. The MAPP you looked at: .mapp
  6. The MAPPFinder GO mappings file: .gmf
Retrieved from "https://xmlpipedb.lmucs.io/biodb/fall2015/index.php?title=Rlegaspi_Week_8&oldid=3326"