转录组入门学习（六）

今天开始差异分析的学习（开心），r什么的最熟悉了～

差异分析

1. 补充上节内容

• STAR + featureCounts = STAR + HTSeq 升级版
• 安装：conda install subread
• 优点：非常快

#featureCount: 基因定量
featureCounts -p -a ../00ref/Araport11_GFF3_genes_transposons.201606.gtf -o out_counts.txt -T 4 -t exon -g gene_id sample*_Aligned.sortedByCoord.out.bam

2. 表达定量结果转换为表达矩阵

需要进行表达矩阵转化的文件

• RSEM 自带脚本

#差异表达
mkdir 06deseq_out

#构建表达矩阵
rsem-generate-data-matrix *_rsem.genes.results > output.matrix

表达矩阵部分结果

• 去除所有样本表达量为0的基因

#删除未监测到表达的基因
awk 'BEGIN{printf "geneid\ta1\ta2\tb1\tb2\n"}{if($2+$3+$4+$5>0)print $0}' output.matrix > deseq2_input.txt
#awk命令为知识盲区，需要补充

3. edgeR

4. DESeq2

• 通过Bioconductor安装（略过）
• 如何使用

#准备工作::
#读取文件
input_data <- read.table("deseq2_input.txt", header = T, row.names = 1)

#取整
input_data <- round(input_data, digits = 0)

#preparing work:
input_data <- as.matrix(input_data)
condition <- factor(c(rep("ctr", 2), rep("exp", 2)))
coldata <- data.frame(row.names = colnames(input_data), condition)

#加载DESeq2包：
library("DESeq2")

#construct deseq2 matrix input:
dds <- DESeqDataSetFromMatrix(countData = input_data, colData = coldata, design = condition)
#conduce different expression analysis
dds <- DESeq(dds)
#实际包含三步操作
#提取结果 
res <- results(dds, alpha = 0.05)
summary(res)
res <- res[order(res$padj), ]
resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized = T)), by = "row.names", sort = F)
names(resdata)[1] <- "Gene"
#head(resdata)

#output the result file
write.table(resdata, file = "diffexpr-results.txt", sep = "\t", quote = F, row.names = F)

#可视化展示
#plotMA(res)

maplot <- function(res, thresh = 0.05, labelsig = T, ...) {
    with(res, plot(baseMean, log2FoldChange, pch = 20, cex = .5, log = "x", ...))
    with(subset(res, padj < thresh), points(baseMean, log2FoldChange, col = "red", pch = 20, cex = 1.5))
}
png("diffexpr-maplot.png", 1500, 1000, pointsize = 20)
maplot(resdata, main = "MA Plot")
dev.off()


install.packages("ggrepel")
library(ggplot2)
library(ggrepel)
resdata$significant <- as.factor(resdata$padj < 0.05 & abs(resdata$log2FoldChange) > 1)
ggplot(data = resdata, aes(x = log2FoldChange, y = -log10(padj), color = significant) +
  geom_point() +
  ylim(0, 8) +
  scale_color_manual(values = c("black", "red")) +
  geom_hline(yintercept = -log10(0.05), lty = 4, lwd = 0.6, alpha = 0.8) +
  geom_vline(xintercept = c(1, -1), lty = 4, lwd = 0.6, alpha = 0.8) +
  theme_bw() +
  theme(panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(colour = "black")) +
  labs(title = "Volcanoplot_biotrainee (by LiangYang)", x = "log2 (fold change)", y = "-log10 (padj)")+
  theme(plot.title = element_text(hjust = 0.5)) +
  geom_text_repel(data = subset(resdata, -log10(padj) >6), aes(label = Gene), col = "black", alpha = 0.8)
)  

diffexpr-maplot.png

• 提取差异基因

awk '{if($3 > 1 && $7 < 0.05) print $0}' diffexpr-results.txt    #上调表达
awk '{if($3 < 1 && $7 < 0.05) print $0}' diffexpr-results.txt    #下调表达