Dear  Giovanni

Interesting! 

One your #1, perhaps you could have used the approximation to the original base quality scores. How strictly are those normalized data validated?   What is the confidence?

#s 2 and 3,  Is it a kind of Phred or BAQ calculation for your samples?

Could you try joint association test for calculating SNPs?  There are some logistic kernel machine tests which can be done for checking the point 2  that you raised. 

There is a workflow in Galaxy which you could use for all the steps which you suggested. 

Best

Prash

Dear Prash,

I have not cheked strickly the distribution fitness, since the analysis is simple and no factorial interactions have to be considered. However, I realized that Chi2 test under normal distribution better found significant difference (less number of significant SNPs, i.e. the most important ones) compared to Chi2 test under binomial distribution (many significant SNPs, even at very small difference in proportipons). In this case, I'm wondering if the choice of data distribution can be based on concept rather than on mathematical models.

Quality score is right a Phred score.

Association tests such as Fisher exact test should be fine, but I found it similar to Chi2 test under binomial distribution, which was not  fine for me (see above). Chi2 test under binomial distribution is right a logistic regression (I believe).

I set up a workflow in Galaxy to obtain a Fastq Tabular from raw Fastq for the statistical analysis in SAS.

Regards,

Giovanni

I see,  nice  Giovanni :-)

Choice of data could certainly be based on models.  I think you could even train them on a network of models.  if you are looking for stringent confidence. 

I suggest you please await for more answers from experts here. 

PCR amplicons from mRNA

Alex,

of course, same conditions for the samples.

As in many NGS experiments, no replicates are used because of high costs. However, some bioinformatics tools are implemented with Z-test for such analysis. As you said, here we can say only that the difference are not due to chance. Nevertheless, significance detection in binomial distribution is very different from normal distribution, so where is the truth?

Therefore, if we go back to the original questions:

1. Can I use the approximation to normal distribution and generalized linear model (proc GENMOD in SAS)?

2. Can I sue the quality score Q (Phred score) as weight for the counts?

3. And now, is Z-test a Chi2 test on proportions under binomial distribution? Or, how can Z-test be done in SAS?

Note that any test you can do with this data is based on the reads within these two samples. So the null hypothesis is actually that the reads of AAT and AAC are taken from samples with the same ratio (AAT:AAC). This is not a very sensible hypothesis - you compare two individuals / sequencing runs, and surely they will differ. So this H0 is false on a priori grounds. No need to test.

Further, the size of the difference is crystal clear from the data. I wonder what more you need to interpret this data.

Finally, if you want a p-value, the normal approximation is ok. Anything is ok, because there are so many reads that any kind of test (even the least powerful and most suboptimal) will give you a p-value very close to zero. I'd simply use the Chi² test, what gives X²=548135 on 1 d.f. and p

Thanks Jochen for your quick and valuable answer.

I understand (and I know) that biological and technical replicates are required for such comparisons, but as you know, high costs limit this approach. Therefore, this is just a preliminary work.

In the mentioned situation, the significance is clear, even withouth a statistics. However, I tested all the 64 codons occurring in the samples. I found that under binomial distribution, even very little differences in proportions yielded significance, but normal distribution was more conservative and yielded, perhaps, more realistic significance . For example, I found all the codons other than AAT with

Hi Giovanni, I didn't follow all the discussion an I hope but suggestions will not be redundant. In NGS sequencing data the probability of a read (or a codon) to be sequenced usually follows a negative binomial distribution instead of a Poisson. You can try to use a G-test instead of a Chisq. Your analysis seems more similar to a bulk segregant analysis (without using bulks though). Try to give a look at the link attached, maybe you could adapted their approach to your case.

Hope this help

http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3207950/

Dear Andrea,

many thanks for the suggestion. Perhaps, G-test could be appropriate in my case, although I cannot understand all the explanations reported in the article, as I'm not a statistician. However, I understand that G-test is an alternative to Chi2-test and Fisher's exact test.

Can you tell me how to do a G-test in SAS?

