search

howrigan / trio_sequence_analysis

3 stars 4 forks source link

The issue of DNM_Burden script. #1

Open LIKUOKUO opened 4 years ago

LIKUOKUO commented 4 years ago

Dear Dr. Dan Howrigan. I read your excellent paper and got questions about the DNM burden, and hope your can help me. I was not familiar with poisson.test(). Why you used *conf.intgrp2_rate as dnm_lowci and dnm_highci in case-control enrichment and mod$estimate/p_all** as dnm_enrichment in Model enrichment. In addition, how to get 95% confidence intervals of DNMs rate in case and control as figure 1 in you paper. Hope you can understand what I mean.

Here are the code in your script: case-control mod <- poisson.test(c(sum(grp1$all),sum(grp2$all)),c(sum(grp1$TRIOS),sum(grp2$TRIOS))) dnm[i] <- sum(grp1$all) dnm_rate[i] <- sum(grp1$all) / sum(grp1$TRIOS) grp2_rate <- sum(grp2$all) / sum(grp2$TRIOS) dnm_lowci[i] <- as.numeric(mod$conf.int[1]grp2_rate) dnm_highci[i] <- as.numeric(mod$conf.int[2]grp2_rate) dnm_pval[i] <- mod$p.value dnm_enrich[i] <- mod$estimate

Model enrichment mod <- poisson.test(sum(grp1$all),sum(grp1$TRIOS),p_all) dnm[i] <- sum(grp1$all) dnm_rate[i] <- sum(grp1$all) / sum(grp1$TRIOS) dnm_lowci[i] <- as.numeric(mod$conf.int[1]) dnm_highci[i] <- as.numeric(mod$conf.int[2]) dnm_pval[i] <- mod$p.value dnm_enrich[i] <- mod$estimate / p_all

howrigan commented 4 years ago

Hi - thanks for reaching out!

For your first question, the confidence intervals returned from the 2-sample poisson.test() are the confidence intervals around the rate ratio of the 1st group vs the 2nd group, and multiplying them by the rate in the 2nd group gives the confidence interval in terms of absolute rates. For example, if the rate in the 1st group is 1.2 and the 2nd group is 0.8 (rate -ratio=1.5), and the confidence intervals returned are 1.3 (lower CI) and 1.7 (upper CI), than the conversion I use says the confidence intervals around 1.2 are 1.04 (1.30.8) and 1.36 (1.70.8). For the 1-sample model comparison, comparing an observed rate to a fixed rate doesn't return the rate ratio, but just the absolute rate, so I divide that by the fixed rate to get the rate ratio.

To get confidence intervals for both samples in figure 1, I just switch group 1 and 2.