Difference between revisions of "Kmill104 Week 9"
(→Calculate the Bonferroni and p value Correction: fixing tenses for this section) |
(→Calculate the Benjamini & Hochberg p value Correction: fixing tenses for this section) |
||
Line 53: | Line 53: | ||
==== Calculate the Benjamini & Hochberg p value Correction ==== | ==== Calculate the Benjamini & Hochberg p value Correction ==== | ||
− | # | + | # I then inserted a new worksheet named "wt_ANOVA_B-H". |
− | # | + | # I copy and pasted the "MasterIndex", "ID", and "Standard Name" columns from my previous worksheet into the first two columns of the new worksheet. |
− | # | + | # I used Paste special > Paste values for the following: I copied my unadjusted p values from my ANOVA worksheet and pasted it into Column D. |
− | # | + | # I then selected all of columns A, B, C, and D. I sorted by ascending values on column D by selecting the sort button from A to Z on the toolbar, and in the window that appeared, clicked sort by column D, smallest to largest. |
− | # | + | # I then typed the header "Rank" in cell E1. This column was used to create a series of numbers in ascending order from 1 to 6189. This gave me the p value rank, smallest to largest. To assign ranks, I typed "1" into cell E2 and "2" into cell E3. I selected both cells E2 and E3, and then double-clicked 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 6189. |
− | # | + | # I used the following to calculate the Benjamini and Hochberg p value corrections. I typed wt_B-H_p-value in cell F1, and typed the following formula in cell F2: <code>=(D2*6189)/E2</code> and pressed enter. I then copied that equation to the entire column. |
− | # | + | # I then typed "wt_B-H_p-value" into cell G1. |
− | # | + | # In G2, I typed the following formula: <code>=IF(F2>1,1,F2)</code> and pressed enter, then copied that equation to the entire column. |
− | # | + | # I then selected columns A through G, and sorted them by my MasterIndex in Column A in ascending order. |
− | # | + | # Finally, I copied column G and use Paste special > Paste values to paste it into the next column on the right of my ANOVA sheet. |
− | |||
− | |||
==== Sanity Check: Number of genes significantly changed ==== | ==== Sanity Check: Number of genes significantly changed ==== |
Revision as of 21:25, 1 April 2024
Purpose
The purpose of this week's journal assignment is to analyze DNA microarray data using several different equations in an Excel format. It is also to learn about the meaning of and importance of p-values, and the ways that we can adjust p-values through the use of additional equations. It is also to keep an organized and detailed electronic notebook that makes our research able to be replicated and reproduced.
Methods/Results
Experimental Design and Getting Ready
The data used in this exercise is publicly available at the NCBI GEO database in record GSE83656.
- I first began by downloading the Excel file for my group's strain.
- I then changed the filename so that it contained my initials to distinguish it from other students' work.
- I then looked at the sheet labeled "Master_Sheet_wt" to find that:
- The strain that me and Dean will analyze is the wild type strain. The name of the original file is BIOL367_S24_microarray-data_wt.xlsx. The name of the file with my saved initials and what I will be working on is BIOL367_S24_microarray-data_wt_KM.xlsx. For the total wt strain data, there are 23 replicates. At time point 15, there are 4 replicates. At time point 30, there are 5 replicates. At time point 60, there are 4 replicates. At time point 90, there are 5 replicates. And at time point 120, there are 5 replicates.
Statistical Analysis Part 1: ANOVA
- I then created a new worksheet, naming it "wt_ANOVA".
- I then copied all data from the "Master_Sheet" worksheet and pasted it into my new worksheet.
- At the top of the first column to the right of your data, I then created five column headers of the form wt_AvgLogFC_(TIME) where (TIME) is 15, 30, 60, 90, and 120.
- In the cell below the wt_AvgLogFC_t15 header, I then typed
=AVERAGE(
- I then highlighted all the data in row 2 associated with t15 (D1:G1), pressed the closing paren key (shift 0), and pressed the "enter" key.
- The cell now contains the average of the log fold change data from the first gene at t=15 minutes.
- I then clicked on this cell and positioned my cursor at the bottom right corner. After the cursor changed to a thin black plus sign, I double clicked and the formula was copied to the remaining genes.
- I then repeated steps (4) through (8) with the t30 (H1:L1), t60 (M1:P1), t90 (Q1:U1), and the t120 data (V1:Z1).
- To the right of the wt_AvgLogFC_t120 calculation, I created the column header wt_ss_HO.
- In the first cell below this header, I typed
=SUMSQ(
- Then, I highlighted all the LogFC data in row 2 (D1:Z1), pressed the closing paren key (shift 0),and pressed the "enter" key.
- In the next empty column to the right of wt_ss_HO, I created the column headers wt_ss_(TIME) as in (3).
- I counted that there are 23 data points for each time point of the wt strain data.
- In the first cell below the header wt_ss_t15, I then typed
=SUMSQ(D2:G2)-COUNTA(D2:G2)*<AA2>^2
and hit enter.- Upon completion of this single computation, I then used the Step (7) trick to copy the formula throughout the column.
- I then repeated this computation for the t30 through t120 data points.
- Below wt_ss_t30, I typed
=SUMSQ(H2:L2)-COUNTA(H2:L2)*<AB2>^2
and hit enter, then copied the formula throughout the column. - Below wt_ss_t60, I typed
=SUMSQ(M2:P2)-COUNTA(M2:P2)*<AC2>^2
and hit enter, then copied the formula throughout the column. - Below wt_ss_t90, I typed
=SUMSQ(Q2:U2)-COUNTA(Q2:U2)*<AD>^2
and hit enter, then copied the formula throughout the column. - Below wt_ss_t120, I typed
=SUMSQ(V2:Z2)-COUNTA(V2:Z2)*<AE2>^2
and hit enter, then copied the formula throughout the column.
- Below wt_ss_t30, I typed
- In the first column to the right of wt_ss_t120, I then created the column header wt_ss_full.
- In the first row below this header, I typed
=sum(AG2:AK2)
and hit enter. - In the next two columns to the right, I then created create the headers wt_Fstat and wt_p-value.
- There are 23 data points for each time point, which I used in the following formula.
- In the first cell of the wt_Fstat column, I typed
=(23-5)/5)*(AF2-AL2)/AL2
and hit enter.- I then copied this to the whole column.
- In the first cell below the wt_p-value header, I typed
=FDIST(AM2,5,23-5)
, and then copied this to the whole column. - Before moving onto the next step, I performed a quick sanity check to see if the computations were done correctly.
- First, I clicked on cell A1 and clicked on the Data tab. I then selected the Filter icon, and used the drop-down arrows that appeared to filter the data.
- I clicked on the drop-down arrow on my wt_p-value column, then selected "Number Filters". In the window that appeared, I set a criterion that filtered my data so that the p value had to be less than 0.05.
- Excel then only displayed the rows that corresponded to data meeting that filtering criterion. A number appeared in the lower left hand corner of the window giving me the number of rows that meet that criterion. These numbers were then used to check our results with each other to make sure that the computations were performed correctly.
Calculate the Bonferroni and p value Correction
- Then, I performed adjustments to the p value to correct for the multiple testing problem. I labeled the next two columns to the right of wt_p-value with the same label, wt_Bonferroni_p-value.
- In the first cell of the first column, I typed the equation
=AN2*6189
, and upon completion of this single computation, copied the formula throughout the entire column. - I then 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 wt_Bonferroni_p-value column:
=IF(AO2>1,1,AO2)
, and then copied the formula throughout the column.
Calculate the Benjamini & Hochberg p value Correction
- I then inserted a new worksheet named "wt_ANOVA_B-H".
- I copy and pasted the "MasterIndex", "ID", and "Standard Name" columns from my previous worksheet into the first two columns of the new worksheet.
- I used Paste special > Paste values for the following: I copied my unadjusted p values from my ANOVA worksheet and pasted it into Column D.
- I then selected all of columns A, B, C, and D. I sorted by ascending values on column D by selecting the sort button from A to Z on the toolbar, and in the window that appeared, clicked sort by column D, smallest to largest.
- I then typed the header "Rank" in cell E1. This column was used to create a series of numbers in ascending order from 1 to 6189. This gave me the p value rank, smallest to largest. To assign ranks, I typed "1" into cell E2 and "2" into cell E3. I selected both cells E2 and E3, and then double-clicked 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 6189.
- I used the following to calculate the Benjamini and Hochberg p value corrections. I typed wt_B-H_p-value in cell F1, and typed the following formula in cell F2:
=(D2*6189)/E2
and pressed enter. I then copied that equation to the entire column. - I then typed "wt_B-H_p-value" into cell G1.
- In G2, I typed the following formula:
=IF(F2>1,1,F2)
and pressed enter, then copied that equation to the entire column. - I then selected columns A through G, and sorted them by my MasterIndex in Column A in ascending order.
- Finally, I copied column G and use Paste special > Paste values to paste it into the next column on the right of my ANOVA sheet.
Sanity Check: Number of genes significantly changed
Before we move on to further analysis of the data, we want to perform a more extensive 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.
- Go to your wt_ANOVA worksheet.
- Select row 1 (the row with your column headers) and select the menu item Data > Filter > Autofilter (The funnel icon on the Data tab). 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 for the unadjusted p value. 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 < 0.05? and what is the percentage (out of 6189)? 2528 genes have p < 0.05. This is 40.85% of the total genes in the dataset.
- How many genes have p < 0.01? and what is the percentage (out of 6189)? 1652 genes have p < 0.01. This is 26.69% of the total genes in the dataset.
- How many genes have p < 0.001? and what is the percentage (out of 6189)? 919 genes have p < 0.001. This is 14.85% of the total genes in the dataset.
- How many genes have p < 0.0001? and what is the percentage (out of 6189)? 496 genes have p < 0.0001. This is 8.01% of the total genes in the dataset.
- Note that it is a good idea to create a new worksheet in your workbook to record the answers to these questions. Then you can write a formula in Excel to automatically calculate the percentage for you.
- 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 by chance less than 5% of the time.
- We have just performed 6189 hypothesis tests. Another way to state what we are seeing with p < 0.05 is that we would expect to see this a gene expression change for at least one of the timepoints by chance in about 5% of our tests, or 309 times. Since we have more than 309 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 6189)? 248 genes are p < 0.05. This is 4.01% of the total genes in the dataset.
- How many genes are p < 0.05 for the Benjamini and Hochberg-corrected p value? and what is the percentage (out of 6189)? 1822 genes are p < 0.05. This is 29.44% of the total genes in the dataset.
- 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.
- We will compare the numbers we get between the wild type strain and the other strains studied, organized as a table. Use this sample PowerPoint slide to see how your table should be formatted. Upload your slide to the wiki.
- Note that since the wild type data is being analyzed by one of the groups in the class, it will be sufficient for this week to supply just the data for your strain. We will do the comparison with wild type at a later date.
- Comparing results with known data: the expression of the gene NSR1 (ID: YGR159C)is known to be induced by cold shock. Find NSR1 in your dataset. What is its unadjusted, Bonferroni-corrected, and B-H-corrected p values?
- Unadjusted: 2.86939E-10
- Bonferroni-corrected: 1.77586E-06
- B-H-corrected: 8.87932E-07
- What is its average Log fold change at each of the timepoints in the experiment?
- Time point 15: 3.279225
- Time point 30: 3.621
- Time point 60: 3.526525
- Time point 90: -2.04985
- Time point 120: -0.60622
- Note that the average Log fold change is what we called "wt_AvgLogFC_(TIME)" in step 3 of the ANOVA analysis. Does NSR1 change expression due to cold shock in this experiment?
- Yes, the unadjusted, Bonferroni-corrected, and B-H-corrected p-values are all less than 0.05, indicating that NSR1 changes expression due to cold shock in this experiment.
- For fun, find MSN1 (from your Week 3 assignment) in the dataset. What is its unadjusted, Bonferroni-corrected, and B-H-corrected p values?
- Unadjusted: 0.563798852
- Bonferroni-corrected: 3489.351093
- B-H-corrected: 0.679258535
- What is its average Log fold change at each of the timepoints in the experiment?
- Time point 15: 0.1076
- Time point 30: -0.46192
- Time point 60: -0.47075
- Time point 90: 0.16805
- Time point 120: -0.18418
- Does your favorite gene change expression due to cold shock in this experiment?
- No, the unadjusted, Bonferroni-corrected, and B-H-corrected values are all greater than 0.05, indicating the gene expression has not changed due to cold shock in the experiment.
Data & Files
Conclusion
During this week's journal assignment, I learned ways to analyze datasets with different equations, and how to use Excel for getting correct and efficient results. I discovered that the meaning of p-values is that , and we can use the Bonferroni and Benjamini and Hochberg p-value correction equations for our obtained p-value data. I also learned that as p-values get smaller, less genes fall in the parameter of being less than that value. I saw that log-fold change is different for the NSR1 and MSN1 genes at each timepoint. Finally, during this week I followed the Week 9 assignment page to keep my own organized and detailed electronic journal that is specific to the wild type strain.
Acknowledgements
I worked under the guidance of Dr. Dahlquist on 3-14-24 and 3-18-24. I consulted with my classmates during class time when I had a question about 21. in the ANOVA section of this journal assignment. I first copied the general procedure from the Week 9 Assignment page, and then adjusted my individual procedure to what exactly I did when following the steps.
Except for what is noted above, this individual journal entry was completed by me and not copied from another source.
Kmill104 (talk) 20:30, 20 March 2024 (PDT)
References
LMU BioDB 2024. (2024). Week 9. Retrieved March 20, 2024, from https://xmlpipedb.cs.lmu.edu/biodb/spring2024/index.php/Week_9
User Page
Assignment Pages
- Week 1
- Week 2
- Week 3
- Week 4
- Week 5
- Week 6
- Week 8
- Week 9
- Week 10
- Week 11
- Week 12
- Week 13
- Week 14
- Week 15
Individual Journal Entry Pages
- Kmill104 Week 1
- Kmill104 Week 2
- MSymond1 KMill104 Week 3
- Monarch Initiative Week 4
- Kmill104 Week 5
- Kmill104 Week 6
- Kmill104 Week 8
- Kmill104 Week 9
- Kmill104 Week 10
- Kmill104 Week 11
- Kmill104 Week 12
- Data Analysts Week 13
- Data Analysts Week 14
- Data Analysts Week 15