Noi/NOTES/2012-1-24

From ZhangLabWiki
Jump to navigation Jump to search

Link to calendar:[[1]]

UCLA (Schizophrenia) data analysis (continued)[edit]

MethylFreq Matrix[edit]

  • I have generated metyhlFreq matrix of total 96 samples which all showed methylFreq Corr > 0.900 by allowing 20% of missing values (min samples = 76)

./allBED2Matrix_DD.pl bed.list 76 0 0.1, output file: UCLA-4batches_methyl_min76_minSTD0.1

  • The matrix contains 79,604 variable CpG sites
  • Directory: /home/nplongth/Noi_scratch/combined_BSPP_UCLA-4batches_20120103/BSPP_UCLA.hg19.mapping-batch4

mQTL analysis (discarded this analysis after mQTL run)[edit]

  • Dr. Zhang suggested that we should start to do mQTL analysis with the SNP data we called by looking at 2M or 5M window.
  • Dinh's link for mQTL analysis on UPenn Data: [[2]]
  • Dinh: /home/dinh/UPenn_Analysis/UPenn48_April12/dbSNP130_hg18_filteredSNPs/

1. All .snp files were filtered with hg19 dbSNP using script written by Dr. Zhang in genome-miner: /home/kunzhang/bin/bisSnpFilterHg19.pl
snpFilter.sh

for f in *.snp
do

/home/kunzhang/bin/bisSnpFilterHg19.pl $f > $f.filtered
done;

2. There are 4 files that were mapped separately since the .fastq files are very large. Dinh helped me wrote script to merge the SNP data.

  • To merge and sort the chromosome position: cat GK* | sort -k1,1 > all.GK.snp.filtered
  • mergeAndSplit.pl: Media:mergeAndSplit-e.txt

