没有生物学重复的转录组数据怎么进行差异分析?

简介: 设置生物学重复这个环节也是你实验设计很重要的一part,设置的好对你下游分析也有利,通常我们做转录组测序,需要的样本量每组至少为3个生物学重复,这个处理起来就很合理,并且现在流行的差异分析软件DEseq2,limma,edgeR等等都是针对有重复的数据去做的,但有时候会不幸碰到样品测序失败不能用,导致每组就给你剩一个重复时候该怎么办,之前我有批数据就是这样,但是办法总比困难多不能放过任何实验数据,搜了搜其实还是有一些方法可以去解决的,在这里介绍下我搜到的几种方法。

假如现在你手头有如下文件(test.txt),只有俩样品RPKM_A (对照) 和RPKM_B (处理), 值为标准化后的RPKM。

d3fd0a8ca7bd1820b3d67d9f2994ee3.png

1. 根据foldchange直接筛选

之前在一篇中文文献中见到有人用这种方法,作者自定义差异基因的标准:至少有一组RPKM值大于5,且满足foldchange(差异倍数)  > 2,我们可以在LInux中直接可以用awk进行过滤,其实Excel、R中也可以操作,根据个人习惯吧。代码如下:


### 上调基因########
# 提取B组大于等于5,A组等于0的基因。
less test.txt | gawk  '{if (($2==0)&&($3>=5)) print $0}'  > up.txt
# 提取A、B俩组至少有一组大于等于5,且B组值/A组值大于等于2
 less test.txt | gawk  '{if (($2!=0)&&($3!=0)) print $0}'|gawk '{if (($2>=5)||($3>=5)) print $0}'|sed '1d'|gawk '{if ($3/$2>=2) print $0}' >> up.txt
### 下调基因#########
# 提取A组大于等于55,B组等于0的基因
 less test.txt | gawk  '{if (($2>=5)&&($3==0)) print $0}'  > down.txt
# 提取A、B俩组至少有一组大于等于5,且A组值/B组值大于等于2
 less test.txt | gawk  '{if (($2!=0)&&($3!=0)) print $0}'|gawk '{if (($2>=5)||($3>=5)) print $0}'|sed '1d'|gawk '{if ($2/$3>=2) print $0}'  >> down.txt

2. edgeR包

这种方法我在提到过,edgeR包可以做无重复的差异分析,不过需要认为指定一个dispersion值(设置BCV值),这样得到的结果比较主观,不同的人就可以有不同的结果。通常如果是实验控制的好的人类数据,那么选择BCV=0.4,比较好的模式生物选择BCV=0.1。参考

代码如下:


library(edgeR)
##跟DESeq2一样,导入数据,预处理(用了cpm函数)
exprSet<- read.table(file = "test.txt", sep = "\t", header = TRUE, row.names = 1, stringsAsFactors = FALSE)
group_list <- factor(c(rep("Contral",1),rep("Treat",1)))
##设置分组信息,并做TMM标准化
exprSet <- DGEList(counts = exprSet, group = group_list)
bcv = 0.1  #设置BCV值
et <- exactTest(exprSet, dispersion=bcv^2)
write.csv(topTags(et, n = nrow(exprSet$counts)), 'result.csv', quote = FALSE)   #输出主要结果

结果文件如下:

44c301e3b821e106e616c1458833b12.png

3. Gfold软件

93eb33b21eee973664308446754c1bb.png

地址:https://zhanglab.tongji.edu.cn/softwares/GFOLD/index.html

Gfold软件应该是做没有生物学重复样本首选的软件,该软件由同济大学开发的,网站 往下拉可以看到该软件的几个功能,其中Example3为鉴定无重复的数据的差异基因。另外,这个软件不支持Windows 版本,是基于linux的一个安装软件,所以我们需要在linux环境下使用

9c54ea802a5bbcdd97b8302988e7536.png

安装GSL

使用Gfold之前必须先安装Gsl,如图下载最新的版本

f9f241cd1e7ce6e019516442b412632.png


# 安装最新版的
wget http://mirrors.ocf.berkeley.edu/gnu/gsl/gsl-latest.tar.gz
tar -zxv -f gsl-latest.tar.gz
.configure --prefix=/home/pub_guest/hekai/soft_ware/GFOLD/gsl-2.6/
make
make check(选做)
make install

安装Gfold

  1. 下载安装包


wget https://zhanglab.tongji.edu.cn/softwares/GFOLD/gfold.V1.1.4.tar.gz
tar -zxv - f gfold.V1.1.4.tar.gz
cd gfold.V1.1.4
  1. 我们打开目录下的README文件,可以需要执行下面俩句。注意 /your/installed/path/ 是你安装GSL的路径,把下面命令行中的替换为相应的路径即可


# 打开你的bashrc文件,再最后添加
export CXXFLAGS="-g -O3 -I/your/installed/path/include -L/your/installed/path/lib" 
export LD_LIBRARY_PATH="/your/installed/path/lib:"$LD_LIBRARY_PATH 
source ~/.bashrc

48b08c0c4c386532d59c53e9f788556.png

  1. 接着输入make, 发现报错了。

95d668f227af299673096ca46e97682.png

  1. 不慌,文档里也有提醒如果报错 直接输入以下命令行即可