I'm not a statistician either unfortunately. I am sure the test are practically the same, but with G-test you have one peer-reviewed reference that could support you approach (like bulk segregant analysis without bulks). I am not confident with SAS unfortunately, but in R there is the likelihood.test (http://www.rforge.net/doc/packages/Deducer/likelihood.test.html). If you need help with R send me a private message here Research Gate and we can arrange something.

I found that tests under binomial distribution are very different from those under normal distribution.

I agree that biological (and perhaps technical) replicates are required for a correct comparison of means.

This is a preliminary work, and in the future more replicates can be included in the same run, although, as you mentioned, the costs for library preparation multiply.

I think you could do a Odds Ratio test

Total 4,343,588

Controls, D-

        E-: AAT 4,308,215 (99,186%)

        E+: AAC 708 (0,092%)

Cases, D+

        E-: AAT 1,397,138 (86,836%)

        E+: AAC 197,248 (12,260%)

Consider AAC as a disease causing exposure, and  AAT as harmless non-exposure.

                E+                       E-

D+    a=197,248       b=4,308,215

D-     c=708               d=1,397,138

For your example SNP, I use this web-page to calculate OR:

http://www.medcalc.org/calc/odds_ratio.php

just plug in the a, b ,c and d, then test.

Odds ratio 90.3486

95 % CI 83.9194 to 97.2703

z statistic 119.581

P < 0.0001

The logged OR follow Normal, so they got a z score and p-value for you, the web-calculator rounded it up, but I think its very very low, because you have a decent sample size.

You might want to control sex, age, race group and do a stratified OR test, usually a even lower p-value will turn up. The web-page calculator cannot do this though. 

I bet you are doing a whole genome screening, you could use Plink1.9 to run stratified odds ratio test for billions of variants, but you will have to study the plink data format and some of its commands.

The other problem hanging is multiple testing, I don't know how many genome variants(mostly SNPs I guess) you are about to test. I hope some of them can at least reach the GWAS 10-8 threshold, good luck.

Dear Xiaoran,

odds ratio seems interesting. Might it be a Z-test? Thank you!

Giovanni

Dear Giovanni,

The odds ratio(OR) itself is not valid for z-test, but the logged odds ratio, Log(OR) is asymptotically normal distributed given large sample size, and the variance of Log(OR), Var(Log(OR)) is 1/a+1/b+1/c+1/d. So yes, the odds ratio test is actually a z-test, the z-score Z=Log(OR)/[Var(Log(OR))]1/2 follows standardized normal distribution, p-value is calculated from the z-score. The logged confidence interval(CI) is also calculated.

When they report back, most package will hide the logged OR but report the original OR, and the logged CI is also "nature exponent-ed" back to usual CI for you.

Dear Xiaoran,

can you tell me how to compute Z-test on ln(OR) in SAS?

SAS can definitely do that for ONE SNP after another, with one procedure, I'll  look that up for you.

The real problem I can think off is your data format and the amount of SNP you are to test, if you are doing the whole genome screening,

I'm not sure how to do a "For Loop" in SAS.

ok, please let me know the SAS procedure doing Z-test (I cannot find it). I will try to manage for the rest. Thank you in advance.

I couldn't find a SAS procedure which just let you put in the four numbers and get a OR test for you.

SAS's FREQ procedure can do the Odds Ratio test. The problem is that you couldn't just plug in the four numbers, but instead you have to make a table of all your samples, listing the exposure(bad SNP)/non-exposure(normal SNP) status as 1/0 in one of the table column, and the disease/non-disease in another column, and than ask SAS to do the OR test.

Here is the link to SAS manual for OR test, under FREQ procedure.

http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_freq_a0000000644.htm

Can you do some R coding? If you have millions SNPs test, SAS is not a good tool.

Send me a little piece of your sample, disease status, a few SNPs, I could work the R code. You could use my email: [email protected]

Yes, I know FREQ procedure. Thank you. I'm not a R user, unfortunately. However, SAS can work with millions of reads as well. Thanks again Xiaoran.

Giovanni

You are welcome, it's good for me to know that SAS can do that much as well, thanks for the information.