Noi/NOTES/2012-4-12

From ZhangLabWiki
Jump to navigation Jump to search
  • Link to calendar: [[1]]
  • To do list today
    • Have finished calculating FDR by BH correction and re-plot manhattan plots again. The results looked the same just change the p-value cutoff at FDR 10%.
    • continue mpo analysis --> convert tstat value to p-value . Then perform multiple testing correction by multtest (R) to get q-value (BH FDR)

ASM analysis of UCLA SZ data set[edit]

  • Continued from: http://genome-tech.ucsd.edu/LabNotes/index.php/Noi/NOTES/2012-3-24
  • I have done multiple testing correction to calculate FDR using multtest function in R
  • Working directory: genome-miner,
  • /home/nplongth/Noi_scratch/ASM_UCLA/combined_ASM_UCLA-4batches_2012_01_24
  • I did multiple testing correction by printing out the SNP-CpG ASM and p-value then use multtest R function to calculate FDR by BH correction.
  • To print out the SNP-CpG ASM position and the corresponding p-value (min p-value 0.05) from the T-test
awk '{print $1"\t"$98}' UCLA_SZ_SeqDepTest_TTest_PVal0.05.txt | sort -g  -k2 > sorted_SNP-CpG_pval_UCLA_SZ_SeqDepTest_TTest_minpval0.05.txt

9124 sorted_SNP-CpG_pval_UCLA_SZ_SeqDepTest_TTest_minpval0.05.txt

  • Count total number of test by editing script to allow all p-value to be counted:
./asmTTest.pl UCLA_MFASM_Qcpg_list.txt > count_total_numberoftests --> Total 73,271 tests

Multiple testing correction[edit]

multtest_R.R

source("http://bioconductor.org/biocLite.R")
biocLite("multtest")
library(multtest)
A=read.table("sorted_SNP-CpG_pval_UCLA_SZ_SeqDepTest_TTest_minpval0.05.txt",header=TRUE)
B=A$p_ttest
p_adj=p.adjust(B, method = "BH", n = 73,271)
summary(p_adj)
write.table(file="UCLA_SZ_SeqDepTest_TTest_pval2qval_p.adjust.BH.txt", p_adj, append = FALSE,row.names = FALSE)
  • I decided to selected the SNP-CpG ASM with FDR < 20% as Dr. Zhang said that it is still acceptable even 5% FDR is ideal. Total I got 8,646 candidate ASM.
  • Sort out the original ASM_TTest matrix based on the p-value
sort -g  -k98 UCLA_SZ_SeqDepTest_TTest_PVal0.05.txt > sorted_UCLA_SZ_SeqDepTest_TTest_PVal0.05.txt
  • Print out the ASM that have FDR <= 10% and 20%
head -8336 sorted_UCLA_SZ_SeqDepTest_TTest_PVal0.05.txt > 10%FDR_sorted_UCLA_SZ_SeqDepTest_TTest.txt (8,335 sites, min p-value ~ 0.01142...)
head -8647 sorted_UCLA_SZ_SeqDepTest_TTest_PVal0.05.txt > 20%FDR_sorted_UCLA_SZ_SeqDepTest_TTest.txt (8,646, min p-value ~ 0.02358...)
  • Extract number of significant ASM
    • SNP-free ASM: 665 sites (8%)
    • SNP on CpG ASM: 7,670 sites (92%)
./getnoSNPCpG_edit.pl 10%FDR_sorted_UCLA_SZ_SeqDepTest_TTest.txt > 10%FDR_noSNPonCpG_UCLA_SeqDep_TTest.txt
./getSNPonCpG.pl 10%FDR_sorted_UCLA_SZ_SeqDepTest_TTest.txt > 10%FDR_SNPonCpG_UCLA_SeqDep_TTest.txt