# If it happens, follow step 1 again. If error remains, try the following command:
g++ -O3 -Wall -g main.cc -o gfold -lgsl -lgslcblas -I/your/installed/path/include -L/your/installed/path/lib
  1. 这是我们就会发现目录下多了一个gfold软件,是可执行状态,我们.gfold -h ,可以看到该软件的帮助文档,证明此时已经安装成功了,我们将其添加为环境变量里方便我们使用。
echo 'export PATH=/your/installed/path:$PATH' >>~/.bashrc
######  **/your/installed/path/** 是你安装**gfold.V1.1.4**的路径
source ~/.bashrc

55835e2900f32114b1d9a194fc61c3b.png

差异分析

  1. 我们需要准备俩个文件Control.txt和Treat.txt,我们看下处理组Control.txt文件都包含哪些,Gfold输入文件规定必须为5列,前俩列可以分别输入你的基因ID号和Symbol号,俩列内容一样其实也不影响,第三列为原始Counts值,第四列为基因长度,最后一列为标准化的RPKM值,同样对照组Treat.txt文件也按照这样准备。
  2. 准备完毕后,直接一条命令计算差异。
gfold diff -s1 Control.txt -s2 Treat.txt -o ControlVSTreat.csv

64d9c2e3bddcaf2db2cf175acb43d9d.png

  1. OK,已经计算完了,我们看下结果文件都有哪些内容。主要一共7列信息,前两列没什么可说,就是gene symbol和gene  name,第三列是GFOLD值,相当于log2(Fold  Change),该值等于0的基因则记为非差异基因,非0的值才是差异基因,E-FDR是基于重复的Empirical  FDR,因此无重复样本的经验FDR均为1。Log2fdc以及后面的RPKM列可以忽略考虑,因为最开始的exon的长度,我们是给定的是一个虚拟的数据。所以真正的确定差异是否显著,主要是看GFOLD值,GFOLD>0,代表case组中高表达,GFOLD<0,代表case组中低表达。后续筛选差异基因如果太多或者太少,我们也可以通过设定GFOLD的阈值来控制,比如设定<-0.3或者大于0.3。

e88f967f80fcadfa1cae158e7a83789.png

除上述之外几种方法,还有几个比较经典的软件可供选择进行无重复数据的差异分析,比如BMC Bioinformatics发表过的一篇文章中提出了一种新的差异基因分析方法利用贝叶斯方法来推断真实基因表达数的  后验分布。其创新型之一该方法包括了由RNA样品浓度决定的覆盖度参数,之二是真实基因表达量后验分布的比较为寻找差异表达基因提供了一个参照。这种方法针对无重复样本的数据是有一定优势的,这里提供一个链接大家有兴趣的话可以去看该博主的讲解,之后有机会也会尝试一下该软件的使用进行比较。


相关文章
|
4天前
|
机器学习/深度学习 数据可视化
数据分享|R语言生存分析模型因果分析:非参数估计、IP加权风险模型、结构嵌套加速失效(AFT)模型分析流行病学随访研究数据
数据分享|R语言生存分析模型因果分析:非参数估计、IP加权风险模型、结构嵌套加速失效(AFT)模型分析流行病学随访研究数据
|
10天前
R语言估计多元标记的潜过程混合效应模型(lcmm)分析心理测试的认知过程
R语言估计多元标记的潜过程混合效应模型(lcmm)分析心理测试的认知过程
31 0
|
2月前
【SPSS】两独立样本的极端反应检验和两配对样本的非参数检验详细操作教程(附案例实战)
【SPSS】两独立样本的极端反应检验和两配对样本的非参数检验详细操作教程(附案例实战)
47 0
|
8月前
|
数据挖掘 Go
文献丨群体转录组分析锁定关键转录因子
文献丨群体转录组分析锁定关键转录因子
|
8月前
|
数据可视化 数据挖掘 Linux
转录组下游分析丨利用limma包进行差异表达分析,结果可视化绘制火山图和热图
转录组下游分析丨利用limma包进行差异表达分析,结果可视化绘制火山图和热图
|
10月前
|
算法 Go
差异分析|DESeq2完成配对样本的差异分析
差异分析|DESeq2完成配对样本的差异分析
301 0
差异分析|DESeq2完成配对样本的差异分析
|
10月前
|
数据可视化 Serverless Go
scRNA分析|单细胞GSVA + limma差异分析-celltype分组?样本分组?
scRNA分析|单细胞GSVA + limma差异分析-celltype分组?样本分组?
709 0
|
10月前
|
Shell Linux 测试技术
CORNAS:一种快速简单鉴定无重复转录组差异基因的方法
还记得上次文章的最后提到CORNAS这种方法吗?最近刚好在Github上看到了这个项目,就花了点时间看了下文档感觉操作也比较简单,这里记录一下使用过程,大家共同学习一下。
106 0
|
10月前
|
数据挖掘 Serverless
Robust火山图:一种含离群值的代谢组数据差异分析方法
代谢组学中差异代谢物的识别仍然是一个巨大的挑战,并在代谢组学数据分析中发挥着突出的作用。由于分析、实验和生物的模糊性,代谢组学数据集经常包含异常值,但目前可用的差异代谢物识别技术对异常值很敏感。作者这里提出了一种基于权重的具有稳健性火山图方法,助于从含有离群值的代谢组数据中更加准确鉴定差异代谢物。
114 0
|
10月前
|
存储 移动开发 Python
从Scanpy的Anndata对象提取信息(适用于单细胞转录组且涉及h5文件读写)
从Scanpy的Anndata对象提取信息(适用于单细胞转录组且涉及h5文件读写)