References：
1.RNA-seq(4)：Hisat2+FeatureCounts+DESeq2流程+作图！
https://pzweuj.github.io/2018/07/18/rna-seq-4.html
2.一个植物转录组项目的实战
http://www.bio-info-trainee.com/2809.html
3.RNA_seq（1）植物转录组差异基因分析简单练习
https://www.jianshu.com/p/7146d5c41294
Part I II 基本是照着做的，脚本会有一丢丢修改，之后吃饱了再弄吧
Part III:Featurecounts —>构建dds object ，获得expression matrix
featureCounts完成生物学定量
featureCounts是一款可以进行生物学定量的软件，它集成在subread的package里了
需要提供gtf格式的注释或者SAF格式的注释;
>gff3='/public/study/mRNAseq/tair/Arabidopsis_thaliana.TAIR10.28.gff3.gz'
>gtf='/public/study/mRNAseq/tair/Arabidopsis_thaliana.TAIR10.28.gtf.gz'
>featureCounts='/trainee/home/yjxiang/tools/subread-1.6.2-source/bin/featureCounts'
>nohup$featureCounts-T 5 -p -t exon -g gene_id -a$gtf-o /trainee/home/yjxiang/practice/subread_workflow/counts_out/counts_id.txt/trainee/home/yjxiang/practice/subread_workflow/mapping_out/*.bam & #挂在后台即便网络不稳也可执行，在提交程序和前台运行之间的选择！
DESeq2差异基因分析
First step:获取表达矩阵和分组信息
> library(DESeq2)
## 数据预处理
> sampleNames<-c('ERR1698194','ERR1698195','ERR1698196','ERR1698197')# 抽取前四列sample
>data <- read.table("count_id.txt", header=TRUE, quote="\t", skip=1)
>names(data)<-c('Geneid','Chr','Start','End','Strand','Length','ERR1698194','ERR1698195','ERR1698196','ERR1698197','ERR1698198','ERR1698199','ERR1698200','ERR1698201','ERR1698202','ERR1698203','ERR1698204','ERR1698205','ERR1698206','ERR1698207','ERR1698208','ERR1698209')#对 第一行重命名
# 前六列分别是Geneid Chr Start End Strand Length# 我们要的是count数，所以从第七列开始
>names(data)[7:10] <- sampleNames
>countData<-as.matrix(data[7:10])#第七列开始是每个sample对应的feature的counts,[前处sampleName命令有误，与此提取数据不match]
#将数据转换为矩阵格式
用featureCounts定量得到的counts_id.txt ，经过格式处理之后得到表达矩阵：countdata:第一列是基因ID，之后的列都是样本ID
每一行代表不同的基因在不同样本中的表达量.
> rownames(countData)<-data$Geneid
> database <- data.frame(name=sampleNames)
#database设置分组信息
>database <- data.frame(name=sampleNames, condition=c('a','a','b','b'))
>rownames(database) <- sampleNames#database是一个二维矩阵，赋予列名
>colnames(countData)<-NULL
##Second step: 构建dds对象
>dds <- DESeqDataSetFromMatrix(countData, colData=database, design= ~ condition)
> dds
class: DESeqDataSet
dim: 33602 4
metadata(1): version
assays(1): counts
rownames(33602): AT1G01010 AT1G01020 ...
ATCG01300 ATCG01310
rowData names(0):
colnames(4): ERR1698194 ERR1698197
ERR1698204 ERR1698207
colData names(2): name condition
> dds <- dds[ rowSums(counts(dds)) > 1, ]
## 使用DESeq函数估计离散度，然后差异分析获得res对象
> dds<-DESeq(dds)#对原始的dds进行标准化
>resultsNames(dds)#查看结果名称
>res <- results(dds)#用results函数提取结果，并赋值给res变量
> summary(res) #查看结果
out of 21129 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up) : 76, 0.36%
LFC < 0 (down) : 115, 0.54%
outliers [1] : 0, 0%
low counts [2] : 3687, 17%
(mean count < 12)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
>write.csv(res, "res_des_output.csv")
>resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)),by="row.names",sort=FALSE)
>write.csv(resdata, "all_des_output.csv", row.names=FALSE)
##Third step:提取结果并绘制火山图
Part VI: Drawing
> #提取差异基因!!!
> res <- res[order(res$padj),]
> resdata <-merge(as.data.frame(res),as.data.frame(counts(dds,normalize=TRUE)),by="row.names",sort=FALSE)
> deseq_res<-data.frame(resdata)
> up_diff_result<-subset(deseq_res,padj < 0.05 & (log2FoldChange > 1)) #提取上调差异表达基因
> down_diff_result<-subset(deseq_res,padj < 0.05 & (log2FoldChange < -1)) #提取下调差异表达基因
> write.csv(up_diff_result,"D:\\R.3.5.3\\上调_diff_results.csv") #输出上调基因
> write.csv(down_diff_result,"D:\\R.3.5.3\\下调_diff_results.csv") #输出下调基因
1.Valcano 火山图 可以非常直观且合理地筛选出在两样本间发生差异表达的基因
R script：
> library(ggplot2)
> # 这里的resdata也可以用res_des_output.csv这个结果重新导入哦。
> # 现在就是用的前面做DESeq的时候的resdata。
> resdata$change <- as.factor(
+ ifelse(
+ resdata$padj<0.01 & abs(resdata$log2FoldChange)>1,
+ ifelse(resdata$log2FoldChange>1, "Up", "Down"),
+ "NoDiff"
+ )
+ )#确定统计学显著性
> valcano <- ggplot(data=resdata, aes(x=log2FoldChange, y=-log10(padj), color=change)) +
+ geom_point(alpha=0.8, size=1) +
+ theme_bw(base_size=15) +
+ theme(
+ panel.grid.minor=element_blank(),
+ panel.grid.major=element_blank()
+ ) +
+ ggtitle("DESeq2 Valcano") +
+ scale_color_manual(name="", values=c("red", "green", "black"), limits=c("Up", "Down", "NoDiff")) +
+ geom_vline(xintercept=c(-1, 1), lty=5, col="gray", lwd=0.5) +
+ geom_hline(yintercept=-log10(0.01), lty=5, col="gray", lwd=0.5)
> valcano
2.PCA 主成分分析，把数据降维后进行分析，pc1和pc2是对整个数据集解释程度贡献率第一和第二的cluster，主成分。
R script:
> # library(ggplot2)
> rld <- rlog(dds)
> pcaData <- plotPCA(rld, intgroup=c("condition", "name"), returnData=T)
> percentVar <- round(100*attr(pcaData, "percentVar"))
> pca <- ggplot(pcaData, aes(PC1, PC2, color=condition, shape=name)) +
+ geom_point(size=3) +
+ ggtitle("DESeq2 PCA") +
+ xlab(paste0("PC1: ", percentVar[1], "% variance")) +
+ ylab(paste0("PC2: ", percentVar[2], "% variance"))
> pca
3.heatmap:热图实现基因表达模式可视化的需求。 4个样本的表达差异并不明显，因为我是从同一批次中随机抽取的无差别处理的4个samples。
Rscript:
> library(pheatmap)
> select <- order(rowMeans(counts(dds, normalized=T)), decreasing=T)[1:1000]
> nt <- normTransform(dds)
> log2.norm.counts <- assay(nt)[select,]
> df <- as.data.frame(colData(dds)[, c("name", "condition")])
> pheatmap(log2.norm.counts, cluster_rows=T, show_rownames=F, cluster_cols=T, annotation_col=df, fontsize=6)
>
