featureCounts是一個(gè)用來統(tǒng)計(jì)count數(shù)的軟件，運(yùn)行的速度飛快魄咕，比之前用的htseq-count快了好多好多。 照例先說一下怎么下載這個(gè)軟件：
wget https://jaist.dl.sourceforge.net/project/subread/subread-1.6.2/subread-1.6.2-Linux-x86_64.tar.gz
tar -zxvf subread-1.6.2-Linux-x86_64.tar.gz
cd subread-1.6.2-Linux-x86_64/bin
./featureCounts -h
然后來說這次的流程蚌父。 照舊用Hisat2來比對(duì)出Bam文件之后哮兰。 使用featureCounts統(tǒng)計(jì)：
featureCounts \
-T 16 \
-p \
-t exon \
-g gene_id \
-a Homo_sapiens.GRCh38.92.chr_patch_hapl_scaff.gtf \
-o all_feature.txt \
1.sort.bam \
2.sort.bam \
3.sort.bam \
4.sort.bam
# -T 使用的線程數(shù)
# -p 如果是paird end 就用
# -t 將exon作為一個(gè)feature
# -g 將gene_id作為一個(gè)feature
# -a 參考的gtf/gff
# -o 輸出文件
# 最后加上bam文件，有幾個(gè)就加幾個(gè)
然后會(huì)得到兩個(gè)文件苟弛，一個(gè)是結(jié)果喝滞，一個(gè)是結(jié)果的summary。 接下來就可以用DESeq2對(duì)結(jié)果進(jìn)行愉快的操作了膏秫。 使用R右遭。 我這次的樣本有6個(gè)。
library(DESeq2)
## 數(shù)據(jù)預(yù)處理
sampleNames <- c("10A_1", "10A_2", "10A_3", "7_1", "7_2", "7_3")
# 第一行是命令信息荔睹，所以跳過
data <- read.table("all_feature.txt", header=TRUE, quote="\t", skip=1)
# 前六列分別是Geneid Chr Start End Strand Length
# 我們要的是count數(shù)狸演，所以從第七列開始
names(data)[7:12] <- sampleNames
countData <- as.matrix(data[7:12])
rownames(countData) <- data$Geneid
database <- data.frame(name=sampleNames, condition=c("10A", "10A", "10A", "7", "7", "7"))
rownames(database) <- sampleNames
## 設(shè)置分組信息并構(gòu)建dds對(duì)象
dds <- DESeqDataSetFromMatrix(countData, colData=database, design= ~ condition)
dds <- dds[ rowSums(counts(dds)) > 1, ]
## 使用DESeq函數(shù)估計(jì)離散度，然后差異分析獲得res對(duì)象
dds <- DESeq(dds)
res <- results(dds)
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)
輸出兩個(gè)文件僻他，一個(gè)只有差異統(tǒng)計(jì)的結(jié)果宵距，一個(gè)包含各個(gè)樣本的結(jié)果。 這樣就完成了DESeq2了吨拗。
接下來是作圖： 一满哪、MA圖 MA圖是拿來展示數(shù)據(jù)表達(dá)是否異常，現(xiàn)在一般都不用了劝篷。
# library(DESeq2)
plotMA(res, main="DESeq2", ylim=c(-2, 2))
MA.png
二哨鸭、火山圖 可以非常直觀且合理地篩選出在兩樣本間發(fā)生差異表達(dá)的基因。
library(ggplot2)
# 這里的resdata也可以用res_des_output.csv這個(gè)結(jié)果重新導(dǎo)入哦娇妓。
# 現(xiàn)在就是用的前面做DESeq的時(shí)候的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=2, col="gray", lwd=0.5) +
geom_hline(yintercept=-log10(0.01), lty=2, col="gray", lwd=0.5)
valcano
valcano.png
三、PCA圖 就是主成分分析哈恰。是把數(shù)據(jù)降維之后進(jìn)行分析只估。PC1和PC2分別是貢獻(xiàn)率第一的主成分和貢獻(xiàn)率第二的主成分志群。
# 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
pca.png
四、熱圖 heatmap 實(shí)現(xiàn)這基因表達(dá)模式可視化的需求蛔钙。 從這里可以看到這6個(gè)樣本的表達(dá)差異锌云。
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)
heatmap.png
