Bklein7 Week 8

From LMU BioDB 2015
Revision as of 17:26, 26 October 2015 by Bklein7 (Talk | contribs) (edited link format)

Jump to: navigation, search

Files Generated in the This Week's Analysis

Links to files below can all be found within the electronic notebook at the point where they were created. For easy access, they are listed here as well.

  1. Analyzed microarray data: File:Merrell Compiled Raw Data Vibrio BK 20151015.xls.
  2. For GenMAPP text file: File:Merrell Compiled Raw Data Vibrio BK 20151015- Tab Delimited.txt.
  3. Exceptions file: File:Merrell Compiled Raw Data Vibrio BK 20151015.EX.txt.
  4. Expression Dataset files:
  5. GO term MAPP: File:Polyribonucleotide nucleotidyltransferase activity.mapp.
  6. MAPPFinder Results:
  7. .gmf file: File:Merrell Compiled Raw Data Vibrio BK 20151025.gmf.

Statistical Analysis of Vibrio cholerae Microarray Data (Part 1)

Normalizing the Log Ratios for the Set of Slides in the Experiment

This section dictates the steps necessary to scale and center the raw microarray data:

  • To begin, I created a new Worksheet in my Excel file entitled "scaled_centered".
  • I went back to the original "compiled_raw_data" worksheet and copied over all the data into the new "scaled_centered" worksheet.
  • I inserted two rows in between the top row of headers and the first data row entitled "Average" (cell A2) "StdDev" (cell A3).
    • I computed the Average log ratio for each chip by inputting the following equation into cell B2 and then pasting it into the rest of the "Average" column: =AVERAGE(B4:B5224).
    • I computed the Standard Deviation of the log ratios on each chip by inputting the following equation into cell B3 and then pasting it into the rest of the "StdDev" column: =STDEV(B4:B5224).
  • I created a new set of headings for the scaled and centered data by copying over the data column headings and then pasting them to the right of the last data column. I edited the names of the columns so that they now read: A1_scaled_centered, A2_scaled_centered, etc.
  • In cell N4 (column with the heading A1_scaled_centered), 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?
  • I copied this equation to the rest of the column and then adapted it for all "_scaled_centered" columns.

Performing Statistical Analysis on the Ratios

This section details the steps necessary to perform statistical analysis on the scaled and centered data produced in the section above:

  • For this section, I created a new worksheet and name it "statistics".
  • I went back to the "scaled_centered" worksheet and copied over both the first column ("ID") and copied the columns that are designated "_scaled_centered".
  • I deleted rows 2 and 3 where it says "Average" and "StDev" so that the data rows with gene IDs were immediately below the header row 1.
  • I created three new columns to the right of the copied ones with the following headers: "Avg_LogFC_A", "Avg_LogFC_B", and "Avg_LogFC_C".
  • I computed the average log fold change for the replicates for each patient by typing the equation:
=AVERAGE(B2:E2)
into cell N2 and then copying/pasting it into the rest of the three columns, adapting it as necessary.
  • I created a new column to compute the average of the averages with the header "Avg_LogFC_all". I then created the equation to compute the average of the three previous averages you calculated and paste it into this entire column.
  • I created a new column with the header "Tstat" to compute the T statistic for each scaled and centered average log ratio. In this column, I entered the following equation and pasted it into the rest of the column: 
=AVERAGE(N2:P2)/(STDEV(N2:P2)/SQRT(3))
  • I created one last column with the heading "Pvalue". In the cell below the label, I entered the following equation and then copied it into the rest of the column: 
=TDIST(ABS(R2),2,2)

Calculating the Bonferroni p value Correction

It is necessary to perform adjustments to the p value to correct for the multiple testing problem. Therefore, I first calculated the Bonferroni p value correction on the p values calculated in the above section:

  • I labelled the next two columns to the right of "Pvalue" with the same label: "Bonferroni_Pvalue".
  • I typed the equation =S2*5221 into the first Bonferroni column and pasted it throughout.
  • Next, I replaced any corrected p value that was greater than 1 with the number 1 by typing the following formula into the first cell below the second Bonferroni_Pvalue header: =IF(T2>1,1,T2).

Calculating the Benjamini & Hochberg p value Correction

