Difference between revisions of "Blitvak Week 8"

Revision as of 01:10, 27 October 2015

Steps taken for Parts 1-2 of the Microarray Data Analysis were edited from:

• Part I
• Part II

Statistical Analysis of Vibrio cholerae Microarray Data (Part 1)

• Merrell_Compiled_Raw_Data_Vibrio.xls was downloaded, saved to the desktop, and renamed with some additional information (initials and the date)

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

The following operations were performed in order to scale and center the microarray data:

• The renamed Excel file was opened and a new Worksheet was inserted with the name scaled_centered
• Everything on the compiled_raw_data worksheet was selected and copied over to scaled_centered (formatting was the same, starting from the left-hand cell, A1)
• Two new rows were inserted between the top row of headers and the first data row in scaled_centered
• In cell A2, Average was typed in; in A3, StdDev was typed in
• The Average log ratio for each chip was computed by typing =AVERAGE(B4:B5224) into cell B2 and pressing enter
• The Standard Deviation of the log ratios on each chip was computed by typing =STDEV(B4:B5224) into cell B3 and pressing enter
• The equations in B2 and B3 were copied and pasted into the empty cells in the rest of the columns (A2 to C4)
• The column headings for all of the data columns were copied and pasted to the right of the last data column; this new set of headers was edited so that they read: A1_scaled_centered, A2_scaled_centered, etc.
• The equation =(B4-B$2)/B$3 was typed into cell N4; the dollar sign symbols were used in front of the "2" and "3" in order to ensure that Excel will not change the reference to that row when that same equation is pasted down the entire column of 5221 genes (this is important because the average and standard deviation is the same for the entire row, and therefore, the reference must stay the same). This equation is the scaling and centering equation.
• The scaling and centering equation was copied and pasted down the entire A1_scaled_centered by clicking the original cell with the equation and double-clicking the bottom right corner of the cell (cursor should change to a black plus sign prior to double-clicking)
• The scaling and centering equation was put in each of the data columns with the _scaled_centered header (was copied and pasted down the entire columns)
• The equation for each column was checked to ensure that it was correct (ex. for A2_scaled_centered, the equation should be =(C4-C$2)/C$3)

Performing statistical analysis on the ratios

Initial statistical analysis

