数据准备
结果可视化
方法1：interaction.plot()函数来展示双因素方差分析的交互效应。
图形解读：无论哪个品种，v1值均是N2显著高于N1，品种来看，a品种显著高于b品种。
方法2：gplots包中的plotmeans()函数来展示交互效应。
图形解读：每个样点都显示了置信区间，对应横坐标上方显示了样本量。本例中，同一品种下，均是N2显著高于N1
方法3：HH包中的interaction2wt()函数来可视化结果。
图形解读：两个箱图（左下角和右上角）是主效应图，左下角图，单就nitrogen看，N2显著高于N1；右上角图，单就variety看，a品种显著高于b品种；两个线图（左上角和右下角）为交互效应图，左上图，nitrogen应用在不同variety下的影响，即不同品种下，都是N2显著高于N1，同时a品种显著高于b品种；右下图，不同品种在不同nitrogen下的影响，即不同nitrogen下，都是a品种显著高于b品种，同时N2显著高于N1。
参考资料：
System.currentTimeMillis() ：返回当前系统的毫秒数，由于取得的是毫秒数，所以在处理UNIX时间戳的时候需要转换成秒也就是：
long epoch = System.currentTimeMillis()/1000
方法：
1、获取当前系统的UNIX时间戳
System.out.println("获取系统毫秒数方法1："+Long.toString(new Date().getTime()))
System.out.println("获取系统毫秒数方法2："+Long.toString(System.currentTimeMillis()))
注意：以上代码获取的都是系统毫秒数，在实际的操作中我们一般都是记录毫秒说以求记录的精度，当处理UNIX时间戳的时候需要把数据进行处理。
2、将UNIX时间戳转换成系统可以处理的时间
System.out.println(""+new java.text.SimpleDateFormat("yyyy MM-dd HH:mm:ss").format(new java.util.Date (1215782027390L)))
输出：2008 07-11 21:13:47
注意：此时处理的数据为系统毫秒不是UNIX时间戳
3、讲时间转换成UNIX时间戳
long epoch = new java.text.SimpleDateFormat ("dd/MM/yyyy HH:mm:ss").parse("09/22/2008 16:33:00").getTime()
注意：
请注意！对与不同的时区处理上有差异，首先要清楚自己所在的时区。
String timezone_info = System.getProperty("user.timezone")
System.out.println("当前的时区:"+timezone_info)
System.out.println("时区信息:"+TimeZone.getDefault())
输出：
当前的时区:Asia/Shanghai
时区信息:sun.util.calendar.ZoneInfo[id="Asia/Shanghai",offset=28800000,dstSavings=0,useDaylight=false,transitions=19,lastRule=null]
处理不同的时区的方法：
SimpleDateFormat sd = new SimpleDateFormat("yyyy-MM-dd HH:mm:ss")
sd.setTimeZone(TimeZone.getTimeZone("GMT+8"))
String strDate = sd.format(new Date(1215782027390L))
System.out.println("正八区当前时间:"+strDate)
输出：
正八区当前时间:2008-07-11 21:13:47
qRT-PCR是一种相对表达定量的方法，他的计算方法有很多，常用的相对定量数据分析方法是KJ Livak（Applied Biosystems）等人在2001年提出的“比较Ct法相对定量”，即：利用ΔCt值差异来推算基因表达差异（Ct目的基因 – Ct内参基因 = ΔCt），该方法的具体计算方法请参见文章： qRT-PCR相对定量计算详解 。
一般在相对定量的最终结果中，样本间的差异是以表达差异倍数（Fold change）来展现的，如下图：
那么样品间基因表达差异倍数多少则可以认为有差异呢？回答此问题，我们需要明确差异该如何去定义！
如何定义差异：
说道差异大家首先想到的肯定是生物学上的差异，例如同一基因在两个样品间的表达差异倍数，一般这个倍数从1.2、1.5、2倍都是可以的（转录组里面一般是按2倍作为筛选指标，我觉得1.2、1.5也是可以接受的）。
另一方面，我们也应考虑随机误差，因为我们无法消除误差，看上去完美的数据也有可能是随机误差造成的，所以，我们在关注生物学差异之外，还要考虑统计学差异。
以上两种差异都是客观上存在的，我们当然是希望数据差异是由实验处理造成的，但随机误差又是客观存在的，所以随机误差发生的概率越小越好。
如何衡量随机误差？
P值（P-value），想必大家都不会陌生，它是用来判定假设检验结果的一个参数，说直白点就是P值代表了一种可能性，衡量的是随机出错的概率。在统计学中，一般要求P值小于0.05；如果P-value=0.05，意味着我们的实验结果有5%的概率是随机误差引起的。
我们经常用到这样的论述p<0.05（显著），可用一颗星号表示“*”，而两颗星“**”代表p<0.01（极显著）；那是不是p<0.01的数据比p<0.05的好，组间的差异也更大呢？答案是否定的！P值衡量的是随机出错的概率，不能衡量差异量变大小，所以我们不能说一个P值<0.01的结果比P值<0.05的结果具有更大的差异，只能说前者出错的概率更低，或者说组间“差异有统计意义”，而不是组间“具有显著的差异”。
P值的计算：
P值的算法有很多种，最常用的是T检验（T-test），亦称student t检验（Student's t test），主要用于样本含量较小（例如n <30），总体标准差σ未知的正态分布。T检验是用t分布理论来推论差异发生的概率，从而比较两个平均数的差异是否显著。在R语言中T检验用的方法为:t.test(),如果数据不符合正态分布，也就是数据当中有较大的离群值时，可选用非参数秩和检验法，如Wilcoxon test，R语言中对应的方法为：wilcox.test()。关于数据类型及检验方法的选择可参考： 差异统计检验如何选择 。
单样品T检验
例： 某鱼塘水的含氧量多年平均值为4.5mg/L，现在该鱼塘设10点采集水样，问该次抽样的水中含氧量与多年平均值是否有显著差异。
#数据
s<-c(4.33,4.62,3.89,4.14,4.78,4.64,4.52,4.55,4.48,4.26)
shapiro.test(s) #如果P>0.05 符合正态分布
t.test(s,mu=4.5) #T检验, 如果 P>0.05 相等
非配对两样本T检验
例： 为了了解某一降血压药物的效果，将28名高血压病患者随机等分到实验组和对照组，实验组采用新降压药物，对照组则用标准药物治疗，测得治疗前后舒张压的差值如下。问新药和标准药的疗效是否不同？
high<-c(134,146,106,119,124,161,107,83,113,129,97,123)
low<-c(70,118,101,85,107,132,94)
x<-c(high,low)
group<-c(rep("high",12),rep("low",7))
#正态性检验，wilcox.test()
shapiro.test(high) #如果P>0.05 符合正态分布
shapiro.test(low) #如果P>0.05 符合正态分布
#方差齐性检验:如果P>0.05 方差齐
bartlett.test(x~group)
#方法二：car包中leveneTest 检验，spss统计软件默认的检验方法
leveneTest(x~group)
#T检验, 如果 P<0.05 存在差异
t.test(high,low,paired=F,var.equal=T) #如果方差不齐，可更改：var.equal=F,
#或者：
t.test(x~group,paired=F,var.equal=T)
配对两样本T检验
例： 为了解DSCT冠状动脉造影和超声心动图检查两种方法测定心脏病患者左室舒张末容积的差别，某医院收集心脏病患者12例，同时分别用两种检测方法测得其大小如下，问两种检测方法的检测结果是否不同？
ds<-c(82.5,85.2,87.6,89.9,89.4,90.1,87.8,87.0,88.5,92.4)
cs<-c(91.7,94.2,93.3,97.0,96.4,91.5,97.2,96.2,98.5,95.8)
#方差齐性检验，car包中leveneTest
leveneTest(ds,cs)
#作差，正态性检验
#差值正态性检验，差值符合正态分布（P>0.05）
d<-ds-cs
shapiro.test(d)
#配对T检验
t.test(ds,cs,paired=T,alternative="two.sided",conf.level=0.95)
统计检验及绘图-ggpubr
ggpubr包既可以做检验，有可以对统计结果进行整理绘图，输出结果比t检验更加友好。
例： 两种基因型（HH、RR）的水稻品种，分别在高氮和低氮条件下，的测FW、DW和PH三种生理指标数据：
women_weight <- c(38.9,61.2,73.3,21.8,63.4,64.6,48.4,48.8,48.5)
men_weight <- c(67.8,60,63.4,76,89.4,73.3,67.3,61.3,62.4)
mydata <- data.frame(
group= rep(c("Woman","Man"),each=9),
weight = c(women_weight, men_weight)
)
#统计检验
com1 <- compare_means( weight~ group , data = mydata, method = "t.test")
#结果P=0.015，小于0.05，具有显著差异：
#.y. group1 group2 p p.adj p.format p.signif method
# weight Man Woman 0.0154 0.015 0.015 * T-test
绘图显示
install.packages("ggpubr")
library(ggpubr)
p<- ggboxplot(mydata, x="group", y ="weight", color ="group", palette ="jco",add="jitter",short.panel.labs =FALSE)
# 添加p值
p+ stat_compare_means(method ="t.test",label.y=100)
# 显示p值但不显示方法
p+ stat_compare_means(aes(label = ..p.format..),method ="t.test",label.x =1.5)
# 只显示显著性水平
p+ stat_compare_means(aes(label = ..p.signif..),method ="t.test",label.x =1.5)
结果图如下：
http://m.study.163.com/provider/400000000234009/index.htm?share=1&shareId=1031484705
更多生物信息课程：
1. 文章越来越难发？是你没发现新思路，基因家族分析发2-4分文章简单快速，学习链接： 基因家族分析实操课程 、 基因家族文献思路解读
2. 转录组数据理解不深入？图表看不懂？点击链接学习深入解读数据结果文件，学习链接： 转录组（有参）结果解读 ； 转录组（无参）结果解读
3. 转录组数据深入挖掘技能-WGCNA，提升你的文章档次，学习链接： WGCNA-加权基因共表达网络分析
4. 转录组数据怎么挖掘？学习链接： 转录组标准分析后的数据挖掘 、 转录组文献解读
5. 微生物16S/ITS/18S分析原理及结果解读 、 OTU网络图绘制 、 cytoscape与网络图绘制课程
6. 生物信息入门到精通必修基础课，学习链接： linux系统使用 、 perl入门到精通 、 perl语言高级 、 R语言画图
7. 医学相关数据挖掘课程，不用做实验也能发文章，学习链接： TCGA-差异基因分析 、 GEO芯片数据挖掘 、 GSEA富集分析课程 、 TCGA临床数据生存分析 、 TCGA-转录因子分析 、 TCGA-ceRNA调控网络分析
8.其他课程链接： 二代测序转录组数据自主分析 、 NCBI数据上传 、 二代测序数据解读 。
