Noi/NOTES/2013-3-26

From ZhangLabWiki
Jump to navigation Jump to search

mQTL analysis on 5M SNPs (UCLA data set)[edit]

  • Actually I finished analysis on these SNPs already, but there were some job exceed the limited walltime and got killed. I repeat analysis again by splitting methylFreq to be smaller (3000 CpGs/job)
  • Working directory on Triton cluster: /oasis/triton/scratch/nplongthongkum/20130219_UCLA_SZ_mQTL/20130324_mQTL_5M.SNPs_maf0.1_repeat
  • Data summary
    • Total tested CpGs: 76,408 autosomal CpGs
    • Total tested SNPs: 5,257,772 SNPs
    • Total number of tests: 445,603,530
- Copy all output onto genome-miner: /home/nplongth/Noi_scratch/DataAnalysis_2013/UCLA_mQTL_5MSNPs_repeat_20130326
- Printout all output and sort out p-value
cat mQTL-maf0.1.xa* > all.mQTL-maf0.1.output (445603530)
- To reduce sorting time, I print out only minimum p-value 0.05 for multiple testing correction (the result will be the same 
since the number of test will be specified in length parameter (n) in the R script.
awk '{if($7<=0.05) print $0}' all.mQTL-maf0.1.output | sort -g -k7 > sorted.mQTL-maf0.1.minp0.05
BH.correction.R
source("http://bioconductor.org/biocLite.R")
biocLite("multtest")
library(multtest)
A=read.table("sorted.mQTL-maf0.1.minp0.05")
B=A$V7
p_adj=p.adjust(B, method = "BH", n=445603530)
dat <- cbind(A, p_adj)
write.table(dat, file="BH.mQTL-maf0.1.minp0.05.output", sep="\t",append=F, row.names=FALSE, col.names=FALSE)