The second p value correction I performed was the Benjamini & Hochberg correction, the methods of which are presented below:

  • I created a new worksheet named "B-H_Pvalue".
  • In this worksheet, I copied over the "ID" column from the previous worksheet into the first column of the new worksheet.
  • Next, I added a new column on the very left and named it "MasterIndex". Under this heading, I added a numbered list from 1 to 5221 (the number of genes on the microarray).
  • I copied over the unadjusted p values from your previous worksheet and pasted them into Column C.
  • I selected all of columns A, B, and C and sorted by ascending values on Column C.
    • This was done by clicking the sort button from A to Z on the toolbar and sorting by column C, smallest to largest.
  • I created a new column and labelled it with the header "Rank" in cell D1 to house another numbered list from 1 to 5221. Because this was done for the sorted columns, these values represented the p value ranks, smallest to largest.
  • I created a two new columns to calculate the Benjamini and Hochberg p value correction and labeled them "B-H_Pvalue" (cell E1 and F1).
    • In the first column, I inputted the following formula and copied it throughout the column: =(C2*5221)/D2.
    • In the second column, I inputted the following formula and copied it throughout: =IF(E2>1,1,E2).
  • Next, I selected columns A through F and sorted them by the MasterIndex in Column A in ascending order.
  • Finally, I copied column F into the next column on the right of your "statistics" sheet.

Preparing the File for GenMAPP

  • I inserted a new worksheet and named it "forGenMAPP".
  • I copied over the entirety of the "statistics" worksheet to this new worksheet.
  • I selected Columns B through Q (all the fold changes) and formatted the cells to show only 2 decimal places.
  • I selected all the columns containing p values and formatted the cells to show only 4 decimal places.
  • I deleted the left-most Bonferroni p value column.
  • I inserted a column to the right of the "ID" column entitled "SystemCode" and filled the entire column with the letter "N".
  • After making the above changes, I saved the file as "Text (Tab-delimited) (*.txt)".

Sanity Check: Number of Genes Significantly Changed

To verify the results of the data analysis performed on the Vibrio cholerae DNA microarray data, I assessed the number of genes that were significantly changed at various p value cut-offs and also compared my data analysis with the published results of Merrell et al. (2002). These analyses were done in the "forGenMapp" tab of the Excel spreadsheet used for the data analysis.

  • Assessing the number of genes significantly changed
    • To answer the questions below, I used custom filter to display only the "Pvalue" data that met specific criterion (e.g. < 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)?
  • ((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. 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)?
  • Assessing the magnitude and direction of significant gene expression changes
    • The "Avg_LogFC_all" indicated the size of the gene expression changes and their direction. Positive values increased relative to the control, and negative values decreased relative to the control.
      • Keeping the (unadjusted) "Pvalue" filter at p < 0.05, I filtered 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, I filtered 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 %)
  • What criteria did Merrell et al. (2002) use to determine a significant gene expression change? How does it compare to our method?
  • For the GenMAPP analysis below, I used 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 the analysis.

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. I reviewed the "forGenMAPP" Excel worksheet to compare these findings to the results of my personal data analysis.

  • What are the fold changes and p values (of these genes)? Are they significantly changed in the analysis?
    • VC0028
    • VC0941
    • VC0869
    • VC0051
    • VC0647
    • VC0468
    • VC2350
    • VCA0583

MAPPFinder Analysis of Vibrio cholerae Microarray Data (Part 2)

Mapping Onto Biological Pathways (GenMAPP & MAPPFinder)

Before beginning the mapping process, it is necessary to load the correct gene database into GenMAPP:

  • I download the 2009 Vibrio cholerae gene database by following this link to the XMLPipeDB SourceForge Download page. My homework partner Veronica downloaded the 2010 version.
  • I saved the above file into the folder C:\GenMAPP 2 Data\Gene Databases and extracted it.
  • Within GenMAPP, I loaded the 2009 gene database by selecting Data > Choose Gene Database and choosing the appropriate file from the directory C:\GenMAPP 2 Data\Gene Databases.

GenMAPP Expression Dataset Manager Procedure

