1. 定义
差异分析就是分析两组数据是否有差异,请注意分析差异的原则是建立在组与组之间.它的目的是比较两个条件(包括种属、表型等)下的基因表达差异,通过一定的统计学方法，从中识别出与条件相关的特异性基因，然后进一步分析这些特异性基因的生物学意义.
基因表达差异分析的第一步是要识别在两个条件下有显著性表达差异的基因,简称差异表达基因.那么怎样才能称得上显著性表达差异?通常的做法是对两组数据的差异倍数进行统计学检验,得到的P value达到某个阈值,则为显著差异.
2. preparation
2.1 输入
差异表达的输入是reads counts矩阵,请注意,这种矩阵是没经过标准化的.即上一节的genes_counts_matrix.
除此之外,输入还需要样本信息表.样本信息表需要描述哪些样本属于同一个分组.这个需要我们自己构建.我们之前在下载测序数据的时候,还手动构建了一个sra_info.tsv文件,可以根据那个修改.
awk '{print $2}' sra_info.tsv | cut -b 1,2,3 > sample_info.txt
awk '{print $2}' sra_info.tsv > temp.txt
#paste 用于合并列
paste sample_info.txt temp.txt >> temp2.txt
rm sample_info.txt
rm temp.txt
mv temp2.txt sample_info.txt
准备上述文件后,我们需要构建constrast表,表示哪两个组进行差异表达分析.
CS2 CS1
CS3 CS1
2.2 软件包
我们一般用DESeq2/edgeR做差异表达分析,这两个包是经过大量文献验证,最为可靠的两个包,不过这两个软件都是用R语言写的.
到目前为止的流程还没学习R语言,因此我们用trinity软件包中,一个run_DE_analysis.pl来执行差异表达分析.
虽然trinity主要用于无参转录组的拼接,但也可以用于差异表达分析.
使用mamba install trinity安装包,我第一次安装的时候没注意,mamba要安装2.3版本的,实在有点老了,而且还出现了大量包降级,最后还安装失败了.这边建议新建一个环境安装..
安装完成后,我们可以在~/miniconda3/envs/trinity/bin/run_DE_analysis.pl运行软件.
3. 差异表达分析流程
我们也可以看下LinuxLinuxLinux内置的文档,输入矩阵不需要标准化的原因是不论是DESeq2/edgeR还是当前这个软件,都会用自己内置的算法标准化.下面解释下参数
-
--matrix|-m : 不需要标准化的counts矩阵
-
--method : 指的是用什么方法进行差异表达分析.
-
--samples_file|s : 样本信息表
-
--min_reps_min_cpm: 这个参数和
edgeR有关,主要用于去除低表达量的基因.min_reps设置重复次数,cpm是edgeR引入的值.老师说这个参数默认意思是基因至少要在两个重复中,表达量大于1,否则就会被过滤.一般用默认.- 关于为什么要过滤,这个博客写了,还标明了原论文Click.
-
--constrasts : 对比组文件.前面是实验组,后面是对照组.
-
--reference_sample: 指定对照组.这里再具体解释一下就是,一般做实验,我们有多个实验组和一个对照组.这里就是根据前面的样本信息表获取有哪些组,再根据这个参数获取对照组,其他组都视为实验组.这时不需要
--constrast参数.
如果我们没有设置--constrasts和--reference_sample,则默认为每个组两两进行比较.
perl ~/miniconda3/envs/trinity/bin/run_DE_analysis.pl --matrix ../23.merge/genes.counts.matrix --method DESeq2 --samples_file sample_info.txt --contrasts contrast.txt
分析完成后,我们会得到一个文件夹,文件夹上的数字实际是指进程编号,没必要理会.文件夹内的内容如下所示.
- 诸如
genes.counts.matrix.Cs2_vs_Cs1.DESeq2.count_matrix的文件,与原来的count矩阵相比,是过滤了表达量过低的基因. genes.counts.matrix.Cs2_vs_Cs1.DESeq2.DE_results.MA_n_Volcano.pdf是火山图.genes.counts.matrix.Cs2_vs_Cs1.DESeq2.DE_results是差异表达分析的结果.genes.counts.matrix.Cs2_vs_Cs1.DESeq2.Rscript是差异分析的代码.
我们主要是看结果文件.格式化得到内容如下:
- baseMeanA: 该基因在样本A中表达量的平均值.
- baseMeanB: 该基因在样本B中表达量的平均值.
- baseMean: 该基因在所有样本中表达量的平均值.
- log2foldchange: 基因间的差异倍数,这里取了log2log2log2,这里计算是log2BaseABaseBlog_2frac{BaseA}{BaseB}log2BaseBBaseA
- LFCSE: 是对于log2Foldchange估计的标准误差估计,效应大小估计有不确定性.
- stat: 是Wald统计量,它是由log2Foldchange除以标准差所得.
- pvalue: 属于统计学的概率,当pvaluepquad valuepvalue < 0.05 表示有显著差异.
- padj: p.adjust。转录组测序的差异表达分析是对大量的基因表达值进行的独立统计假设检验，存在假阳性问题，因此引入Padj对显著性P值(P.adjust)进行校正.Padj是对P-value的再判断,筛选更为严格.
p-value和padj想后面再研究下,现在先放着
4. 差异表达分析的原理
4.1 question one : 为什么要做生物学重复
首先阐述什么是生物学重复,当我们想测序一个人的基因序列时,如果对这个人测多次,这时技术重复.如果我们再找3个人测同一个基因的序列,这就是生物学重复.也就是说,我们增加了样本数.
我们做差异表达分析是为了找出差异表达基因.如果第一个样本测到的KCNA3KCNA3KCNA3在癌症样本1和正常样本1之间差不多有两倍差异的表达量.但是第二个样本测到的表达量又差不多,这会让我们得不出正确的结论,因此需要多测几个样本减小测序中的一些问题带来的影响.
目前最少能接受生物学重复的个数就是3个.
4.2 question tow : 如何筛选表达基因
一般是从多个样本构造的分布来进行比较.比如上图的样本可以构成一个癌症的分布和一个正常的分布.这个说得比较通俗易懂Click.
4.3 question three : 是否生物学重复越多,鉴定到的差异基因就越多?
正确.如果样本过少,p-value的值可能会偏大,就无法判断是否是差异表达基因.但是如果我们不断地增加样本,差异表达的基因太多了,人为分析是不现实的,因此实际判断往往会增加log2foldchange的条件,让差异表达基因少一点,只研究差异表达明显的基因.
4.4 question four : DESeq2 鉴定到 500 个差异基因,换 edgeR 却有 2000 个,哪个是对的？
这个问题是没有对错的.每个软件的严格程度不同.
4.5 question five : 没有生物学重复可以做差异表达分析吗?
正常来说是不可以的,但是edgeR对此情况也进行了实现.我们要做无生物学重复的分析,只能先告诉edgeR下图参数的值(组内分歧度).如果不知道,按经验可以从0.1开始设起,然后根据差异基因的数量再调整.但是无论结果如何,这种结果可靠性不高.
5. 差异表达分析结果可视化
最终我们获得的是差异基因的列表.有时我们会根据列表绘制火山图.横坐标是log2foldchange,纵坐标是−log10(p−value)-log_{10}(p-value)−log10(p−value),越往上p值越小,因此研究左上和右上的基因.
除了火山图,有时也是有HeatmapHeatmapHeatmap.横行代表基因,纵行代表样本.越红表达量越高.
