Noi/NOTES/2012-4-12
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)
- Continued from :http://genome-tech.ucsd.edu/LabNotes/index.php/Noi/NOTES/2012-3-24
- continue ASM analysis --> find FDR by printing out SNP-CpG ASM and p-value then perform multiple testing correction. Also need to find total number of tests.
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