Noi/NOTES/2012-4-11: Difference between revisions
Jump to navigation
Jump to search
>Noi No edit summary |
>Noi |
||
Line 2: | Line 2: | ||
== mQTL analysis of UCLA SZ data set == | == mQTL analysis of UCLA SZ data set == | ||
* 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. | * 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. 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 . | * 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. 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 | * Working directory: genome-miner | ||
* /home/nplongth/Noi_scratch/Data_analysis/UCLA_SZ_data_analysis_2012_04_05/mQTL_UCLA_SZ | * /home/nplongth/Noi_scratch/Data_analysis/UCLA_SZ_data_analysis_2012_04_05/mQTL_UCLA_SZ |
Revision as of 09:36, 12 April 2012
- Link to calendar: [[1]]
mQTL analysis of UCLA SZ data set
- 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. 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)
- 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)