今天開(kāi)始差異分析的學(xué)習(xí)（開(kāi)心）芒珠，r什么的最熟悉了～
差異分析
1. 補(bǔ)充上節(jié)內(nèi)容
- STAR + featureCounts = STAR + HTSeq 升級(jí)版
- 安裝：conda install subread
- 優(yōu)點(diǎn)：非呈厶迹快
#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. 表達(dá)定量結(jié)果轉(zhuǎn)換為表達(dá)矩陣
需要進(jìn)行表達(dá)矩陣轉(zhuǎn)化的文件
- RSEM 自帶腳本
#差異表達(dá)
mkdir 06deseq_out
#構(gòu)建表達(dá)矩陣
rsem-generate-data-matrix *_rsem.genes.results > output.matrix
表達(dá)矩陣部分結(jié)果
- 去除所有樣本表達(dá)量為0的基因
#刪除未監(jiān)測(cè)到表達(dá)的基因
awk 'BEGIN{printf "geneid\ta1\ta2\tb1\tb2\n"}{if($2+$3+$4+$5>0)print $0}' output.matrix > deseq2_input.txt
#awk命令為知識(shí)盲區(qū)，需要補(bǔ)充
3. edgeR
4. DESeq2
- 通過(guò)Bioconductor安裝（略過(guò)）
- 如何使用
#準(zhǔn)備工作::
#讀取文件
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)
#實(shí)際包含三步操作
#提取結(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)表達(dá)
awk '{if($3 < 1 && $7 < 0.05) print $0}' diffexpr-results.txt #下調(diào)表達(dá)
