edgeR对表达量矩阵进行基因表达差异分析
通过前期的基因定量，我们通过**RSEM**或者Kallisto这些工具得到了每个基因或者转录本的表达量矩阵。输出的矩阵包含了基因以及在各个样本中的表达量信息，这些信息将为我们后期判定基因表达量变化趋势以及变化的显著性提供原始参考。
Trinity工具包同样也提供了一系列的脚本用以处理这些基因表达量变化差异显著性分析，在Trinity中主要使用的是edgeR包进行处理。
在此之前我们需要注意一点：对于基因表达量分析由于是要对比不同样本，因此需要每个基因在不同样本中的表达量均有记录（哪怕是无表达量的0也可）。因此，Trinity建议在计算基因表达量时采用的参考转录本库应该是一个库，即先通过所有样本原始数据进行拼接而得到统一的转录本库，进而依照该转录本库进行回帖计算表达量。
在这个工作之前需要两个数据：
1. 基因表达的counts.matrix 文件
2. 生物学重复的表文件（由于后期比较的时候采用的配对显著性验证，因此会将每个样本之间进行比较，而存在生物学重复的类型需要通过该文件指明生物学重复情况）
此外，由于Trinity调用了多个R包，因此需要在进行运算前安装这些R包
> source("http://bioconductor.org/biocLite.R") #由于包都是Bioconductor上的，因此需要安装这个东西
> biocLite('edgeR')
> biocLite('limma')
> biocLite('DESeq2')
> biocLite('ctc')
> biocLite('Biobase')
> install.packages('gplots')
> install.packages('ape')
可以看到Trinity提供的这个脚本是可以支持edgeR、DESeq2、voom/limma以及ROTS的，也就是说后期需要指定使用哪个软件进行分析。
edgeR软件可以支持没有生物学重复的实验数据。通过指定dispersion参数进行设定，该参数很重要需要参考edgeR的使用说明进行设定。我们的数据是有生物学重复的因此采用以下脚本进行：
$cat samples.isoforms.txt
J J1.rsem
J J2.rsem
J J3.rsem
SF SF1.rsem
SF SF2.rsem
SF SF3.rsem
SM SM1.rsem
SM SM2.rsem
SM SM3.rsem
#查看样品信息三个处理三个生物学重复
$run_DE_analysis.pl #调用脚本(trinity自带)
--matrix isoforms.counts.matrix #采用的表达矩阵(脚本align_and_estimate_abundance.pl生成的)
--method edgeR #采用的比对方式
--samples_file samples.isoforms.txt
--output DEEdgeRout #结果输出到一个单独的文件夹
运行完了过后我们就去看看结果：
(bioinforspace) [spider 17:47:55 e[35;1m]/data1/spider/ytbiosoft/data/trinity.all/DEEdgeRout
$ll
total 19148
-rw-rw-r-- 1 spider spider 1258 Apr 12 09:37 isoforms.counts.matrix.J_vs_SF.J.vs.SF.EdgeR.Rscript
-rw-rw-r-- 1 spider spider 3719202 Apr 12 09:37 isoforms.counts.matrix.J_vs_SF.edgeR.DE_results
-rw-rw-r-- 1 spider spider 672993 Apr 12 09:37 isoforms.counts.matrix.J_vs_SF.edgeR.DE_results.MA_n_Volcano.pdf
-rw-rw-r-- 1 spider spider 1865558 Apr 12 09:37 isoforms.counts.matrix.J_vs_SF.edgeR.count_matrix
-rw-rw-r-- 1 spider spider 1258 Apr 12 09:37 isoforms.counts.matrix.J_vs_SM.J.vs.SM.EdgeR.Rscript
-rw-rw-r-- 1 spider spider 3956875 Apr 12 09:38 isoforms.counts.matrix.J_vs_SM.edgeR.DE_results
-rw-rw-r-- 1 spider spider 743109 Apr 12 09:38 isoforms.counts.matrix.J_vs_SM.edgeR.DE_results.MA_n_Volcano.pdf
-rw-rw-r-- 1 spider spider 2024265 Apr 12 09:38 isoforms.counts.matrix.J_vs_SM.edgeR.count_matrix
-rw-rw-r-- 1 spider spider 1264 Apr 12 09:38 isoforms.counts.matrix.SF_vs_SM.SF.vs.SM.EdgeR.Rscript
-rw-rw-r-- 1 spider spider 3946943 Apr 12 09:38 isoforms.counts.matrix.SF_vs_SM.edgeR.DE_results
-rw-rw-r-- 1 spider spider 723375 Apr 12 09:38 isoforms.counts.matrix.SF_vs_SM.edgeR.DE_results.MA_n_Volcano.pdf
-rw-rw-r-- 1 spider spider 1919916 Apr 12 09:38 isoforms.counts.matrix.SF_vs_SM.edgeR.count_matrix
一共12个文件共4类，包括：
1. Rscript（包括分析的脚本，通过这个脚本的修改可以调整参数和输出的图片）
2. cout_matrix （提取出来的两两比对的数据）
3. DE_results （这个是比对结果）
4. PDF （可视化的图片）
在此我们看一下这四类文件：
Rscript：脚本就不用看了
cout_matrix ：
DE_results：
PDF：
（MA plot请参考：https://www.jianshu.com/p/cdfac0bfb733）
接下来对差异数据进行注释和处理：
首先绘制差异基因的热图以及样本分布
调用analyze_diff_expr.pl脚本：
获取处理脚本的帮助文件
$analyze_diff_expr.pl
####################################################################################
Required:
--matrix|m <string> TMM.EXPR.matrix
Optional:
-P <float> p-value cutoff for FDR (default: 0.001)
-C <float> min abs(log2(a/b)) fold change (default: 2 (meaning 2^(2) or 4-fold).
--output <float> prefix for output file (default: "diffExpr.P${Pvalue}_C${C})
Misc:
--samples|s <string> sample-to-replicate mappings (provided to run_DE_analysis.pl)
--max_DE_genes_per_comparison <int> extract only up to the top number of DE features within each pairwise comparison.
This is useful when you have massive numbers of DE features but still want to make
useful heatmaps and other plots with more manageable numbers of data points.
--order_columns_by_samples_file instead of clustering samples or replicates hierarchically based on gene expression patterns,
order columns according to order in the --samples file.
--max_genes_clust <int> default: 10000 (if more than that, heatmaps are not generated, since too time consuming)
--examine_GO_enrichment run GO enrichment analysis
--GO_annots <string> GO annotations file
--gene_lengths <string> lengths of genes file
--include_GOplot optional: will generate inputs to GOplot and attempt to make a preliminary pdf plot/report for it.
##############################################################
跑如下代码：
```
$analyze_diff_expr.pl --matrix isoforms.TMM.EXPR.matrix --samples ../samples.isoforms.txt -P 0.001 -C 2
得到下面8个文件
-rw-rw-r-- 1 spider spider 4560 Apr 12 10:47 diffExpr.P0.001_C2.matrix
-rw-rw-r-- 1 spider spider 5023 Apr 12 10:47 diffExpr.P0.001_C2.matrix.R
-rw-rw-r-- 1 spider spider 54737 Apr 12 10:47 diffExpr.P0.001_C2.matrix.RData
-rw-rw-r-- 1 spider spider 10253 Apr 12 10:47 diffExpr.P0.001_C2.matrix.log2.centered.dat
-rw-rw-r-- 1 spider spider 10877 Apr 12 10:47 diffExpr.P0.001_C2.matrix.log2.centered.genes_vs_samples_heatmap.pdf
-rw-rw-r-- 1 spider spider 6620 Apr 12 10:47 diffExpr.P0.001_C2.matrix.log2.dat
-rw-rw-r-- 1 spider spider 1496 Apr 12 10:47 diffExpr.P0.001_C2.matrix.log2.sample_cor.dat
-rw-rw-r-- 1 spider spider 6821 Apr 12 10:47 diffExpr.P0.001_C2.matrix.log2.sample_cor_matrix.pdf
在这些结果包括
1. 两两比较的样本信息（.samples）
2. 样本比较后的基因表达情况上调还是下调的文件以及整合在一起的基因（.subset）
3. 差异基因在不同样本中分布数量（DE_feature_counts.P0.001_C2.matrix）
4. 差异基因表达量矩阵（diffExpr.P0.001_C2.matrix）
5. 差异基因对数转换后的表达矩阵 （diffExpr.P0.001_C2.matrix.log2.centered.dat）
6. 出的两个图（PDF）
看看两个PDF的图
我的图可能效果不大好，可能是样本的问题。可以看自己的图。
另外，还包括一个R脚本和一个Rdata（这两个文件可以通过修改和调用实现对后期出图的控制）
要提示一点就是注意对于R脚本修改后可以通过shell在所在文件夹进行直接运行后得到对应的图片
$Rscript diffExpr.P0.001_C2.matrix.R
通过edgeR这个R包对之前得到的表达量矩阵进行了差异表达验证，其中两个参数很重要一个是FDR和FoldChange，通过设定这两个值可以筛选出想要的差异基因。
最后还可以用脚本define_clusters_by_cutting_tree.pl做一个聚类分析的图(数据就是上一步生成的.matrix.RData文件)
$define_clusters_by_cutting_tree.pl --Ptree 60 -R diffExpr.P0.001_C2.matrix.RData
最后生成一个diffExpr.P0.001_C2.matrix.RData.clusters_fixed_P_60文件：
-rw-rw-r-- 1 spider spider 816 Apr 12 11:28 __tmp_plot_clusters.R
-rw-rw-r-- 1 spider spider 7384 Apr 12 11:28 my_cluster_plots.pdf
-rw-rw-r-- 1 spider spider 5407 Apr 12 11:28 subcluster_1_log2_medianCentered_fpkm.matrix
-rw-rw-r-- 1 spider spider 4924 Apr 12 11:28 subcluster_2_log2_medianCentered_fpkm.matrix
可以查看一下PDF文件
完毕！
欢迎交流哦：909474045@qq.com
