这一章主要是介绍RNA-seq的基本知识，以及简单的差异分析
RNA-seq
我们分析RNA-seq数据的时候，往往按照线图的步骤进行分析
欧老师在书中就RNA-seq是否去除dup，提出了争议，欧老师倾向于不去。
但我认为，具体情况得看你的QC报告以及你的实验目的了
目前在上游分析中：
1.mapping常用的软件有tophat（逐渐被淘汰了），hisat，STAR（这个速度很快），其中STAR我比较推荐，因为STAR用的是k-mer进行mapping的，并且错误率比较低，直接输出排好序的bam文件
2.计数上的主流软件有htseq-count，cufflinks里面的featureCount，steingtie，salmon（比较精确）
当然R里面也可以计数
library(Rsubread)
library(SCLCBam) ##bam文件
fl <- dir(getPathBamFolder(), "bam$", full.names = TRUE)
cnt <- featureCounts(files=fl, annot.inbuilt = "mm10")
##选取注释
3差异分析的R包主流的有三个edgeR，DESeq2和limma，当然也可以自己建立混合线性模型计算
差异表达分析
我们先来看edgeR
#载入edgeR。
library(edgeR)
# 载入数据
load(url("http://qiubio.com/bioconductor/RNA-seq/ds1.Rdata"))
#提取分组信息
grp <- as.factor(substr(colnames(counts), 1, 2))
#绘制样品间相关情况
o <- order(grp)
pairs(log2(1+counts[,o[1:7]]), pch=".",lower.panel=NULL)
#建立对象
d <- DGEList(counts=counts, group=grp)
d <- calcNormFactors(d)
#以cpm作为计数单位
cps <- cpm(d) ## count per million
k <- rowSums(cps>=1) > 2
d <- d[k,]
#MDS降维
cols <- as.numeric(d$samples$group)
plotMDS(d,col=cols)
接下来我饿么看下生物学重复所带的变异
plotBCV(d)
我们说的生物学重复是经过相同处理的不同生物学个体
这个表示生物学重复间的CV值与表达量之间的关系，每个点表示每个基因在每个生物学重复的表达
其中，表达量低的，在（同样处理）的生物学重复中的变异会增大，导致数据差异很大
DESeq2
library(DESeq2)
##我们载入与上一例相同的数据。
load(url("http://qiubio.com/bioconductor/RNA-seq/ds1.Rdata"))
ls()
head(counts[, 1:7], 3)
grp <- as.factor(substr(colnames(counts), 1, 2)) ##substr截取列名的前两个字母
dds <- DESeqDataSetFromMatrix(counts, colData=data.frame(grp), design=formula(~-1+grp)) ##这里的colData必须是一个DataFrame或者data.frame。每一行都对应着counts中的一列。design中的公式和limma中的方法一致。
design(dds)
## ~-1 + grp
dds <- DESeq(dds) ##一步到位，不需要太多步骤
res <- results(dds) ##得到结果
res
mcols(res, use.names=TRUE) ##结果中各例的意义
plotMA(res)
plotDispEsts(dds)
write.csv(as.data.frame(res), file="results.csv") ##输出结果
limma
library(limma)
v <- voom(d, mm, plot=TRUE) ##这里的d就是edgeR中的d，下面就是标准的limma流程了。
vf <- lmFit(v,mm) # 'mm' 前文定义了
cf <- contrasts.fit(vf,con) # 'con' 前文定义了
cf <- eBayes(cf)
topTable(cf)