./mergeAndSplit.pl < all-GK0091-003.snp.filtered > GK0091-003.fwd.pileup.snp.filtered
./mergeAndSplit.pl < all-GK0277-002.snp.filtered > GK0277-002.fwd.pileup.snp.filtered
./mergeAndSplit.pl < all-GK0277-005.snp.filtered > GK0277-005.fwd.pileup.snp.filtered
./mergeAndSplit.pl < all-GK0123-002.snp.filtered > GK0123-002.fwd.pileup.snp.filtered
3. Generated tped and tfam files: GetTPED.txtMedia:GetTPED.txt
ls UCLA_filtered_SNP/*snp.filtered > list_snp
awk '{print "123"$0"123\t"$0;}' list_snp | sed 's/123UCLA_filtered_SNP\///g' | sed 's/.fwd.pileup.snp.filtered123//g' > filtered.snp_list

  • Still have some problem when running plink. The error message said about the number of column didn't match with the expected number of column. Need to figure out why it expected the wrong number.
  • I test with smaller number of samples in the matrix to see if the script work.
  • Finally, I could figure out that because the header of sample was not list properly since there were no tabs between some sample IDs.
  • Dinh removed the - sign from the sample IDs and I edited the whole methylFreq matrix before split the matrix and repeated running mQTL again.

split -l 10000 UCLA-4batches_methyl_min76_minSTD0.1-edit methylSplit_UCLA.

  • Total line of methylfreq matrix: 79605, after split, there were 8 files
  1. methylSplit_UCLA.aa
  2. methylSplit_UCLA.ah
  3. methylSplit_UCLA.ac
  4. methylSplit_UCLA.ad
  5. methylSplit_UCLA.ae
  6. methylSplit_UCLA.af
  7. methylSplit_UCLA.ag
  8. methylSplit_UCLA.ah

../mQTL_2MB_plink.pl ../UCLA.hg19snp ../methylSplit_UCLA.aa > UCLA.mQTL_aa.txt --> output in folder mQTL-aa
../mQTL_2MB_plink.pl ../UCLA.hg19snp ../methylSplit_UCLA.ab > UCLA.mQTL_ab.txt --> output in folder mQTL-ab
../mQTL_2MB_plink.pl ../UCLA.hg19snp ../methylSplit_UCLA.ac > UCLA.mQTL_ac.txt --> output in folder mQTL-ac
../mQTL_2MB_plink.pl ../UCLA.hg19snp ../methylSplit_UCLA.ad > UCLA.mQTL_ad.txt --> output in folder mQTL-ad
../mQTL_2MB_plink.pl ../UCLA.hg19snp ../methylSplit_UCLA.ae > UCLA.mQTL_ae.txt --> output in folder mQTL-ae
../mQTL_2MB_plink.pl ../UCLA.hg19snp ../methylSplit_UCLA.af > UCLA.mQTL_af.txt --> output in folder mQTL-af
../mQTL_2MB_plink.pl ../UCLA.hg19snp ../methylSplit_UCLA.ag > UCLA.mQTL_ag.txt --> output in folder mQTL-ag
../mQTL_2MB_plink.pl ../UCLA.hg19snp ../methylSplit_UCLA.ah > UCLA.mQTL_ah.txt --> output in folder mQTL-ah
also run the whole matrix since I think it should not take too much time. ../mQTL_2MB_plink.pl ../UCLA.hg19snp ../UCLA-4batches_methyl_min76_minSTD0.1-edit > UCLA.mQTL.txt --> output in folder mQTL-wholematrix

  • Did the same analysis, but increasing the window to 5MB.

../../mQTL_5MB_plink.pl ../../UCLA.hg19snp ../../methylSplit_UCLA.aa > 5MB.UCLA.mQTL_aa.txt
../../mQTL_5MB_plink.pl ../../UCLA.hg19snp ../../methylSplit_UCLA.ab > 5MB.UCLA.mQTL_ab.txt
../../mQTL_5MB_plink.pl ../../UCLA.hg19snp ../../methylSplit_UCLA.ac > 5MB.UCLA.mQTL_ac.txt
../../mQTL_5MB_plink.pl ../../UCLA.hg19snp ../../methylSplit_UCLA.ad > 5MB.UCLA.mQTL_ad.txt
../../mQTL_5MB_plink.pl ../../UCLA.hg19snp ../../methylSplit_UCLA.ae > 5MB.UCLA.mQTL_ae.txt
../../mQTL_5MB_plink.pl ../../UCLA.hg19snp ../../methylSplit_UCLA.af > 5MB.UCLA.mQTL_af.txt
../../mQTL_5MB_plink.pl ../../UCLA.hg19snp ../../methylSplit_UCLA.ag > 5MB.UCLA.mQTL_ag.txt
../../mQTL_5MB_plink.pl ../../UCLA.hg19snp ../../methylSplit_UCLA.ah > 5MB.UCLA.mQTL_ah.txt
../../mQTL_5MB_plink.pl ../../UCLA.hg19snp ../../UCLA-4batches_methyl_min76_minSTD0.1-edit > 5MB.UCLA.mQTL.txt
cat 5MB.UCLA.mQTL_* > merge.5MB.UCLA.mQTL.aa-ah --> mv merge.5MB.UCLA.mQTL.aa-ah merge.5MB.UCLA.mQTL.aa-ah-tmp --> awk '{if ($0 ~ /chr/) print $0;}' merge.5MB.UCLA.mQTL.aa-ah-tmp > merge.5MB.UCLA.mQTL.aa-ah

  • Then sorted out the result based on the p-value.
===================================================================================================================================================================================[edit]

2012_01_30 I will discard the above analysis since there were some thing wrong with the script and data mislabeling correction between methylFreq and SNP data that were not in agreement. Here, I'll add all details of the labeling of the mislabeled data between the two data. Generate tped and tfam files ./GetTPED.pl UCLA.hg19SNP.tfam UCLA.hg19SNP.tped < filtered.snp_list

  • Correct data mislabeling in UCLA.hg19SNP.tfam
    • GK0210-001 --> GK0210-002
    • GK0210-002 --> GK0210-001
    • GK0050-005a --> GK0277-002
    • GK0277-002 --> GK0344-004
    • GK0344-004 --> GK0050-005a
  • Check data labeling in methylFreq matrix
    • GK0210-001 --> GK0210-002
    • GK0210-002 -> GK0210-001
    • GK0050-005a --> GK0277-002
    • GK0277-002 --> GK0344-004
    • GK0344-004 --> GK0050-005a
  • Therefore all data labelings were correct.

Start running mQTL (2MB and 5MB window)[edit]

2MB window : in "/home/nplongth/Noi_scratch/combined_BSPP_UCLA-4batches_20120103/BSPP_UCLA.hg19.mapping-batch4/BSPP_UCLA.hg19.SNP-batch4/2MB_mQTL_UCLA/" ../../mQTL_2MB_plink.pl ../../UCLA.hg19SNP ../../methylSplit_UCLA.aa > 2MB.UCLA.mQTL_aa.txt
../../mQTL_2MB_plink.pl ../../UCLA.hg19SNP ../../methylSplit_UCLA.ab > 2MB.UCLA.mQTL_ab.txt
../../mQTL_2MB_plink.pl ../../UCLA.hg19SNP ../../methylSplit_UCLA.ac > 2MB.UCLA.mQTL_ac.txt
../../mQTL_2MB_plink.pl ../../UCLA.hg19SNP ../../methylSplit_UCLA.ad > 2MB.UCLA.mQTL_ad.txt
../../mQTL_2MB_plink.pl ../../UCLA.hg19SNP ../../methylSplit_UCLA.ae > 2MB.UCLA.mQTL_ae.txt
../../mQTL_2MB_plink.pl ../../UCLA.hg19SNP ../../methylSplit_UCLA.af > 2MB.UCLA.mQTL_af.txt
../../mQTL_2MB_plink.pl ../../UCLA.hg19SNP ../../methylSplit_UCLA.ag > 2MB.UCLA.mQTL_ag.txt
../../mQTL_2MB_plink.pl ../../UCLA.hg19SNP ../../methylSplit_UCLA.ah > 2MB.UCLA.mQTL_ah.txt
../../mQTL_2MB_plink.pl ../../UCLA.hg19SNP ../../UCLA-4batches_methyl_min76_minSTD0.1-edit > 2MB.UCLA.mQTL.txt --> for the whole methylFreq marix

awk '{if ($0 ~ /chr/) print $0;}' 2MB.UCLA.mQTL.txt > 2MB.UCLA.mQTL-SNP
  • To remove CpG site that associated with SNP on it site (not all SNP-containing CPGs were removed):
sed s'/:/\t/g' 2MB.UCLA.mQTL-SNP | awk '{if ($2 == $5 || $2 == $5-1) print $0;}' > CpG-disrupted-SNP-2MB_mQTL (after sed ":" to "\t")

I ended up getting 12,932 site of CpG-SNP association (the minimum P-value 0.05), and after filtered out the SNP-disrupting CpG, I got 12,518 sites.
I performed Manhatton plot for every associated site within 2MB. In addition, I also plotted based on the distance of CpG-SNP in different ranges (2kb, 2-10kb, 10-50kb, and 50kb-2MB) to see the distribution of the associated CpG-SNP with the very low P-value.

  • 2kb: 1,434 CpG-SNP association
  • 2-10KB: 666 CpG-SNP association
  • 10-50kb: 1,446 CpG-SNP association
  • 50kb-2MB: 8,976 CpG-SNP association
    • Note: The manhattan plots before getting the P-value cutoff were removed
  • Plotting manhattan plot:

File:2MB mht UCLA-SZ.png
File:2kb mht UCLA-SZ.png File:2-10kb mht UCLA-SZ.png File:10-50 mht UCLA-SZ.png File:50kb-2Mb mht UCLA-SZ.png

  • From the Manhattan plots, I could observe the very low P-value associated-SNPs fell within 2kb of the associated-CpG sites rather than the long distance associated-CpG.
  • Even though the large fraction of CpG-SNP association fell in between 50kb and 2MB, the P-value seemed to be higher than the sites closed to CpG. I need to do shuffling to find the FDR and get the P-value cut-off to identify the significant association. Finally, Dr. Zhang suggested to used the number of CpG-SNP pairs and P-value to get the P-value cutoff as described below.
  • Then I plotted methylation frequency vs. genotypes by choosing the associated-SNPs that showed different P-values from high to low.

File:UCLA rs4862337 P 1.66E-107.png File:UCLA rs7786917 P 1.37E-60.png File:UCLA rs1716355 P 1.09E-40.png File:UCLA rs4581329 P 1.14E-10.png

  • Dr. Zhang had the questions:
    • How many SNP-CpG pairs have been testes for association?
    • And if the genome-wide threshold for significance with Bonferroni correction is the p_value (such as 0.01) divided by one million ==> 1E-8. So how many CpG-SNP pairs passed this threshold, and also the number from the breakdown based of different distance of CpG-SNP.
  • After extracting the number from mQTL.log file: I got total 582,597 CpG-SNP pairs have been tested. Script to count the # of CpG-SNP pair have been tested from the log file of mQTL analysis(mQTL.log): extracSNPnumber.txt
  • Therefore, the threshold should be 0.05/582,597 = 8.58E-8 and 1,393 paired passed the threshold (in 2MB window).
    • 2kb: 605 pairs passed
    • 2-10kb: 234 pairs passed
    • 10-50kb: 308 pairs passed
    • 50kb-2MB: 246 pairs passed
  • Then Dr. Zhang suggested to identify the CpG associated with SNP from mQTL analysis if some of them exhibit heritable patterns.
    • I extract the CpG-SNP associated sites by comparing to the heritable CpG sites derive from mpo analysis either heritable with or without CpG on SNPs.
    • Note that when I got CpG-SNP associated sites, I removed only the CpG that associated with SNP on its own site. This means that not all CpG on SNPs were removed from the list.
  • Directory: /home/nplongth/Noi_scratch/combined_BSPP_UCLA-4batches_20120103/BSPP_UCLA.hg19.mapping-batch4/BSPP_UCLA.hg19.SNP-batch4/2MB_mQTL_UCLA/mpo-mQTL_analysis
  • extract_heri_CpG-SNP_r2.pl: Media:extract_heri_CpG-SNP_r2.txt
  • sorted_passedP_2MB_UCLA.mQTL.txt (CpG-SNP associated sites passed the P-Value cutoff as described above): sorted_passedP_2MB_UCLA.mQTL.txt]
  • UCLA-RO-mpo_out-heritable.txt: UCLA-RO-mpo_out-heritable.txt
    • Note: Header = chr_pos valid_trios pearson_r2 slope t_test heritable
  • UCLA-RO-mpo_out-heritable-noSNP.txt: UCLA-RO-mpo_out-heritable-noSNP.txt
  • heriatble_CpG-SNP_assoc_UCLA-SZ_r2.txt: heriatble_CpG-SNP_assoc_UCLA-SZ_r2.txt
  • heriatble_CpG-SNP_assoc_noSNP_UCLA-SZ_r2.txt: heriatble_CpG-SNP_assoc_noSNP_UCLA-SZ_r2.txt
  • Note: Header = output from mQTL analysis and the last column is pearson_r2 from mpo analysis
./extract_heri_CpG-SNP_r2.pl sorted_passedP_2MB_UCLA.mQTL < UCLA-RO-mpo_out-heritable.txt > heriatble_CpG-SNP_assoc_UCLA-SZ_r2.txt --> got 610 sites
./extract_heri_CpG-SNP_r2.pl sorted_passedP_2MB_UCLA.mQTL < UCLA-RO-mpo_out-heritable-noSNP.txt > heriatble_CpG-SNP_assoc_noSNP_UCLA-SZ_r2.txt --> got 104 sites
  • I got 104 of CpG associated with SNP (from mQTL) are heritable CpGs, and all are not CpG on SNPs, also got 506 sites are CpG on SNP sites.
===================================================================================================================================================================================[edit]

5MB window: in /home/nplongth/Noi_scratch/combined_BSPP_UCLA-4batches_20120103/BSPP_UCLA.hg19.mapping-batch4/BSPP_UCLA.hg19.SNP-batch4/5MB_mQTL_UCLA

  • Haven't finished summarizing the data

../../mQTL_5MB_plink.pl ../../UCLA.hg19SNP ../../methylSplit_UCLA.aa > 5MB.UCLA.mQTL_aa.txt
../../mQTL_5MB_plink.pl ../../UCLA.hg19SNP ../../methylSplit_UCLA.ab > 5MB.UCLA.mQTL_ab.txt
../../mQTL_5MB_plink.pl ../../UCLA.hg19SNP ../../methylSplit_UCLA.ac > 5MB.UCLA.mQTL_ac.txt
../../mQTL_5MB_plink.pl ../../UCLA.hg19SNP ../../methylSplit_UCLA.ad > 5MB.UCLA.mQTL_ad.txt
../../mQTL_5MB_plink.pl ../../UCLA.hg19SNP ../../methylSplit_UCLA.ae > 5MB.UCLA.mQTL_ae.txt
../../mQTL_5MB_plink.pl ../../UCLA.hg19SNP ../../methylSplit_UCLA.af > 5MB.UCLA.mQTL_af.txt
../../mQTL_5MB_plink.pl ../../UCLA.hg19SNP ../../methylSplit_UCLA.ag > 5MB.UCLA.mQTL_ag.txt
../../mQTL_5MB_plink.pl ../../UCLA.hg19SNP ../../methylSplit_UCLA.ah > 5MB.UCLA.mQTL_ah.txt
../../mQTL_5MB_plink.pl ../../UCLA.hg19SNP ../../UCLA-4batches_methyl_min76_minSTD0.1-edit > 5MB.UCLA.mQTL.txt --> for the whole methylFreq marix