今天開始差異分析的學習（開心），r什么的最熟悉了～

差異分析

1. 補充上節(jié)內(nèi)容

• STAR + featureCounts = STAR + HTSeq 升級版
• 安裝：conda install subread
• 優(yōu)點：非?？?/li>

#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. 表達定量結(jié)果轉(zhuǎn)換為表達矩陣

需要進行表達矩陣轉(zhuǎn)化的文件

• RSEM 自帶腳本

#差異表達
mkdir 06deseq_out

#構(gòu)建表達矩陣
rsem-generate-data-matrix *_rsem.genes.results > output.matrix

表達矩陣部分結(jié)果

• 去除所有樣本表達量為0的基因

#刪除未監(jiān)測到表達的基因
awk 'BEGIN{printf "geneid\ta1\ta2\tb1\tb2\n"}{if($2+$3+$4+$5>0)print $0}' output.matrix > deseq2_input.txt
#awk命令為知識盲區(qū)，需要補充

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)
#實際包含三步操作
#提取結(jié)果 
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    #上調(diào)表達
awk '{if($3 < 1 && $7 < 0.05) print $0}' diffexpr-results.txt    #下調(diào)表達