Once the appropriate gene database has been loaded into GenMAPP, the expression dataset can be uploaded and configured:

  • I opened the Expression Dataset Manger from the Data drop-down list in GenMAPP.
  • I selected New Dataset from the Expression Datasets menu and choose the tab-delimited text file formatted for GenMAPP (.txt).
  • Upon specifying that all data was numerical, the Expression Dataset Manager converted my data to .gex file. This process took approximately one minute to complete. In addition to converting the data to a .gex file, an exceptions file (.EX.txt) was also produced, as 772 errors were reportedly detected in the raw data.
    • The .gex file generated can be found here: File:Merrell Compiled Raw Data Vibrio BK 20151015.gex
    • The .EX file generated can be found here: File:Merrell Compiled Raw Data Vibrio BK 20151015.EX.txt
    • 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.
      • GenMAPP errors message.PNG
      • The above screenshot shows the error message I received after using the Expression Dataset Manager to convert my raw data file. A reported 772 errors were detected in my raw data.
      • Upon opening the .EX.txt file, I found the error message "Gene not found in OrderedLocusNames or any related system." repeated several times.
        • To determine if this was the only error message, I uploaded the exceptions file to my LMU CMSI directory using the command scp /Users/brandonklein/Desktop/Merrell_Compiled_Raw_Data_Vibrio_BK_20121015.EX.txt bklein7@my.cs.lmu.edu:~bklein7. In the command line, I used the following command sequence to determine how often the "Gene not found..." error message was repeated: grep "Gene not found"Merrell_Compiled_Raw_Data_Vibrio_BK_20121015.EX.txt | wc. This yielded the output 772 23932 145791, confirming that the above error message was responsible for all 772 reported errors.
      • 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?
        • With the 2009 Vibrio cholerae gene database loaded into GenMAPP, I encountered 772 errors when converting the analyzed Merrell et al. microarray data. When converting the exact same data with the 2010 Vibrio cholera gene database loaded into GenMapp, Veronica encountered 121 errors. All error messages for both of us were the same: "Gene not found in OrderedLocusNames or any related system". Therefore, I had more errors. I believe this occured because the 2009 version of the Vibrio cholerae database was less developed/complete than the more recent 2010 version. Thus, it likely had less Gene listings, and therefore less of the Ordered Locus Names used by Merrell matched to the database.
  • I customized the new Expression Dataset by creating a Color Sets= with instructions to GenMAPP for displaying data on MAPPs. The new Color Set was entitled "LogFoldChange".
    • First, I created a criterion for this color set to label genes that demonstrated a significant increase in their expression.
      • I specified the Gene value as "Avg_LogFC_all" for the Vibrio dataset.
      • I activated the Criteria Builder by clicking the New button and named the criterion "Increased".
      • I selected the color for this criterion as red using the color box.
      • I stated the criterion as follows and added it to the Criteria List: [Avg_LogFC_all] > 0.25 AND [Pvalue] < 0.05
    • Second, I created a criterion for this color set to label genes that demonstrated a significant decrease in their expression.
      • I specified the Gene value as "Avg_LogFC_all" for the Vibrio dataset.
      • I activated the Criteria Builder by clicking the New button and named the criterion "Decreased".
      • I selected the color for this criterion as green using the color box.
      • I stated the criterion as follows and added it to the Criteria List: [Avg_LogFC_all] < -0.25 AND [Pvalue] < 0.05
  • Upon entering these color sets, I savedthe entire Expression Dataset by selecting Save from the Expression Dataset menu.

MAPPFinder Procedure

