对day49，day50两次学习的总结。
拿到reads表达矩阵后，在自己电脑上用DESeq2这个R包进行数据分析。
一、先要准备表达矩阵和分组表格
1，表达矩阵
count之后的txt文件进行整理，只留下geneid和四个样本列。
其中有一些行不完整，有NA的数据，在后面读入的时候一定要用fill=TURE参数。
2，分组表格
可以用txt设置分组表格，但是可能是用WPS建立表格，转成txt有格式问题，一直报错：“Error in make.names(col.names, unique = TRUE) :
invalid multibyte string at '<ff><fe>s'
In addition: Warning messages:
1: In read.table(file = file, header = header, sep = sep, quote = quote, :
line 1 appears to contain embedded nulls”
最后解决办法是，转成scv格式，用read.csv 命令读入。group_list一列的列名必须是“group_list”。
二、在R中编辑，依此运行下面的脚本：
#step1: read counts
rm(list = ls()) #首先清空R中的变量
setwd("/Users/meraner/Documents/bioinfo")
a=read.table('/Users/meraner/Documents/bioinfo/counts2.txt', fill=T, header=T, row.names=1, stringsAsFactors = F)
a = na.omit(a) #含有NA值的行删除
index1=c(rowSums(a)>0)
exprSet=a[index1,]
write.csv(exprSet,'filter_reads_matrix.csv' ) #过滤后的数据保存
说明
由于count文件中可能会有很一些行里面有缺省NA，就是没比对成功的。所以用fill=TURE，即可读入，这些数据显示为NA。在Mac电脑的R里面需要用fill=T啊！不然会报错： object 'TURE' not 。
接着可以删除有NA的行。即用no.omit命令把NA的行去除，避免后面分析中的干扰。这些NA数据的基因名通常都很奇怪，很长一串后续分析有点儿烦，去除更好。
index1为逻辑型。就是计算每一行的和，得到的数大于零就是ture，等于零就是false，再把ture的那些行单独拎出来建立一个expreSet表达框，就是过滤过后的数据，用来进行后面的分析。
#step2: set group & build dds
library(DESeq2)
coldata <- read.csv ('/Users/meraner/Documents/bioinfo/group.csv', row.names=1, stringsAsFactors =T)
group_list=coldata[,1] #取这一列的数据为变量group_list
colnames(exprSet)<-rownames(coldata) #把表达矩阵的列名重设为简单格式且和coldata里的样本名保持一致
dds <- DESeqDataSetFromMatrix(countData = exprSet, colData = coldata, design = ~group_list) #这个命令创建DEseqDataSet（dds）。
说明
创建一个数据框coldata，只有一列内容，即分组信息。行名为样本名，顺序必须和前面表达矩阵从左到右的顺序一致。
建立一个变量group_list（名字必须是这个），里面内容就是每个样本的分组信息，顺序必须和前面表达矩阵从左到右的顺序一致。这个变量还必须是一个因子型的。所以用了stringsAsFactors =T。
countData用于说明数据来源，colData用于说明不同组数据的实验操作类型，design用于声明自变量，即谁和谁进行对比
#step3: DEG
dds2 <- DESeq(dds)
res <- results(dds2,contrast = c("group_list","siSUZ12","con"))
resOrdered<- res[order(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")
说明
用DESeq命令对dds数据进行运算及分析,共三步：文库大小估计；离散程度估计；统计检验。生成结果文件是按照siSUZ12组vs. con组进行比较的。#order函数是给出从小到大排序后的位置(默认升序), 将结果文件按照res文件中的padj这一列进行降序排列，其中$符号表示res中的padj这一列。再按照自定义阈值提取差异基因。
#step4：normalization matrix
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包中的函数，可以把阵提取出来。
#step5: visulization
#fig1
hist(exprSet_new)
#fig2
tiff(filename = "exprBOX.tiff",width = 1500, height = 1500, units = "px", res = 150)
par(cex = 1.2) #指定符号的大小
par(mar=c(10, 5, 5, 5))
n.sample=ncol(exprSet) #样本数
if(n.sample>40) par(cex = 0.5)
n.sample=ncol(exprSet)#样本数
cols <- rainbow(n.sample*1.2)
boxplot(exprSet_new, col = cols,main="expression value",las=2)
dev.off()
#fig3
sigGene= read.csv('treat_vs_con_sig.csv',row.names=1)
for(i in rownames(sigGene)){
genename=paste(i,sep="") #基因名
pdf(file=paste(i,'_counts.pdf',sep="")) #以pdf格式输出
plotCounts(dds2, genename, intgroup = "group_list")
dev.off()
}
#fig4
tiff(filename = "PCA.tiff", width = 1500, height = 1500, units = "px", res = 150)
plotPCA(rld, intgroup="group_list")
dev.off()
说明
hist() 用来做直方图（柱形图），应该用归一化的数据来看。
boxplot是用来看每个样本的整体表达水平。
plotCounts是用来看单个感兴趣基因的表达差异
plotPCA用来做PCA聚类分析
生成的图片都保存在工作目录中，就是前面setwd的那个路径。
