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

From ZhangLabWiki
Jump to navigation Jump to search
>Noi
No edit summary
>Noi
No edit summary
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%).
* 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.6E-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)

Revision as of 16:43, 20 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 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

  • 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.
* If I allow 10% FDR, there are 3,202 CpG-SNP associations with the minimum p-value 0.0005448
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 5.6E-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).
  • 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 ~ 0.0005448). 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
2kb 1,433 922
2-10kb 665 380
10-50kb 1,443 622
50kb-2MB 8,972 864
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