EBSeq/R/PostFC.R

   1 PostFC=function(EBoutput) {
   2         GeneRealMeanC1=unlist(EBoutput$C1Mean)
   3         GeneRealMeanC2=unlist(EBoutput$C2Mean)
   4         GeneRealMean=(GeneRealMeanC1+GeneRealMeanC2)/2
   5
   6         GeneRealFC=GeneRealMeanC1/GeneRealMeanC2
   7
   8         GeneR=unlist(EBoutput$RList)
   9         GeneR[GeneR<=0 | is.na(GeneR)]=GeneRealMean[GeneR<=0 | is.na(GeneR)]*.99/.01
  10
  11         GeneAlpha=EBoutput[[1]][nrow(EBoutput[[1]]),]
  12         GeneBeta=unlist(sapply(1:length(EBoutput$C1Mean),function(i)rep(EBoutput[[2]][nrow(EBoutput[[1]]),i],length(EBoutput$C1Mean[[i]]))))
  13         GeneBeta=as.vector(GeneBeta)
  14         # Post alpha = alpha + r_C1 * 3
  15         # Post beta = beta + Mean_C1 * 3
  16         # Post Mean of q in C1 P_q_C1= P_a/ (P_a + P_b)
  17         # Post FC = (1-p_q_c1)/p_q_c1 /( (1-p_q_c2)/p_q_c2)
  18
  19         GenePostAlpha=GeneAlpha+3*GeneR
  20         GenePostBetaC1=GeneBeta+3*GeneRealMeanC1
  21         GenePostBetaC2=GeneBeta+3*GeneRealMeanC2
  22         GenePostQC1=GenePostAlpha/(GenePostAlpha+GenePostBetaC1)
  23         GenePostQC2=GenePostAlpha/(GenePostAlpha+GenePostBetaC2)
  24
  25         GenePostFC=((1-GenePostQC1)/(1-GenePostQC2))*(GenePostQC2/GenePostQC1)
  26         Out=list(GenePostFC=GenePostFC, GeneRealFC=GeneRealFC)
  27
  28 }