day49 转录组 差异基因表达分析
学习资料：B站视频：生物技能树 第八季转录组测序数据分析
得到表达矩阵之后就要考虑如何知道到差异表达基因（DEG），这个要在R中操作，具体是需要airway、DESeq2、edgeR这三个包。
airway这个包是数据包，但我觉得没啥用，它那些个命令只能用于它自己的数据，用自己从GEO上搞来的数据不香么？鸡肋。
DESeq2和EdgeR都是做基因差异表达分析的，主要用于RNA-Seq数据，有些研究者用来处理类似的ChIP-Seq, shRNA、质谱数据等，原理都是在不同处理因素的组别之间找差异。
一、安装R包
1，执行下面命令：
install.packages("BiocManager") #先安装biocmanager
library(BiocManager)
BiocManager::install ('airway')#这个包安装很费劲，报错，但实际上也没用。后悔安装之。
BiocManager::install ('DESeq2')
BiocManager::install ('edgeR')
Update all/some/none? [a/s/n]: 回答a。
2，安装airwa的时候出现报错：
image.png
解决办法一——失败：
remove.packages("cluster", lib="C:/Program Files/R/R-4.1.3/library")
remove.packages("MASS", lib="C:/Program Files/R/R-4.1.3/library")
remove.packages("Matrix", lib="C:/Program Files/R/R-4.1.3/library")
remove.packages("mgcv", lib="C:/Program Files/R/R-4.1.3/library")
remove.packages("nlme",lib="C:/Program Files/R/R-4.1.3/library")
remove.packages("survival", lib="C:/Program Files/R/R-4.1.3/library")
重新BiocManager::install ('airway')，依然报错。
解决办法二——成功：
下载该软件到本地，再从右下角的install按钮安装，如图：
image.png
3，卸载R包
由于后来发现airway没用，就练习了一下卸载。
remove.packages('airway')
二、文件准备
1，表达矩阵
通过前面计数后得到的txt文件，在excel里修整成为如下表，保留列名为样本名，行名为geneid。
image.png
2，分组准备
建立txt文件
image.png
三、读入并对表达矩阵进行过滤
rm(list = ls()) #首先清空R中的变量
setwd("D:/work/NGSanalysis/rnaseq")
a=read.table('D:/work/NGSanalysis/rnaseq/counts2.txt', header=T, row.names=1, stringsAsFactors = F) #行名，列名都设定了
index1=c(rowSums(a)>0) #删除都是零的基因删除
exprSet=a[index1,] #建立一个表达矩阵
write.csv(exprSet,'filter_reads_matrix.csv' )#过滤后的数据
四、差异分析：接上面步骤
library(DESeq2)
group=read.table('D:/work/NGSanalysis/rnaseq/group.txt', row.names=1,stringsAsFactors =T )#读入分组文件，为因子，有行名（样本名）
group_list=group[,1] #取这一列的数据为变量group_list
coldata <- data.frame(row.names =colnames(exprSet), group_list) #创建一个数据框coldata，只有一列内容是分组信息，行名为样本名，
dds <- DESeqDataSetFromMatrix(countData = exprSet, colData = coldata, design = ~group_list) #这个命令创建DEseqDataSet（dds）。countData用于说明数据来源，colData用于说明不同组数据的实验操作类型，design用于声明自变量，即谁和谁进行对比
dds2 <- DESeq(dds) #用于对dds数据进行运算及分析,共三步：文库大小估计；离散程度估计；统计检验
res <- results(dds2,contrast = c("group_list","siSUZ12","con")) #生成结果文件,siSUZ12组vs. con组。
resOrdered<- res[order(res$padj),] #order函数是给出从小到大排序后的位置(默认升序), 将结果文件按照res文件中的padj这一列进行降序排列，其中$符号表示res中的padj这一列
DEG=as.data.frame(resOrdered)
write.csv(DEG,file = "treat_vs_con.csv") #在工作目录中生成结果文件
#按照自定义阈值提取差异基因
DEGsig <-subset(res, padj < 0.001 & abs(log2FoldChange) > 1)
write.csv(DEGsig,file= "treat_vs_con_sig.csv")
image.png
五、用DESeq2归一化的表达矩阵
rld <- rlogTransformation(dds2) ## 得到经过normlization的表达矩阵。
exprSet_new=assay(rld) #生成矩阵文件
write.csv(exprSet_new,'Normalized_matrix.csv' ) #导出数据，保存为csv格式。
说明
rlogTransformation是DESeq2包中的一个命令（函数），在help文档中查看，它能把原始含有reads的表达数值，根据整个样本和每个基因的长度等，计算出相对表达。并表示为log2格式，因此最后导出的数据里会有负数。比如-1.7，就是2^（-1.7）=0.3。
assay是SummarizedExperiment包中的函数，可以把阵提取出来。