• A new worksheet was inserted with the name statistics
• The first ID column in the scaling_centering worksheet was copied and pasted into the first column of statistics
• The columns that are designated with _scaled_centered were copied and pasted into the new worksheet (starting from B1); for the pasting, "Paste Special" was required in order to just paste the numerical results into the new worksheet ("Values" was selected from the "Paste Special" window)
• Rows 2 and 3 (corresponding to Average and StDev) were deleted
• The headers Avg_LogFC_A, Avg_LogFC_B, and Avg_LogFC_C were typed into the next three empty columns to the right (immediately adjacent to the last _scaled_centered column)
• The average log fold change for the replicates for each patient was computed by typing =AVERAGE(B2:E2) into cell N2. The equation was copied and pasted down the entire column.
• The equation for the average log fold change was created for patients B and C; this equation was copied and pasted down their respective columns (for patient B the equation was =AVERAGE(F2:I2), for patient C the equation was =AVERAGE(J2:M2))
• In the first cell that corresponds to the next empty column, the header Avg_LogFC_all was typed
• An equation that can compute the average of the three previously calculated averages was created (=AVERAGE(N2:P2)); this equation was pasted into this entire column (Avg_LogFC_all)
• A new column was inserted next to Avg_LogFC_all. This column was given the label/header of Tstat (purpose of this column is to compute a T statistic that will inform whether
• The equation =AVERAGE(N2:P2)/(STDEV(N2:P2)/SQRT(3)) was entered into Tstat and pasted into all of the rows within that column
• The top cell in the next column was labeled with Pvalue; the equation =TDIST(ABS(R2),2,2) was entered in the cell below the label and copied and pasted into all of the rows in that column. The first "2" in that equation is the degrees of freedom (there are 2 degrees of freedom since the number of replicates, which is 3, minus 1 is 2)

Calculating the Bonferroni p-value Correction

• Adjustments to the p-value were performed with the purpose of correcting for the multiple testing problem. The next two columns, to the right, in statistics were both labeled with Bonferroni_Pvalue
• The equation =S2*5221 was typed into the first cell under the first Bonferroni_Pvalue header; the formula was copied down the entire column
• Any corrected p-value that is greater than 1 was replaced with the number 1 by typing =IF(T2>1,1,T2) into the first cell below the second Bonferroni_Pvalue header; the formula was copied throughout the entire column

Calculating the Benjamini & Hochberg p-value Correction

• A new worksheet named B-H_Pvalue was inserted
• The ID column from the previous worksheet was copied and pasted into the first column of the new worksheet
• A new column was inserted to the very left and labeled as MasterIndex (purpose is to create a numerical index of genes)
  • In the MasterIndex column, a "1" was typed into cell A2 and a "2" was typed into cell A3
  • Both cells were selected and the bottom-right corner (where the cursor becomes a thin black plus sign) was double-clicked. This filled the entire column with the numbers 1 to 5221 (# of genes)
• Using Paste special > Paste values, the unadjusted p-values from the previous worksheet were copied and pasted into column C of this worksheet
• Columns A, B, and C were all selected and sort by ascending values was performed on Column C (sort button on toolbar -> custom sort -> Sort by Pvalue, Sort on Values, Order Smallest to Largest)
• The header Rank was typed into cell D1. "1" was typed into cell D2 and "2" was typed into cell D3; both cells were selected and the double-clicking of the lower right corner was employed in order to fill the column with a series of numbers from 1 to 5221.
• B-H_Pvalue was typed into cell E1 and the formula =(C2*5221)/D2 was typed into cell E2; the equation was copied down the entire column
• B-H_Pvalue was typed into cell F1 and the formula =IF(E2>1,1,E2) was typed into cell F2; the equation was copied down the entire column
• Columns A through F were selected and the columns were sorted by the MasterIndex in column A in ascending order
• Column F was copied and the values were pasted via Paste special into the next column on the right

Preparing file for GenMAPP

• A new worksheet was inserted with the name forGenMAPP
• Everything in the statistics worksheet was selected and copied over to forGenMAPP via Past special > values
• Columns B through Q were selected and the number of decimal places was set to 2 via Format > Cells > number tab, set to 2 decimal places
• All of the columns containing p-values were selected and the number of decimal places was set to 4
• The left-most Bonferroni p-value column was deleted (the one with an "if" statement was kept)
• A column to the right of ID column was inserted and given the header SystemCode. The entire column was filled with the letter "N".
• While on the forGenMAPP worksheet, the file was saved as "Text (Tab-delimited) (*.txt)"
• The resulting file was checked and opened via notepad

Sanity Check: Number of genes significantly changed

• The spreadsheet was opened and the forGenMAPP worksheet was selected
• Cell A1 was clicked and the and the autofilter was turned on via Data > Filter > Autofilter
• The drop-down arrow on the Pvalue column was clicked. "Number filters" was selected, then "Less than...", and then "0.05" was typed into the window that appeared, in order to filter the "Pvalue" column so that only p-values that are less than 0.05 appear
  • 948 genes out of 5221 were found to have a p-value < 0.05, which is 18.16% of genes
• The Pvalue column was then filtered so that only p-values that are less than 0.01 appear
  • 235 genes out of 5221 were found to have a p-value < 0.01, which is 4.50% of genes
• The Pvalue column was then filtered so that only p-values that are less than 0.001 appear
  • 24 genes out of 5221 were found to have a p-value < 0.001, which is 0.46% of genes
• The Pvalue column was then filtered so that only p-values that are less than 0.0001 appear
  • 2 genes out of 5221 were found to have a p-value < 0.0001, which is 0.038% of genes
• Bonferroni_Pvalue was then filtered in order to determine the genes that are p < 0.05 for the Bonferroni-corrected p-value
  • 6 genes out of 5221 were found to have a p < 0.05 for the Bonferroni-corrected p-value, which is 0.115% of genes
• B-H_Pvalue was then filtered in order to determine the genes that are p < 0.05 for the Benjamini and Hochberg-corrected p value
  • 0 genes out of 5221 were found to have a p < 0.05 for the Benjamini and Hochberg-corrected p value, which is 0% of genes
• Avg_LogFC_all was then filtered in order to only show the genes with an average log fold change greater than zero (while keeping the p-value filter at less than 0.05)
  • 352 genes were found to have an average log fold change greater than zero, which is 6.74% of genes
• Avg_LogFC_all was then filtered in order to only show the genes with an average log fold change less than zero (while keeping the p-value filter at less than 0.05)
  • 596 genes were found to have an average log fold change less than zero, which is 11.42% of genes
• Avg_LogFC_all was then filtered in order to only show the genes with an average log fold change > 0.25 (while keeping the p-value filter at less than 0.05)
  • 339 genes were found to have an average log fold change > 0.25, which is 6.49% of genes
• Avg_LogFC_all was then filtered in order to only show the genes with an average log fold change < -0.25 (while keeping the p-value filter at less than 0.05)
  • 579 genes were found to have an average log fold change < -0.25, which is 11.09% of genes
• What criteria did Merrell et al. (2002) use to determine a significant gene expression change? How does it compare to our method?
  • Merrell et al. employed the Statistical Analysis for Microarrays (SAM) program with the intensity ratios in order to "identify significant differences in gene expression"; they conducted a two-class SAM analysis, with the in vitro strain being class I and the individual stool samples being class II. Merrell et al. selected genes with statistically significant changes in expression (which was at least two fold) in each patient sample and this individual stool sample data was used to identify genes that were significantly changed (in expression) in all three samples. The method used by Merrel et al. is somewhat similar to the method used in this investigation but it involved the use of an other computer program (SAM) and the selection of genes that had at least two fold changes in expression. This method used in this investigation primarily involved the use of p-values (that are less than 0.05) in order to identify significant changes in gene expression; similar to what was done by Merrel et al., this method also involved the use of the data from all three patients in order to find changes in gene expression that are significant among all three samples.

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. These genes were looked up in the spreadsheet and their fold changes, p-values, and significance was noted:

• VC0028 (2 entries were found)
  • Fold Change: first entry = 1.65, second entry = 1.27
  • P-Value: first entry = 0.0474, second entry = 0.0692
  • Significance: first entry = statistically significant, second entry = not statistically significant
• VC0941 (2 entries were found)
  • Fold Change: first entry = 0.09, second entry = -0.28
  • P-Value: first entry = 0.6759, second entry = 0.1636
  • Significance: first entry = not statistically significant, second entry = not statistically significant
• VC0869 (5 entries were found)
  • Fold Change (nth entry): 1 = 1.59, 2 = 1.95, 3 = 2.20, 4 = 1.50, 5 = 2.12
  • P-Value (nth entry): 1 = 0.0463, 2 = 0.0227, 3 = 0.0020, 4 = 0.0174, 5 = 0.0200
  • Significance (nth entry): 1 = significant, 2 = significant, 3 = significant, 4 = significant, 5 = significant
• VC0051 (2 entries were found)
  • Fold Change: first entry = 1.92, second entry = 1.89
  • P-Value: first entry = 0.0139, second entry = 0.0160
  • Significance: first entry = statistically significant, second entry = statistically significant
• VC0468
  • Fold Change: -0.17
  • P-Value: 0.3350
  • Significance: not statistically significant
• VC2350
  • Fold Change: -2.40
  • P-Value: 0.0130
  • Significance: statistically significant
• VCA0583
  • Fold Change: 1.06
  • P-Value: 0.1011
  • Significance: not statistically significant

Statistical Analysis of Vibrio cholerae Microarray Data (Part 2)

• GenMAPP was launched and the 2009 Vibrio Cholera database was downloaded and loaded into the program (placed into C:\GenMAPP 2 Data\Gene Databases)
• The data Menu from the main Drafting Board window was selected and then Expression Dataset Manager was chosen. In the Expression Dataset Manager window, New Dataset was selected and then the tab-delimited text file that was formatted for GenMAPP (.txt)
• Expression Dataset Manager was allowed to convert the data and create a new converted dataset

Error Analysis

• After conversion, it was found that 772 errors were detected in the raw data by genMAPP using the 2009 database
• My partner, Anindita V., found that 121 errors were detected in the raw data by genMAPP using the 2010 database
• The .EX.txt file generated by GenMAPP was opened (which contains error messages, along with the raw data) and analyzed
  • It was found that the error message for the unprocessed genes was Gene not found in OrderedLocusNames or any related system.
  • Compared to my partner's results (121 errors with the 2010 database), I had many more errors (772) using the 2009 database. Given that the error message for the 772 errors was Gene not found in OrderedLocusNames or any related system., it appears that the old database covers less genes than the new one (it appears that GenMAPP is giving an error because some of the genes that are present in the data set are not present in the 2009 Database)

Creating Color Sets

• Increased and Decreased LogFoldChange color sets were created in GenMAPP by going to the Expression Dataset Manager and filling in these fields in the Color Set area: name for the Color Set, the gene value, and the criteria that determines how a gene object is colored (on the MAPP)
  • The name of the color set was set as LogFoldChange and Avg_LogFC_all was used as the Gene Value
  • For the Increased criterion (increased LogFoldChange), the name was set as Increased, red was used as the color, and the criterion was [Avg_LogFC_all] > 0.25 AND [Pvalue] < 0.05
  • For the decreased criterion (decreased LogFoldChange), the name was set as Decreased, green was used as the color, and the criterion was [Avg_LogFC_all] < -0.25 AND [Pvalue] < 0.05
  • The whole Expression Dataset was saved and the Expression Dataset Manager was exited

MAPPFinder Procedure

• Assigned criterion: Increased, using 2009 Database
• The MAPPFinder program was launched within GenMAPP (Tools > MAPPFinder)
• "Calculate New Results" was clicked in the window that appeared by launching MAPPFinder
• For "Find File", the Expression Dataset file (with a .gex extension) was selected, and OK was clicked
• The LogFoldChange color set was selected and the Increased criterion was selected (to filter the data)
• The boxes corresponding to "Gene Ontology" and "p value" were checked
• "Browse" button was clicked to add a name to the file that will be created
• "Run MAPPFinder" was clicked and the program was allowed to complete its analysis
• "Show Ranked List" was clicked to see a list of the most significant Gene Ontology terms
• Top 10 Ranked GO Terms

1. biopolymer biosynthetic process
2. macromolecule biosynthetic process
3. macromolecule metabolic process
4. localization
5. transporter activity
6. cellular biopolymer biosynthetic process
7. cellular macromolecule metabolic process
8. transport
9. establishment of localization
10. biopolymer metabolic process

• The top 10 GO terms, found by my partner Anu using the 2010 Database are:
  • 1. branched chain family amino acid metabolic process
    2. branched chain family amino acid biosynthetic process
    3. IMP metabolic process
    4. IMP biosynthetic process
    5. purine ribonucleoside monophosphate metabolic process
    6. purine nucleoside monophosphate metabolic process
    7. purine nucleoside monophosphate biosynthetic process
    8. purine ribonucleoside monophosphate biosynthetic process
    9. arginine metabolic process
    10. cellular nitrogen compound biosynthetic process
  • The top 10 GO terms, compared between the 2009 and 2010 database, appear to mostly involve metabolic/biosynthetic processes. The 2010 database seems to have top GO terms that are more specific (than the 2009 database) and most of these terms involve purine metabolism and metabolic pathways tied to amino acids (and other compounds). Unlike the 2010 database, the 2009 database had top GO terms that involved transport (transporter activity) and movement (localization). It is suspected that the variation for the top GO terms is due to differences between the 2009 and 2010 databases with respect to the number of genes covered/represented. The 2009 database had many more errors, and thus, many more genes that were not taken into account by GenMAPP; it is possible that the genes represented by the GO terms related to the amino acid/purine metabolic pathways are not present in the GenMAPP analysis with the 2009 database (potentially, these genes are the ones that are tied to errors).
• In the main MAPPFinder Browser window, "Collapse the Tree" was clicked on and these genes were searched for (one by one): VC0028, VC0941, VC0869, VC0051, VC0647, VC0468, VC2350, and VCA0583. The ID for the genes was put in the gene ID search field and "OrderedLocusNames" was selected from the drop-down menu to the right of the search field. GeneID search button was clicked in order to commence the search for the gene. These are the genes mentioned by Merrell et al. in their paper.
• VC0028: "No MAPPs or GO terms could be found for that OrderedLocusNames ID."
• VC0941: "No MAPPs or GO terms could be found for that OrderedLocusNames ID."
• VC0869: "No MAPPs or GO terms could be found for that OrderedLocusNames ID."
• VC0051: "No MAPPs or GO terms could be found for that OrderedLocusNames ID."
• VC0647: 3'-5'-exoribonuclease activity, transferase activity, nucleotidyltransferase activity, polyribonucleotide nucleotidyltransferase activity
• VC0468: "No MAPPs or GO terms could be found for that OrderedLocusNames ID."
• VC2350: "No MAPPs or GO terms could be found for that OrderedLocusNames ID."
• VCA0583: transport, outer membrane-bounded periplasmic space, transporter activity
  • My partner, Anu, had very different results for these genes. The results were exactly the same for VC0647 and very similar for VCA0583; it is very likely that these genes are both covered in the same way or a similar way in the two databases. Anu found GO terms associated with the rest of the genes (using the 2010 database), however, using the 2009 database, MAPPfinder returned the prompt "No MAPPs or GO terms could be found for that OrderedLocusNames ID.". The prompt for most of the genes, along with the high error count, suggests that the 2009 database did not include these genes (that were mentioned by Merrell et al.). Judging by the search results, it appears that the 2010 database includes all of these genes while the 2009 database only includes a few.

GO Term Investigation for VCA0583, with UniProt

• outer membrane-bounded periplasmic space, corresponding to VCA0583, was clicked and a MAPP listing all of the genes associated with that GO term opened up.
• In the window that appeared, Q9KP40_VIBCH was selected
• Q9KP40_VIBCH had a significant increase of expression (1.67 LogFoldChange)
• Using the UniProt database, details regarding Q9KP40_VIBCH were found
  • According to UniProt, this gene codes for Thiamin ABC transporter/periplasmic thiamin-binding protein
  • Periplasm was searched up; it was found to be a gel-like matrix existing between the inner cytoplasmic membrane and the outer membrane of a bacterium
  • ABC transporter was searched up; it was found to be a class of transmembrane proteins that use ATP (via hydrolysis of the phosphoanhydride bonds) in order to carry out certain processes/functions (such as the movement of substrates across membranes)
  • Thiamin was found to be a vitamin (vitamin B1), that is utilized by many organisms, including gram-negative bacteria like Vibrio cholerae
  • Given this UniProt result, it appears that the gene codes for a protein that binds to thiamin and employs ATP to move it across membranes (is a transmembrane protein that works via active-transport)

Working with the Results File

• The results file, in the format "XXX-CriterionX-GO.txt" was copied and pasted in the folder housing all of the files related to this investigation (XXX = name given to file, CriterionX refers to the Criterion #; Criterion0 refers to the Increased criterion, Criterion1 refers to the Decreased criterion, if the Increased criterion was done first)
• The "CriterionGO" file was opened via Excel in order to show the results in a tabular format
• The top rows of the opened file (in Excel) were examined; these are rows of information that give background information as to how MAPPFinder made the calculations
  • 
• The list in the file was filtered in order to show the top GO terms represented in the data for the "Increased" criterion
• A cell was clicked in the row of headers and the autofilter was turned on (Data menu -> Filter > Autofilter)
• The two custom filters Z Score (in column N) greater than 2 and PermuteP (in column O) less than 0.05 were set for the columns
• The custom filters Number Changed (in column I) greater than or equal to 5 AND less than 100 and Percent Changed (in column L) greater than or equal to 26.5 was also applied to filter down the results to 20 entries
• The top 20 GO terms in the file were searched in MAPPFinder in order to find if any relationships existed between the terms. Terms that had a relationship (direct parent/child) with another term on the list were highlighted; the term that had a relationship with the highlighted term was highlighted with the same color
• Unfamiliar GO terms included:

• Modified Excel File/Data

• For GenMAPP

Working with: 2009 Vibrio Cholera database

• 772 errors were detected in the raw data by genMAPP using the 2009 database
• 121 errors were detected in the raw data by genMAPP using the 2010 database by my partner, Anindita V.

Results of 10/20 Work Session

.gex file for GenMAPP

Results of 10/22 Work Session

• Top 10 Ranked GO Terms

1. macromolecule metabolic process
2. localization
3. transporter activity
4. cellular macromolecule metabolic process
5. transport
6. establishment of localization
7. cell projection organization
8. cellular biopolymer metabolic process
9. macromolecule biosynthetic process
10. biopolymer metabolic process

• Analysis of the .EX.txt file produced by GenMAPP, via Excel, revealed that the 772 errors were: Gene not found in OrderedLocusNames or any related system; this suggests that the 2009 database did not

increased expression

Brandon Litvak
BIOL 367, Fall 2015