Noi/NOTES/2012-1-24
UCLA (Schizophrenia) data analysis (continued)
mQTL analysis
- Dr. Zhang suggested that we should start to do mQTL analysis with the SNP data we called by look at 2M or 5M window.
- Dinh's link for mQTL analysis on UPenn Data: [[1]]
- 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
- methylSplit_UCLA.aa
- methylSplit_UCLA.ah
- methylSplit_UCLA.ac
- methylSplit_UCLA.ad
- methylSplit_UCLA.ae
- methylSplit_UCLA.af
- methylSplit_UCLA.ag
- 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.
=========
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.
- 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)
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
- Print out only the SNP sites
awk '{if ($0 ~ /chr/) print $0;}' 2MB.UCLA.mQTL_aa.txt > 2MB.UCLA.mQTL-SNP_aa.txt
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
../../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