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 49: Line 49:
  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)
* 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.
* 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.
{| {{table}} border = 1
| align="center" style="background:#f0f0f0;"|'''mQTL results'''
| align="center" style="background:#f0f0f0;"|'''# of associated CpG-SNP pairs,p-value < 0.05'''
| align="center" style="background:#f0f0f0;"|'''# of associated CpG-SNP pairs, 10%FDR'''
|-
| 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
|}

Revision as of 20:30, 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 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)
  • 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%).
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)
  • 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.
mQTL results # of associated CpG-SNP pairs,p-value < 0.05 # of associated CpG-SNP pairs, 10%FDR
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