mQTL analysis on UCLA-SZ data set (no CpG-SNP)
- Details will be filled later
- Filter out SNP from methylMatrix (Note: the sample IDs were modified to be compatible with tfam file and plink)
./remove_snps.pl ./UCLA-4batches_methyl_min76_minSTD0.1-edit /t4T/GenomeDB/HsGenome/snp134_snv.txt > UCLA-4batches_methyl-noSNP_min76_minSTD0.1
61664 UCLA-4batches_methyl-noSNP_min76_minSTD0.1
- Split the methylMatrix in to 10,000 CpG sites each
split -10000 UCLA-4batches_methyl-noSNP_min76_minSTD0.1 UCLA_noSNP_
10000 UCLA_noSNP_aa
10001 UCLA_noSNP_ab
10001 UCLA_noSNP_ac
10001 UCLA_noSNP_ad
10001 UCLA_noSNP_ae
10001 UCLA_noSNP_af
1666 UCLA_noSNP_ag
nohup ../mQTL_2MB_plink_allP.pl ../UCLA.hg19SNP ../UCLA_noSNP_aa > 2MB_noNSP_aa
nohup ../mQTL_2MB_plink_allP.pl ../UCLA.hg19SNP ../UCLA_noSNP_ab > 2MB_noNSP_ab
nohup ../mQTL_2MB_plink_allP.pl ../UCLA.hg19SNP ../UCLA_noSNP_ac > 2MB_noNSP_ac
nohup ../mQTL_2MB_plink_allP.pl ../UCLA.hg19SNP ../UCLA_noSNP_ad > 2MB_noNSP_ad
nohup ../mQTL_2MB_plink_allP.pl ../UCLA.hg19SNP ../UCLA_noSNP_ae > 2MB_noNSP_ae
nohup ../mQTL_2MB_plink_allP.pl ../UCLA.hg19SNP ../UCLA_noSNP_af > 2MB_noNSP_af
nohup ../mQTL_2MB_plink_allP.pl ../UCLA.hg19SNP ../UCLA_noSNP_ag > 2MB_noNSP_ag
- Concatenate the out put files to the same file
cat mQTL_noSNP_a*/2MB_noNSP_a* > 2MB_mQTL-noSNP_rawoutput
awk '{if ($1 ~ /chr/) print $0}' 2MB_mQTL-noSNP_rawoutput > 2MB_mQTL-noSNP_output
433,947 2MB_mQTL-noSNP_output
- Sort out adj. p-val and count the number of significant associated CpG-SNP roughly at different p-value
sort column7 (adj p_val)
sort -g -k7 2MB_mQTL-noSNP_output | awk '{if ($7<=0.05) print $0}' | wc -l --> 6,418
sort -g -k7 2MB_mQTL-noSNP_output | awk '{if ($7<=0.01) print $0}' | wc -l --> 2,292
sort -g -k7 2MB_mQTL-noSNP_output | awk '{if ($7<=0.001) print $0}' | wc -l --> 1,026
sort -g -k7 2MB_mQTL-noSNP_output | awk '{if ($7<=0.0001) print $0}' | wc -l --> 698
sort -g -k7 2MB_mQTL-noSNP_output | awk '{if ($7<=0.00001) print $0}' | wc -l --> 526
- Print out the site at min p-value 0.05 before multiple test correction
sort -g -k7 2MB_mQTL-noSNP_output | awk '{if ($7<=0.05) print $0}' > sorted_2MB_mQTL-noSNP_minPVal0.05.txt
6418 sorted_2MB_mQTL-noSNP_minPVal0.05.txt
- Concatenate QVal to the mQTL file
paste -d"\t" sorted_2MB_mQTL-noSNP_minPVal0.05.txt sorted_2MB_mQTL-noSNP_QVal.txt.txt > sorted_2MB_mQTL-noSNP_minPVal0.05_QVal.txt.txt
awk '{if ($8 <=0.1) print $0}' sorted_2MB_mQTL-noSNP_minPVal0.05_QVal.txt.txt > 10%FDR_2MB_mQTL-noSNP
741 10%FDR_2MB_mQTL-noSNP
p-value cut off: 0.0001636 or 1.636E-4
- Look at # of associated CpG
awk '{print $1}' 10%FDR_2MB_mQTL-noSNP | sort | uniq -u | wc -l 295
awk '{print $1}' 10%FDR_2MB_mQTL-noSNP | sort | uniq -d | wc -l 121
295+121X = 741 --> X ~ (741-295)/121 ~3.69 this show >1 CpGs associated with the same SNP
Total single CpG: 295+121 = 416
- Look at # of associated SNP
awk '{print $2}' 10%FDR_2MB_mQTL-noSNP | sort | uniq -u | wc -l --> 195
awk '{print $2}' 10%FDR_2MB_mQTL-noSNP | sort | uniq -d | wc -l --> 109
195+109Y = 741 --> Y ~ (741-195)/109 ~ 5 this also suggest >1 SNP associated with the same CpG
Total single SNP 195+109 = 304
- Separate the CpG-SNP association based on the distance between CpG and SNP
- Note: Dr. Zhang suggested to break down the CpG-SNP distance with equal length
sed s'/:/\t/g' 10%FDR_2MB_mQTL-noSNP > 10%FDR_2MB_mQTL-noSNP_split_pos
Run sh print_distance.sh
- Note again: P-value of 10%FDR is p <= 1.636E−4
Associated CpG-SNP distance
|
Number of assiciated CpG-SNP
|
%
|
10%FDR_10kb |
378 |
50.94%
|
10%FDR_10-20kb |
27 |
3.64%
|
10%FDR_20-30kb |
21 |
2.83%
|
10%FDR_30-40kb |
36 |
4.85%
|
10%FDR_40-50kb |
2 |
0.27%
|
10%FDR_50kb-2Mb |
278 |
37.47%
|
Total in 2Mb |
742 |
|
Associated CpG-SNP distance
|
Number of assiciated CpG-SNP
|
%
|
minpval0.05_10kb |
862 |
13.43%
|
minpval0.05_10-20kb |
124 |
1.93%
|
minpval0.05_20-30kb |
133 |
2.07%
|
minpval0.05_30-40kb |
108 |
1.68%
|
minpval0.05_40-50kb |
75 |
1.17%
|
minpval0.05_50kb-2Mb |
5,116 |
79.71%
|
Total in 2Mb |
6,418 |
|
- Get sequences for manhattan plots
Run get_manht_seq.sh
Comparing no CpG-SNP mQTL to the whole CpG and SNP mQTL result
- Analyzing the association of CpG and SNP with 10%FDR derived from each analysis
'
|
no CpG-SNP
|
SNP-free CpG and CpG-SNP
|
% of no CpG-SNP to all significant associations
|
# of association (10%FDR) |
742 |
3,202 |
23.17%
|
significant p-value |
1.636E-04 |
5.448E-04 |
|
- From the table above, I could see the different of significant p-value from each test. The reason is that the number of tests were different between the test of the whole methylation matrix and the whose SNPs have been filtered out. Therefore the BH correction which deal with the number of test in calculation provide the q-value with at the different significant p-value. I turned out that the significant p-value of no SNP analysis was lower than the whole CpG. This implied that there might be more no CpG-SNP between 1.636E-04 to 5.448E-04 p-value of the 2MB_UCLA_mQTL_3202_10%FDR.txt file (Total significant CpG-SNP association 3,202 sites, 10%FDR)
- I then checked the number of CpG-SNP association with the significant p-value 5.448E-04 in no SNP-CpG file and got total 900 sites (but have ~27.4% FDR by the test with out CpG-SNP.
- Note: I can count the # of sites based on p-value because this value was not changed, and just used it as the reference. This account for ~28% of the association I got from all association.
- Then, I printed out another files of CpG-SNP association with no CpG-SNP significant p-value 0.0005613
awk '{if ($7<=0.0005613) print $0}' sorted_2MB_mQTL-noSNP_minPVal0.05_QVal.txt.txt > sorted_2MB_mQTL-noSNP_minPVal0.0005613_QVal.txt
- I also checked if those sites exit in the 10%FDR mQTL result derived from the previous analysis (on the whole mathylMatrix)
./get_overlapping_CpG_SNP.pl ./2MB_UCLA_mQTL_3202_10%FDR.txt <./sorted_2MB_mQTL-noSNP_minPVal0.0005613_QVal.txt > test_mQTL
- I checked the overlapped sites on both files again, and there were consistent. This helped to confirm that the analysis of this time and the previous time were likely to be correct.
head sorted_2MB_mQTL-noSNP_minPVal
chr8:17533620 rs17125083 8 17533468 0.8961 1.331e-40 2.262e-39
chr8:17533620 rs2301541 8 17533568 0.9029 3.106e-40 2.64e-39
chr8:17533604 rs17125083 8 17533468 0.8839 1.099e-38 1.867e-37
chr19:58791188 rs260500 19 58791213 0.889 5.069e-38 6.083e-37
chr8:17533604 rs2301541 8 17533568 0.8859 1.473e-37 1.252e-36
chr19:58791159 rs260500 19 58791213 0.8844 2.421e-37 2.906e-36
head sorted_2MB_mQTL-noSNP_minPVal0.0005613_QVal.txt
chr8:17533620 rs17125083 8 17533468 0.8961 1.331e-40 2.262e-39 5.9861604e-34
chr8:17533620 rs2301541 8 17533568 0.9029 3.106e-40 2.64e-39 5.9861604e-34
chr8:17533604 rs17125083 8 17533468 0.8839 1.099e-38 1.867e-37 2.82226299666667e-32
chr19:58791188 rs260500 19 58791213 0.889 5.069e-38 6.083e-37 6.8965556275e-32
chr8:17533604 rs2301541 8 17533568 0.8859 1.473e-37 1.252e-36 1.135556488e-31
chr19:58791159 rs260500 19 58791213 0.8844 2.421e-37 2.906e-36 2.19643713666667e-31
- Also break down CpG-SNP distance not only at equal length, but also at a very short distance like within 2kb and 5kb. The table below implied the cis regulation of genetic factor on methylation status especially within 2kb since the number of association not largely increase from 2-5kb (393-327 = 66 association) or 2-10kb (426-327 = 99 association)
CpG-SNP distance
|
number of association
|
%
|
2kb-noSNP |
327 |
36.33%
|
5kb-noSNP |
393 |
43.67%
|
10kb-noSNP |
426 |
47.33%
|
10-20kb-noSNP |
34 |
3.78%
|
20-30kb-noSNP |
24 |
2.67%
|
30-40kb-noSNP |
41 |
4.56%
|
40-50kb-noSNP |
4 |
0.44%
|
50kb-2Mb-noSNP |
371 |
41.22%
|