RNA-seq下游分析

这段时间，忙于工作和学习，没有时间更新文章，下面直接给大家带来一段RNA-seq下游分析的R语言代码，该代码主要是利用featureCounts软件对各样本数据统计的结果展开下游分析，代码如下：

##接featureCounts的基因表达水平的分析结果suppressMessages(library(DESeq2))setwd("F:/RNA_learning/")rm(list=ls())count_tab<-read.table("F:/RNA_learning/5.counts/counts",header=T)##删除count_tab的1、2列#count_tab<-count_tab[,-c(1,2)]sampleNames <-c('SRR6791080','SRR6791081','SRR6791082',                'SRR6791083','SRR6791084','SRR6791085')#给7-12列修改列名names(count_tab)[7:12]<-sampleNamescountMatrix<-as.matrix(count_tab[7:12])rownames(countMatrix)<-count_tab$Geneidhead(countMatrix)save(countMatrix,file='expr.Rdata')#deal是处理后样本，control是野生型，基因表达应该是control比上deal，要把control放在前面group_list <- c('control','control','control','deal','deal','deal')colData<-data.frame(sampleNames,group_list)colData$group_list = factor(colData$group_list,c('control','deal'))##colData用于存放样本信息数据，design是实验设计，表示counts文件中每个基因与colData中变量的依赖关系。dds<-DESeqDataSetFromMatrix(countData = countMatrix,                            colData = colData,                            design = ~ group_list)#过滤掉那些 count 结果为0的数据，这些基因没有表达dds <- dds[rowSums(counts(dds)) > 1,]dds2<-DESeq(dds)resultsNames(dds2) # lists the coefficientsres <- results(dds2, name="group_list_deal_vs_control")res <- res[order(res$padj),]resDF = as.data.frame(res)##PCAsuppressMessages(library(ggplot2))##PCA中归一化方法vsd <- vst(dds, blind=FALSE)pcaData <- plotPCA(vsd, intgroup=c("group_list"), returnData=TRUE)percentVar <- round(100 * attr(pcaData, "percentVar"))plotFig <- ggplot(pcaData, aes(PC1, PC2, color=group_list)) +  geom_point(size=3) +  xlab(paste0("PC1: ",percentVar[1],"% variance")) +  ylab(paste0("PC2: ",percentVar[2],"% variance")) +   coord_fixed()##保存图片ggsave(plotFig, filename = "yeast_DESeq2_PCA.pdf")##MApdf("yeast_DESeq2_MAplot.pdf")plotMA(res,ylim=c(-2,2))dev.off##筛选差异基因##按行筛选非NA的数据resDF_noNA <- resDF[complete.cases(resDF), ]##筛选padj<=0.05且log2FoldChange>2或<-2的行resDFfil = resDF_noNA[resDF_noNA$padj <= 0.05 &                         (resDF_noNA$log2FoldChange > 2 | resDF_noNA$log2FoldChange < -2)]##热图pheatmapnr_resDF=na.omit(resDF)suppressMessages(library("pheatmap"))##筛选差异基因,或者直接按paj.value升序排列，按自己需求选取#第一种方法select <- order(rowMeans(counts(dds2,normalized=TRUE)),decreasing=TRUE)[1:1000]#第二种方法select1 <- order(rowVars(counts(dds2,normalized=TRUE)),decreasing=TRUE)[1:1000]#归一化nt <- normTransform(dds2) # defaults to log2(x+1)#提取归一化后的基因log2.norm.counts <- assay(nt)[select,]df <- as.data.frame(colData(dds2))# pdf('heatmap1000.pdf',width = 6, height = 7)pheatmap(log2.norm.counts, cluster_rows=TRUE, show_rownames=FALSE,         cluster_cols=TRUE, annotation_col=df)