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

From ZhangLabWiki
Jump to navigation Jump to search
>Noi
No edit summary
>Noi
No edit summary
Line 168: Line 168:
  awk '{print $2}' sorted_2MB_mQTL-noSNP_minPVal0.0005613_QVal.txt | sort | uniq -d | wc -l --> 134
  awk '{print $2}' sorted_2MB_mQTL-noSNP_minPVal0.0005613_QVal.txt | sort | uniq -d | wc -l --> 134
  236+134Y = 900 --> Y ~ (900-236)/134 ~ 5 this suggests that > 1 CpGs associated with the same SNP (134 of individual SNPs associated with more than one CpGs)
  236+134Y = 900 --> Y ~ (900-236)/134 ~ 5 this suggests that > 1 CpGs associated with the same SNP (134 of individual SNPs associated with more than one CpGs)
* Print out the list of unique CpGs and SNPs and duplicated CpGs and SNPs
awk '{print $1}' sorted_2MB_mQTL-noSNP_minPVal0.0005613_QVal.txt | sort | uniq -u > uniqCpG
awk '{print $1}' sorted_2MB_mQTL-noSNP_minPVal0.0005613_QVal.txt | sort | uniq -d > duplicatedCpG
awk '{print $2}' sorted_2MB_mQTL-noSNP_minPVal0.0005613_QVal.txt | sort | uniq -u > uniqSNP
awk '{print $2}' sorted_2MB_mQTL-noSNP_minPVal0.0005613_QVal.txt | sort | uniq -d > duplicatedSNP

Revision as of 20:37, 27 April 2012

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
  • Run mQTL (no CpG-SNP)
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
  • Number of test
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 SNPs associated with the same CpG
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 CpGs associated with the same SNP
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%
  • Look at # of associated CpG
900 sorted_2MB_mQTL-noSNP_minPVal0.0005613_QVal.txt
awk '{print $1}' sorted_2MB_mQTL-noSNP_minPVal0.0005613_QVal.txt | sort | uniq -u | wc -l --> 361
awk '{print $1}' sorted_2MB_mQTL-noSNP_minPVal0.0005613_QVal.txt | sort | uniq -d | wc -l --> 155
361+155X = 900 --> X ~ (900-361)/155 ~ 3.5 this shows that >1 SNPs associated with the same CpG (155 of individual CpGs associated with more than one SNPs)

Total single CpG: 361+155 = 516

  • Look at # of associated CpG
awk '{print $2}' sorted_2MB_mQTL-noSNP_minPVal0.0005613_QVal.txt | sort | uniq -u | wc -l --> 236
awk '{print $2}' sorted_2MB_mQTL-noSNP_minPVal0.0005613_QVal.txt | sort | uniq -d | wc -l --> 134
236+134Y = 900 --> Y ~ (900-236)/134 ~ 5 this suggests that > 1 CpGs associated with the same SNP (134 of individual SNPs associated with more than one CpGs)
  • Print out the list of unique CpGs and SNPs and duplicated CpGs and SNPs
awk '{print $1}' sorted_2MB_mQTL-noSNP_minPVal0.0005613_QVal.txt | sort | uniq -u > uniqCpG
awk '{print $1}' sorted_2MB_mQTL-noSNP_minPVal0.0005613_QVal.txt | sort | uniq -d > duplicatedCpG
awk '{print $2}' sorted_2MB_mQTL-noSNP_minPVal0.0005613_QVal.txt | sort | uniq -u > uniqSNP
awk '{print $2}' sorted_2MB_mQTL-noSNP_minPVal0.0005613_QVal.txt | sort | uniq -d > duplicatedSNP