Noi/NOTES/2012-4-11: Difference between revisions
Jump to navigation
Jump to search
>Noi No edit summary |
>Noi No edit summary |
||
(9 intermediate revisions by the same user not shown) | |||
Line 39: | Line 39: | ||
Note: Both commands gave the same results of q-value or FDR. Length could be set longer than the length of input p-value. | Note: Both commands gave the same results of q-value or FDR. Length could be set longer than the length of input p-value. | ||
* Links: http://rss.acs.unt.edu/Rdoc/library/stats/html/p.adjust.html, http://www.cbs.umn.edu/courses/katagiri/lecture_notes/081103Lec_Statistics3.pdf | * Links: http://rss.acs.unt.edu/Rdoc/library/stats/html/p.adjust.html, http://www.cbs.umn.edu/courses/katagiri/lecture_notes/081103Lec_Statistics3.pdf | ||
* If I allow 10% FDR, there are 3,202 CpG-SNP associations with the minimum p-value | * If I allow 10% FDR, there are 3,202 CpG-SNP associations with the minimum p-value 5.6E-4 | ||
head -3203 sorted_adjustedPVal_2MB_mQTL_output.txt > 2MB_UCLA_mQTL_3202_10%FDR.txt | head -3203 sorted_adjustedPVal_2MB_mQTL_output.txt > 2MB_UCLA_mQTL_3202_10%FDR.txt | ||
Line 45: | Line 45: | ||
sed s'/:/\t/g' 2MB_UCLA_mQTL_3202_10%FDR.txt | awk '{if ($2 == $5 || $2 == $5-1) print $0;}' > 10%FDR_CpG-disrupted-SNP-2MB_UCLA_mQTL (414 CpG-SNP pairs) | sed s'/:/\t/g' 2MB_UCLA_mQTL_3202_10%FDR.txt | awk '{if ($2 == $5 || $2 == $5-1) print $0;}' > 10%FDR_CpG-disrupted-SNP-2MB_UCLA_mQTL (414 CpG-SNP pairs) | ||
sed s'/:/\t/g' 2MB_UCLA_mQTL_3202_10%FDR.txt | awk '{if ($2 != $5 && $2 != $5-1) print $0;}' > 10%FDR_noCpG-disrupted-SNP-2MB_UCLA_mQTL (2789 CpG-SNP pairs) | sed s'/:/\t/g' 2MB_UCLA_mQTL_3202_10%FDR.txt | awk '{if ($2 != $5 && $2 != $5-1) print $0;}' > 10%FDR_noCpG-disrupted-SNP-2MB_UCLA_mQTL (2789 CpG-SNP pairs) | ||
* For the sites that have adjusted p-value (from mQTL analysis) < 0.05, I listed all those sites in separate file and also classified CpG-SNP pair that showed association of CpG and SNP on their own sites. It turned out that all CpG that associate with SNP on their site are the association that pass the threshold (FDR 10%, significant p-value 5. | * For the sites that have adjusted p-value (from mQTL analysis) < 0.05, I listed all those sites in separate file and also classified CpG-SNP pair that showed association of CpG and SNP on their own sites. It turned out that all CpG that associate with SNP on their site are the association that pass the threshold '''(FDR 10%, significant p-value 0.0005448 or 5.448E-4).''' | ||
head -12928 sorted_adjustedPVal_2MB_mQTL_output.txt > sorted_adjustedPVal0.05_2MB_mQTL.txt | head -12928 sorted_adjustedPVal_2MB_mQTL_output.txt > sorted_adjustedPVal0.05_2MB_mQTL.txt | ||
sed s'/:/\t/g' sorted_adjustedPVal0.05_2MB_mQTL.txt |awk '{if ($2 == $5 || $2 == $5-1) print $0;}' > adjPVal0.05_CpG-disrupted-SNP-2MB_UCLA_mQTL (414) | sed s'/:/\t/g' sorted_adjustedPVal0.05_2MB_mQTL.txt |awk '{if ($2 == $5 || $2 == $5-1) print $0;}' > adjPVal0.05_CpG-disrupted-SNP-2MB_UCLA_mQTL (414) | ||
sed s'/:/\t/g' sorted_adjustedPVal0.05_2MB_mQTL.txt |awk '{if ($2 != $5 && $2 != $5-1) print $0;}' > adjPVal0.05_no-CpG-disrupted-SNP-2MB_UCLA_mQTL (12,514) | sed s'/:/\t/g' sorted_adjustedPVal0.05_2MB_mQTL.txt |awk '{if ($2 != $5 && $2 != $5-1) print $0;}' > adjPVal0.05_no-CpG-disrupted-SNP-2MB_UCLA_mQTL (12,514) | ||
* Note: Actually, I no need to repeat mQTL again since I just find another command that I just put the number of test (p.adjust) instead of providing all tests ( mt.rawp2adjp). | * Note: Actually, I no need to repeat mQTL again since I just find another command that I just put the number of test (p.adjust) instead of providing all tests ( mt.rawp2adjp). | ||
* I plotted manhattan plots by separating the CpG-SNP based on the distance between CpG and SNP pairs. Note that I used all the SNP sites that showed association and have adjust p-value (from mQTL analysis) < 0.05. Then use the suggestive line to indicate the threshold at 10%FDR (equivalent to adjust p-value ~ | * Note again: '''FDR 10%, significant p-value 0.0005448 or 5.448E-4.''' | ||
* I plotted manhattan plots by separating the CpG-SNP based on the distance between CpG and SNP pairs. Note that I used all the SNP sites that showed association and have adjust p-value (from mQTL analysis) < 0.05. Then use the suggestive line to indicate the threshold at 10%FDR (equivalent to adjust p-value ~5.448E-4 ). Also The number of CpG-SNP pairs with 10%FDR have been listed on the plots. All the plots should look like the plots I previously did, except that the threshold line was change since I used 10% FDR as the bottom line. | |||
{| {{table}} border = 1 | {| {{table}} border = 1 | ||
| align="center" style="background:#f0f0f0;"|'''CpG-SNP distance''' | | align="center" style="background:#f0f0f0;"|'''CpG-SNP distance''' | ||
Line 56: | Line 57: | ||
| align="center" style="background:#f0f0f0;"|'''# of association, 10%FDR <br>after multiple testing correction''' | | align="center" style="background:#f0f0f0;"|'''# of association, 10%FDR <br>after multiple testing correction''' | ||
|- | |- | ||
| CpG-disrupted SNP||414||414 | | CpG-disrupted SNP||414||414 (12.9%)<br> the rest below is 87.1% | ||
|- | |- | ||
| 2kb||1,433||922 | | 2kb||1,433||922 (28.8%) | ||
|- | |- | ||
| 2-10kb||665||380 | | 2-10kb||665||380 (11.9%) | ||
|- | |- | ||
| 10-50kb||1,443||622 | | 10-50kb||1,443||622 (19.4%) | ||
|- | |- | ||
| 50kb-2MB||8,972||864 | | 50kb-2MB||8,972||864 ((27.0%) | ||
|- | |- | ||
| Total||12,927||3,202 | | Total||12,927||3,202 | ||
Line 74: | Line 75: | ||
[[File:adjPVal0.05_no-CpG-disrupted_10-50KB_mh.pdf]] | [[File:adjPVal0.05_no-CpG-disrupted_10-50KB_mh.pdf]] | ||
[[File:adjPVal0.05_no-CpG-disrupted_50KB-2MB_mh.pdf]] | [[File:adjPVal0.05_no-CpG-disrupted_50KB-2MB_mh.pdf]] | ||
* Continued on: http://genome-tech.ucsd.edu/LabNotes/index.php/Noi/NOTES/2012-4-26 | |||
* Combined the sorted (p-value) output of mQTL and the q-value from multtest | |||
paste -d"\t" sorted_adjustedPVal_2MB_mQTL_output.txt 2MB_UCLA_mQTL_p.adjust_pval_BH.txt > sorted_adjustedPVal_2MB_mQTL_pval_Qval.txt |
Latest revision as of 23:37, 26 April 2012
- Link to calendar: [[1]]
mQTL analysis of UCLA SZ data set[edit]
- Because there are many issues have been discussed and commented about the result from mQTL analysis including p-value cutoff, and the way I characterized and classified CpG sites especially CpG-SNP.
- From the previous analysis: http://genome-tech.ucsd.edu/LabNotes/index.php/Noi/NOTES/2012-1-24, I used the script to run mQTL and allowed only the associated CpG-SNP pairs that have adjusted p-val < 0.05. Total I got 12,932 pairs of CpG-SNP association from total 582,597 CpG-SNP pairs have been tested. After calculating p-value cutoff using Bonferroni correction (p-value/total tested CpG/SNP), there were 1,393 CpG-SNP pairs passed the threshold with p-value 8.58E-8 . By using Bonferroni correction, this method is too stringent, so many true positive signals might be removed even it allows low false positive. Since Dr. Zhang suggested to use multtest function of R to calculate FDR and selected the significant CpG-SNP association by the FDR value. I edited the mQTL script to allow all the CpG-SNP and adjusted p-valued to be printed out. Since to calculate FDR by Benjamini-Hochberg (BH) FDR, all the tested samples and p-value must be listed. The reason that I used Benjamini-Hochberg for multiple test correction are it is not as strict as Bonferroni even still conservative, and it allows the estimation of type I error among the significant results we found . It's more powerful than the other methods.
- Working directory: genome-miner
- /home/nplongth/Noi_scratch/Data_analysis/UCLA_SZ_data_analysis_2012_04_05/mQTL_UCLA_SZ
I generated 8 folders of 2MB.UCLA.mQTL_aa - 2MB.UCLA.mQTL_ah nohup ../mQTL_2MB_plink_allP.pl ../UCLA.hg19SNP ../methylSplit_UCLA.aa >2MB.UCLA.mQTL_aa nohup ../mQTL_2MB_plink_allP.pl ../UCLA.hg19SNP ../methylSplit_UCLA.ab >2MB.UCLA.mQTL_ab nohup ../mQTL_2MB_plink_allP.pl ../UCLA.hg19SNP ../methylSplit_UCLA.ac >2MB.UCLA.mQTL_ac nohup ../mQTL_2MB_plink_allP.pl ../UCLA.hg19SNP ../methylSplit_UCLA.ad >2MB.UCLA.mQTL_ad nohup ../mQTL_2MB_plink_allP.pl ../UCLA.hg19SNP ../methylSplit_UCLA.ae >2MB.UCLA.mQTL_ae nohup ../mQTL_2MB_plink_allP.pl ../UCLA.hg19SNP ../methylSplit_UCLA.af >2MB.UCLA.mQTL_af nohup ../mQTL_2MB_plink_allP.pl ../UCLA.hg19SNP ../methylSplit_UCLA.ag >2MB.UCLA.mQTL_ag nohup ../mQTL_2MB_plink_allP.pl ../UCLA.hg19SNP ../methylSplit_UCLA.ah >2MB.UCLA.mQTL_ah - Concatenate all output from mQTL cat UCLA*/2MB* > combined_2MB awk '{if ($1 ~ /chr/) print $0}' combined_2MB > 2MB_mQTL_output.txt header: CpG_pos SNP_ID Chr(SNP) pos(SNP) R p_val adjp_val total number of tests: 582,597 sort -g -k7 2MB_mQTL_output.txt > sorted_adjustedPVal_2MB_mQTL_output.txt if count only the sites that have adjusted p-value < 0.05: awk '{if ($7 < 0.05) print $0}' | wc -l --> 12,927 CpG-SNP pairs ($7 = adjusted p-value)
Multiple testing correction[edit]
- R "multtest"
source("http://bioconductor.org/biocLite.R") biocLite("multtest") library(multtest) A=read.table("sorted_adjustedPVal_2MB_mQTL_output.txt",header=TRUE) B=A$adjp_val p_adj <- mt.rawp2adjp(B, proc="BH", alpha = 0.05) summary(p_adj) write.table(file="2MB_UCLA_mQTL_mt.rawp2adjp_pval_BH.txt", p_adj$adjp, append = FALSE,row.names = FALSE)
#an alternative command: p_adj=p.adjust(B, method = "BH", n = length(B)) A=read.table("sorted_adjustedPVal_2MB_mQTL_output.txt",header=TRUE) B=A$adjp_val p_adj=p.adjust(B, method = "BH", n = length(B)) write.table(file="2MB_UCLA_mQTL_p.adjust_pval_BH.txt", p_adj, append = FALSE,row.names = FALSE) Note: Both commands gave the same results of q-value or FDR. Length could be set longer than the length of input p-value.
- Links: http://rss.acs.unt.edu/Rdoc/library/stats/html/p.adjust.html, http://www.cbs.umn.edu/courses/katagiri/lecture_notes/081103Lec_Statistics3.pdf
* If I allow 10% FDR, there are 3,202 CpG-SNP associations with the minimum p-value 5.6E-4 head -3203 sorted_adjustedPVal_2MB_mQTL_output.txt > 2MB_UCLA_mQTL_3202_10%FDR.txt
- Count CpG-SNP pair that showed association of CpG and SNP on their own sites
sed s'/:/\t/g' 2MB_UCLA_mQTL_3202_10%FDR.txt | awk '{if ($2 == $5 || $2 == $5-1) print $0;}' > 10%FDR_CpG-disrupted-SNP-2MB_UCLA_mQTL (414 CpG-SNP pairs) sed s'/:/\t/g' 2MB_UCLA_mQTL_3202_10%FDR.txt | awk '{if ($2 != $5 && $2 != $5-1) print $0;}' > 10%FDR_noCpG-disrupted-SNP-2MB_UCLA_mQTL (2789 CpG-SNP pairs)
- For the sites that have adjusted p-value (from mQTL analysis) < 0.05, I listed all those sites in separate file and also classified CpG-SNP pair that showed association of CpG and SNP on their own sites. It turned out that all CpG that associate with SNP on their site are the association that pass the threshold (FDR 10%, significant p-value 0.0005448 or 5.448E-4).
head -12928 sorted_adjustedPVal_2MB_mQTL_output.txt > sorted_adjustedPVal0.05_2MB_mQTL.txt sed s'/:/\t/g' sorted_adjustedPVal0.05_2MB_mQTL.txt |awk '{if ($2 == $5 || $2 == $5-1) print $0;}' > adjPVal0.05_CpG-disrupted-SNP-2MB_UCLA_mQTL (414) sed s'/:/\t/g' sorted_adjustedPVal0.05_2MB_mQTL.txt |awk '{if ($2 != $5 && $2 != $5-1) print $0;}' > adjPVal0.05_no-CpG-disrupted-SNP-2MB_UCLA_mQTL (12,514)
- Note: Actually, I no need to repeat mQTL again since I just find another command that I just put the number of test (p.adjust) instead of providing all tests ( mt.rawp2adjp).
- Note again: FDR 10%, significant p-value 0.0005448 or 5.448E-4.
- I plotted manhattan plots by separating the CpG-SNP based on the distance between CpG and SNP pairs. Note that I used all the SNP sites that showed association and have adjust p-value (from mQTL analysis) < 0.05. Then use the suggestive line to indicate the threshold at 10%FDR (equivalent to adjust p-value ~5.448E-4 ). Also The number of CpG-SNP pairs with 10%FDR have been listed on the plots. All the plots should look like the plots I previously did, except that the threshold line was change since I used 10% FDR as the bottom line.
CpG-SNP distance | # of association, p-value < 0.05 before multiple testing correction |
# of association, 10%FDR after multiple testing correction |
CpG-disrupted SNP | 414 | 414 (12.9%) the rest below is 87.1% |
2kb | 1,433 | 922 (28.8%) |
2-10kb | 665 | 380 (11.9%) |
10-50kb | 1,443 | 622 (19.4%) |
50kb-2MB | 8,972 | 864 ((27.0%) |
Total | 12,927 | 3,202 |
File:AdjustedPVal0.05 2MB UCLA mQTL mh.txt.pdf File:AdjPVal0.05 CpG-disrupted-SNP-2MB UCLA mQTL mh.pdf File:AdjPVal0.05 no-CpG-disrupted 2KB mh.pdf File:AdjPVal0.05 no-CpG-disrupted 2-10KB mh.pdf File:AdjPVal0.05 no-CpG-disrupted 10-50KB mh.pdf File:AdjPVal0.05 no-CpG-disrupted 50KB-2MB mh.pdf
- Continued on: http://genome-tech.ucsd.edu/LabNotes/index.php/Noi/NOTES/2012-4-26
- Combined the sorted (p-value) output of mQTL and the q-value from multtest
paste -d"\t" sorted_adjustedPVal_2MB_mQTL_output.txt 2MB_UCLA_mQTL_p.adjust_pval_BH.txt > sorted_adjustedPVal_2MB_mQTL_pval_Qval.txt