Noi/NOTES/2012-4-11: Difference between revisions

From ZhangLabWiki
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)