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

简介: 设置生物学重复这个环节也是你实验设计很重要的一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样品浓度决定的覆盖度参数,之二是真实基因表达量后验分布的比较为寻找差异表达基因提供了一个参照。这种方法针对无重复样本的数据是有一定优势的,这里提供一个链接大家有兴趣的话可以去看该博主的讲解,之后有机会也会尝试一下该软件的使用进行比较。


相关文章
|
搜索推荐 Linux Python
VET:一个基于R语言的VCF数据提取工具,支持按基因ID、物理位置、样品名称提取指定变异信息
VET:一个基于R语言的VCF数据提取工具,支持按基因ID、物理位置、样品名称提取指定变异信息
|
6月前
|
数据可视化
R语言生态学进化树推断物种分化历史:分类单元数与时间关系、支系图可视化
R语言生态学进化树推断物种分化历史:分类单元数与时间关系、支系图可视化
R语言生态学进化树推断物种分化历史:分类单元数与时间关系、支系图可视化
|
6月前
|
机器学习/深度学习 人工智能 运维
人工智能平台PAI 操作报错合集之请问Alink的算法中的序列异常检测组件,是对数据进行分组后分别在每个组中执行异常检测,而不是将数据看作时序数据进行异常检测吧
阿里云人工智能平台PAI (Platform for Artificial Intelligence) 是阿里云推出的一套全面、易用的机器学习和深度学习平台,旨在帮助企业、开发者和数据科学家快速构建、训练、部署和管理人工智能模型。在使用阿里云人工智能平台PAI进行操作时,可能会遇到各种类型的错误。以下列举了一些常见的报错情况及其可能的原因和解决方法。
|
6月前
|
机器学习/深度学习 数据可视化
数据分享|R语言生存分析模型因果分析:非参数估计、IP加权风险模型、结构嵌套加速失效(AFT)模型分析流行病学随访研究数据
数据分享|R语言生存分析模型因果分析:非参数估计、IP加权风险模型、结构嵌套加速失效(AFT)模型分析流行病学随访研究数据
|
6月前
R语言估计多元标记的潜过程混合效应模型(lcmm)分析心理测试的认知过程
R语言估计多元标记的潜过程混合效应模型(lcmm)分析心理测试的认知过程
|
数据挖掘 Go
文献丨群体转录组分析锁定关键转录因子
文献丨群体转录组分析锁定关键转录因子
|
算法 Go
差异分析|DESeq2完成配对样本的差异分析
差异分析|DESeq2完成配对样本的差异分析
421 0
差异分析|DESeq2完成配对样本的差异分析
|
JSON 算法 数据格式
优化cv2.findContours()函数提取的目标边界点,使语义分割进行远监督辅助标注
可以看到cv2.findContours()函数可以将目标的所有边界点都进行导出来,但是他的点存在一个问题,太过密集,如果我们想将语义分割的结果重新导出成labelme格式的json文件进行修正时,这就会存在点太密集没有办法进行修改,这里展示一个示例:没有对导出的结果进行修正,在labelme中的效果图。
212 0
|
数据可视化 Serverless Go
scRNA分析|单细胞GSVA + limma差异分析-celltype分组?样本分组?
scRNA分析|单细胞GSVA + limma差异分析-celltype分组?样本分组?
996 0
|
Shell Linux 测试技术
CORNAS:一种快速简单鉴定无重复转录组差异基因的方法
还记得上次文章的最后提到CORNAS这种方法吗?最近刚好在Github上看到了这个项目,就花了点时间看了下文档感觉操作也比较简单,这里记录一下使用过程,大家共同学习一下。
183 0