I used MAPPFinder to Note: Veronica and I both used the "Decreased" criterion.

  • I launched the MAPPFinder program from within GenMAPP and ensured that the 2009 Gene Database was still loaded into GenMAPP.
  • I clicked on the button "Calculate New Results" followed by "Find File", at which point I chose my .gex file.
  • Veronica and I both chose to filter the data with the "Decreased" criterion present within the LogFoldChange Color Set.
  • I checked the boxes next to "Gene Ontology" and "p value", specified the results file, and then clicked "Run MAPPFinder".
    • This analysis took several minutes to complete.
  • I clicked on the menu item "Show Ranked List" to see a list of the most significant Gene Ontology terms.
    • List the top 10 Gene Ontology terms.
      • TOP10.PNG
      1. protein folding
      2. chorismate metabolic process
      3. aromatic amino acid family biosynthetic process
      4. unfolded protein binding
      5. cytoplasm
      6. intracellular part
      7. localization
      8. nucleotide catabolic process
      9. locomotion
      10. aromatic amino acid family 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?
      • Upon comparing my results with Veronica (who used the 2010 Vibrio database), I found that our list of top 10 Gene Ontology terms were entirely different. ... (different genes in database + different maps?)
  • I used MAPPFinder to find the Gene Ontology term(s) with which the following genes mentioned by Merrell et al. (2002) were associated with: VC0028, VC0941, VC0869, VC0051, VC0647, VC0468, VC2350, and VCA0583. This was done by typing the identifiers for these genes into the MAPPFinder browser gene ID search field, choosing "OrderedLocusNames" from the drop-down menu to the right of the search field, and clicking on the GeneID Search button. The GO term(s) that were associated with the genes were highlighted in blue. List the GO terms associated with each of those genes. (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?
    • 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: mRNA catabolic process, RNA processing, cytoplasm, RNA binding, 3'-5'-exoribonuclease activity, transferase activity, nucleotidyltransferase activity, and 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, and transporter activity.
  • VC0647 Investigation: List 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.
    • I clicked on the GO term "polyribonucleotide nucleotidyltransferase activity". The MAPP this produced only contained the gene "PNP_VIBCH", which I found in the [UniProt Database]. The expression of this gene, referenced by the gene ID VC0647, significantly decreased in the microarray experiment.
    • I Double-clicked on the above gene box to find out more about the gene PNP_VIBCH. Click on the links to find out the function of this gene.
      • According to the [UniProt Database], the function of this gene is to degrade mRNA. It does this by coding for the production of the protein polyribonucleotide nucleotidyltransferase, which "catalyzes the phosphorolysis of single-stranded polyribonucleotides processively in the 3'- to 5'-direction".
    • The MAPP created while investigating polyribonucleotide nucleotidyltransferase activity can be found here: File:Polyribonucleotide nucleotidyltransferase activity.mapp.
  • Overall MAPPFinder results file: File:Merrell Compiled Raw Data Vibrio BK 20151015-Criterion1-GO.txt.
  • I opened a copy of the .txt file listed above in Excel to filter the results.
    • There were rows of information that gave background information on how MAPPFinder made the calculations at the top of the Excel file. Compare this information with your partner who used a different version of the Vibrio Gene Database. Which numbers are different? Why are they different?
      • ... TO DO ...
    • I used the following filters to show the top 20 GO terms represented in my data for both the "Increased" and "Decreased" criteria:
Z Score (in column N) greater than 2
PermuteP (in column O) less than 0.05
Number Changed (in column I) greater than or equal to 5 AND less than 100
Percent Changed (in column L) greater than or equal to 30%
  • 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.
  • 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.
    • TO DO
    • Top 20 GO Terms
      • group 1: protein folding / unfolded protein binding
      • group 1: cis-trans isomerase activity / peptidyl-prolyl cis-trans isomerase activity
        • decreased (~1.5) -> promotes proper protein folding (helps chaperone) & prevents protein aggregation
      • zinc ion binding
        • decreased -> protein involved in glycolytic process
        • increased pepsidase activity again
      • group 2: sugar:hydrogen symporter activity / solute:hydrogen symporter activity / cation:sugar symporter activity
        • decreased sugar-dependent active transport into membrane
      • group 3: protein-N(PI)-phosphohistidine-sugar phosphotransferase activity / sugar transmembrane transporter activity / carbohydrate transmembrane transporter activity
        • same as above
      • group 4: phosphoenolpyruvate-dependent sugar phosphotransferase system / carbohydrate transport
        • (decreased) same as above... + intake of sugars from external environment decreased
      • group 5: translation regulator activity / translation factor activity, nucleic acid binding
        • decreased stimulation of aa-tRNA binding -> decreased protein production
      • endonuclease activity
        • decreased hydrolysis of nucleic acids
      • glucose metabolic process
      • decreased glycolytic processes
      • aromatic compound biosynthetic process
        • decreased chorismate biosynthesis (pathway)
      • group 6: intracellular protein transport / intracellular transport
        • decreased intracellular protein transport
  • This is a response to stress present in the human stomach (low pH)? instead of refolding denatured proteins (low pH), they spend less energy on refolding proteins and instead hydrolyze the denatured proteins. In addition, sugars are conserved to survive in a nutrient-poor environment (less glycolysis, sugar-related active transport, and chorismate production) and there is less transport from the external environment. Mutants that express these genes at their altered levels can become opportunistic human pathogens.
  • 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.

Links

Assignments Pages

Individual Journal Entries

Shared Journal Entries

Retrieved from "https://xmlpipedb.lmucs.io/biodb/fall2015/index.php?title=Bklein7_Week_8&oldid=3287"