Noi/NOTES/2012-4-26

From ZhangLabWiki
Revision as of 22:38, 26 April 2012 by >Noi
Jump to navigation Jump to search

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
294+121X = 741 --> X ~ (741-295)/121 ~3.69 this show >1 CpGs associated with the same SNP
Total single CpG: 294+121 = 415
*Look at # of associated  SNP
awk '{print $2}' 10%FDR_2MB_mQTL-noSNP | sort | uniq -u | wc -l  --> 195 (subtract header line)
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