这一讲是miRNA-seq数据分析的分水岭，前面的5讲说的是读文献下载数据比对然后计算表达量，属于常规的流程分析，一般在公司测序之后都可以拿到分析结果，或者文献也会给出下载结果。但是单纯的分析一个样本意义不大，一般来说，我们做研究都是针对于不同状态下的miRNA表达量差异分析，然后做注释，功能分析，网络分析，这才是重点，也是难点。我这里就直接拿文献处理好的miRNA表达量来展示如何做下游分析，首先就是差异分析啦：根据文献，我们可以知道样本的分类情况是:

GSM1470353: control-CM, experiment1; Homo sapiens; miRNA-Seq   SRR1542714

GSM1470354: ET1-CM, experiment1; Homo sapiens; miRNA-Seq  SRR1542715

GSM1470355: control-CM, experiment2; Homo sapiens; miRNA-SeqSRR1542716

GSM1470356: ET1-CM, experiment2; Homo sapiens; miRNA-Seq SRR1542717

GSM1470357: control-CM, experiment3; Homo sapiens; miRNA-Seq SRR1542718

GSM1470358: ET1-CM, experiment3; Homo sapiens; miRNA-Seq SRR1542719

可以看到是6个样本的测序数据，分成两组，就是ET1刺激了CM细胞系前后对比而已！

### step8: differential expression analysis by R package for miRNA expression patterns:
## 文章里面提到的结果是：
MicroRNA sequencing revealed over 250 known and 34 predicted novel miRNAs to be differentially expressed between ET-1 stimulated and unstimulated control hiPSC-CMs.
## (FDR < 0.1 and 1.5 fold change)
rm(list=ls())
setwd('J:\\miRNA_test\\paper_results')  ##把从GEO里面下载的文献结果放在这里
sampleIDs=c()
groupList=c()
allFiles=list.files(pattern = '.txt')
i=allFiles[1]
sampleID=strsplit(i,"_")[[1]][1]
treat=strsplit(i,"_")[[1]][4]
dat=read.table(i,stringsAsFactors = F)
colnames(dat)=c('miRNA',sampleID)
groupList=c(groupList,treat)
for (i in allFiles[-1]){
sampleID=strsplit(i,"_")[[1]][1]
treat=strsplit(i,"_")[[1]][4]
a=read.table(i,stringsAsFactors = F)
colnames(a)=c('miRNA',sampleID)
dat=merge(dat,a,by='miRNA')
groupList=c(groupList,treat)
}

### 上面的代码只是为了把6个独立的表达文件给合并成一个表达矩阵
## we need to filter the low expression level miRNA
exprSet=dat[,-1]
rownames(exprSet)=dat[,1]
suppressMessages(library(DESeq2))
exprSet=ceiling(exprSet)
(colData <- data.frame(row.names=colnames(exprSet), groupList=groupList))

## DESeq2就是这么简单的用
dds <- DESeqDataSetFromMatrix(countData = exprSet,
colData = colData,
design = ~ groupList)
dds <- DESeq(dds)
png("qc_dispersions.png", 1000, 1000, pointsize=20)
plotDispEsts(dds, main="Dispersion plot")
dev.off()
res <- results(dds)
## 画一些图，相当于做QC吧
png("RAWvsNORM.png")
rld <- rlogTransformation(dds)
exprSet_new=assay(rld)
par(cex = 0.7)
n.sample=ncol(exprSet)
if(n.sample>40) par(cex = 0.5)
cols <- rainbow(n.sample*1.2)
par(mfrow=c(2,2))
boxplot(exprSet,  col = cols,main="expression value",las=2)
boxplot(exprSet_new, col = cols,main="expression value",las=2)
hist(exprSet[,1])
hist(exprSet_new[,1])
dev.off()library(RColorBrewer)
(mycols <- brewer.pal(8, "Dark2")[1:length(unique(groupList))])

# Sample distance heatmap
sampleDists <- as.matrix(dist(t(exprSet_new)))
#install.packages("gplots",repos = "http://cran.us.r-project.org")
library(gplots)
png("qc-heatmap-samples.png", w=1000, h=1000, pointsize=20)
heatmap.2(as.matrix(sampleDists), key=F, trace="none",
col=colorpanel(100, "black", "white"),
ColSideColors=mycols[groupList], RowSideColors=mycols[groupList],
margin=c(10, 10), main="Sample Distance Matrix")
dev.off()

png("MA.png")
DESeq2::plotMA(res, main="DESeq2", ylim=c(-2,2))
dev.off()
## 重点就是这里啦，得到了差异分析的结果
resOrdered <- res[order(res$padj),]
resOrdered=as.data.frame(resOrdered)
write.csv(resOrdered,"deseq2.results.csv",quote = F)

##下面也是一些图，主要是看看样本之间的差异情况
library(limma)
plotMDS(log(counts(dds, normalized=TRUE) + 1))
plotMDS(log(counts(dds, normalized=TRUE) + 1) - log(t( t(assays(dds)[["mu"]]) / sizeFactors(dds) ) + 1))
plotMDS( assays(dds)[["counts"]] )  ## raw count
plotMDS( assays(dds)[["mu"]] ) ##- fitted values.