R语言read.table的问题

Python010

R语言read.table的问题,第1张

建议你明确地设定 header 参数。按照惯例,首行只有对应列的字段而没有行标签对应的字段。因此,它会比余下的行少一个字段。

(如果需要在 R 里面看到这一行,设置 header = TRUE。)如果要读取的文件里面有行标签的头字段(可能是空的),以下面的方式读取!

今天就先来聊聊如何看差异表达基因数据,火山图,聚类图又怎么看。1差异基因筛选方法那差异基因是如何筛选出来的呢?差异基因的筛选方法有很多,包括倍数法、T检验、F检验及SAM等。

下面简单介绍一下GCBI上用的倍数法和SAM法。

倍数法适用于没有生物学重复的样本,其计算基因在两个条件下表达水平的比值,确定比值的阈值,将绝对值大于此阈值的基因判断为差异基因。

SAM算法适用于有生物学重复的样本,通过对分母增加一个常量 T 检验过程减小了假阳性发生的概率。文献中报道,相较于其他算法,SAM算法更为稳定,筛选出的结果也更为准确。2差异基因数据解读经过合适的差异基因方法筛选出的差异基因,结果一般分为两部分,数据+图形。

数据结果展示如下图所示(两分组)众多参数中,重点看三个。p-value或q-value没有做生物学重复请跳过这一步。

p-value或q-value是统计学检验变量,代表差异显著性,一般p-value或q-value小于0.05代表具有显著性差异,但可根据具体情况适当调整。

因为p-value或q-value衡量地是某个基因假阳性的概率,如果p-value或q-value越低,那么挑选该基因出现假阳性的概率就越低,可验证性就越高。

两者具体的计算方法具体如下:那p-value、q-value同时存在时看哪个呢?

SAM法只有q-value。当两者同时存在时,可根据具体情况具体分析。

差异筛选是一个典型的多重假设检验过程,对于多重假设检验,单次检验中差异显著基因的假阳性率(p-value较小)可能会较大,而q-value和FDR值较常见的BH校正方法得到的FDR值而言,改进了其对假阳性估计的保守性。

即q-value相比于p-value更加严格,当差异基因结果较少时,可以退而求其次看p-value。Fold ChangeFold Change表示实验组比上对照组的差异表达倍数,一般表达相差2倍以上是有意义的,放宽要求1.5倍或者1.2倍也可以接受。

看表达倍数的同时还需结合基因表达丰度,信号值太低的基因会在后续的验证实验中检测不到。3差异基因图表解读在差异结果的图形展示结果中,主要是火山图和聚类图。火山图火山图只针对两分组且有生物学重复的情况。

如何看火山图呢?火山图可反映总体基因的表达情况,横坐标代表log2(Fold Change),纵坐标表示-log10(P值),每个点代表一个基因,颜色用以区分基因是否差异表达,图中橙色的点代表差异表达基因,蓝色的点代表没有差异表达的基因。聚类图聚类图可以衡量样本或基因之间表达的相似性。

如上图所示的聚类图中,横坐标代表样本聚类,一列代表一个样本,聚类基于样本间基因表达的相似性,样本间基因表达越接近,靠的越近,以此类推。

纵坐标代表基因聚类,一行代表一个基因,聚类基于基因在样本中表达的相似性,基因在样本中表达越接近,靠的越近,以此类推。

色阶代表基因表达丰度,越红代表上调得越明显,越绿代表下调得越明显。

如何做聚类图请戳往期推送做个聚类图只需1分钟

差异基因有了,如何挑选潜在基因进行实验验证呢?

关键还在于感兴趣点在哪了。粗略的看,可以先看KEGG或者GO功能分类,看差异基因具体富集在哪些通路或功能。

比如关注的是细胞内酸合成关键酶,可以重点看酸合成和碳流相关通路。具体如何看KEGG或者GO功能分类,请听下回分解。

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数据上传 、 二代测序数据解读 。