featureCounts是一个用来统计count数的软件，运行的速度飞快，比之前用的htseq-count快了好多好多。 照例先说一下怎么下载这个软件：
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来比对出Bam文件之后。 使用featureCounts统计：
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 使用的线程数 # -p 如果是paird end 就用 # -t 将exon作为一个feature # -g 将gene_id作为一个feature # -a 参考的gtf/gff # -o 输出文件 # 最后加上bam文件，有几个就加几个
然后会得到两个文件，一个是结果，一个是结果的summary。 接下来就可以用DESeq2对结果进行愉快的操作了。 使用R。 我这次的样本有6个。
library(DESeq2) ## 数据预处理 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数，所以从第七列开始 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 ## 设置分组信息并构建dds对象 dds <- DESeqDataSetFromMatrix(countData, colData=database, design= ~ condition) dds <- dds[ rowSums(counts(dds)) > 1, ] ## 使用DESeq函数估计离散度，然后差异分析获得res对象 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)
输出两个文件，一个只有差异统计的结果，一个包含各个样本的结果。 这样就完成了DESeq2了。
接下来是作图： 一、MA图 MA图是拿来展示数据表达是否异常，现在一般都不用了。
# library(DESeq2) plotMA(res, main="DESeq2", ylim=c(-2, 2))
二、火山图 可以非常直观且合理地筛选出在两样本间发生差异表达的基因。
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=2, col="gray", lwd=0.5) + geom_hline(yintercept=-log10(0.01), lty=2, col="gray", lwd=0.5) valcano
三、PCA图 就是主成分分析。是把数据降维之后进行分析。PC1和PC2分别是贡献率第一的主成分和贡献率第二的主成分。
# 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
四、热图 heatmap 实现这基因表达模式可视化的需求。 从这里可以看到这6个样本的表达差异。
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)
