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

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

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

图片.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)   #输出主要结果

结果文件如下:

image-20200928151036170

3. Gfold软件

image-20200928143759190

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

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

image-20200928144553377

安装GSL

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

image-20200928145510085
# 安装最新版的
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
image-20200928152732420
  1. 接着输入make, 发现报错了。
image-20200928154823119
  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
    
image-20200928155428074

差异分析

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

    gfold diff -s1 Control.txt -s2 Treat.txt -o ControlVSTreat.csv
    
    image-20200928162344142
  2. 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。

    image-20200928162708344

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

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 195,980评论 5 462
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 82,422评论 2 373
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 143,130评论 0 325
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 52,553评论 1 267
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 61,408评论 5 358
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 46,326评论 1 273
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 36,720评论 3 386
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 35,373评论 0 254
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 39,678评论 1 294
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,722评论 2 312
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,486评论 1 326
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,335评论 3 313
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 37,738评论 3 299
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,009评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,283评论 1 251
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 41,692评论 2 342
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 40,893评论 2 335