R语言miRNA差异表达分析

Posted 2021-04-18 易学统计

     上篇介绍的是miRNA的数据下载和数据整理，本章介绍miRNA的差异表达分析，跟前面mRNA的表达分析大致类似，为节省篇幅，同样的注释不再累述了，大家可以参看。

01

样本分组

###样本分组load(file ='expr_df.Rda')metadata <- data.frame(names(expr_df)[-1])for (i in 1:length(metadata[,1])) { num <- as.numeric(substring(metadata[i,1],14,15)) if (num %in% seq(1,9)) {metadata[i,2] <- "tumor"} if (num %in% seq(10,29)) {metadata[i,2] <- "normal"}}#01-09表示肿瘤样本，10-19表示正常对照组织#加个名称names(metadata) <- c("TCGA_id","sample")#将sample转化成factormetadata$sample <- as.factor(metadata$sample)#我们可总结一下,有4例normal，179例肿瘤组织metadata %>% dplyr::group_by(sample) %>% summarise(n())save(metadata,file = "metadata.Rda")###保存分组信息

02

数据预处理

load(file ='expr_df.Rda')mycounts <- expr_dfmycounts=as.matrix(mycounts) ##DGEList对象要求矩阵格式rownames(mycounts)=mycounts[,1]exp=mycounts[,2:ncol(mycounts)]dimnames=list(rownames(exp),colnames(exp))data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)data=avereps(data) ##如果有重名的基因，就自动取平均值，方便后面的分析##一般miRNA没有重名的。data=data[rowMeans(data)>0,] ##删除在所有的样品上基本都不表达或者是低表达。##在做基因的时候有6万多数量大，可以范围限制严格点##而在miRNA条数才1881个，因此这里我们限制条件改成大于0即可。

     预处理前是1881个，删除在所有样品上基本不表达的，剩下个1449个。

   也可以采用cpm校准化的基因表达量作为过滤低表达基因的标准。

03

生成DGElist

group <- as.vector(metadata$sample) #按照癌症和正常样品数目修改design <- model.matrix(~group) ##变成一个实验设计对象y <- DGEList(counts=data,group=group) ##变成edgr包可以识别的数据结构，data需要是矩阵格式

    edgeR这个包将数据存储在DGEList对象中，需要指定的参数包括counts 和 group。其中表达矩阵存储在$counts中，样本分组信息等存储在$samples中。

    如果采用cpm过滤基因，R代码操作如下：

##也可以采用cpm校准化的基因表达量作为过滤低表达基因的标准。# cpm normalizationcountsPerMillion <- cpm(y)countCheck <- countsPerMillion > 0keep <- which(rowSums(countCheck) >= 1)degs.keep <- y[keep,]dim(degs.keep)## [1] 1449 183

04

 差异表达分析

      ①构建完DGEList对象，做因子校正

y <- calcNormFactors(y) ##对因子进行校正

      ②差异分析 dispersion(离散值)的估计

y <- estimateCommonDisp(y) ##估计变异系数 估计正常组内的变异系数y <- estimateTrendedDisp(y)y <- estimateTagwiseDisp(y) ##估计变异系数  估计癌症组内的变异系数

       ③：显著性差异检验：

et <- exactTest(y,pair = c("normal","tumor")) #检验两组时间的变异系数## ##类似配对T检验ordered_tags <- topTags(et, n=100000)##P值的校正 采用BH方法

     这个执行完差异表达的工作已经完成，后面的工作就是基于这个表格做一些基因的筛选和图片的绘制。

05

差异表达基因筛选

     ① 提取数据

allDiff=ordered_tags$tableallDiff=allDiff[is.na(allDiff$FDR)==FALSE,]diff=allDiffnewData=y$pseudo.countswrite.table(diff,file="allDiff.xls",sep="\t",quote=F)

       ② 筛选

foldChange=1 #差异在1倍以上padj=0.05 ##校正后的P值小于0.05
diffSig = diff[(diff$FDR < padj & (diff$logFC>foldChange | diff$logFC<(-foldChange))),]write.table(diffSig, file="diffSig.xls",sep="\t",quote=F)## 7 4###显示上调基因diffUp = diff[(diff$FDR < padj & (diff$logFC>foldChange)),]write.table(diffUp, file="up.xls",sep="\t",quote=F)## 0 4##显示下调基因diffDown = diff[(diff$FDR < padj & (diff$logFC<(-foldChange))),]write.table(diffDown, file="down.xls",sep="\t",quote=F)## 7   4###输出所有基因校正后的表达normalizeExp=rbind(id=colnames(newData),newData)write.table(normalizeExp,file="normalizeExp.txt",sep="\t",quote=F,col.names=F)
#输出差异基因校正后的表达diffExp=rbind(id=colnames(newData),newData[rownames(diffSig),]) # ##662 182write.table(diffExp,file="diffmRNAExp.txt",sep="\t",quote=F,col.names=F) #输出差异基因校正后的表达值（diffmRNAExp.txt）heatmapData <- newData[rownames(diffSig),] ##7 183

06

火山图绘制

pdf(file="vol.pdf")xMax=max(-log10(allDiff$FDR))+1yMax=12plot(-log10(allDiff$FDR), allDiff$logFC, xlab="-log10(FDR)",ylab="logFC", main="Volcano", xlim=c(0,xMax),ylim=c(-yMax,yMax),yaxs="i",pch=20, cex=0.4)diffSub=allDiff[allDiff$FDR<padj & allDiff$logFC>foldChange,]points(-log10(diffSub$FDR), diffSub$logFC, pch=20, col="red",cex=0.4)diffSub=allDiff[allDiff$FDR<padj & allDiff$logFC<(-foldChange),]points(-log10(diffSub$FDR), diffSub$logFC, pch=20, col="green",cex=0.4)abline(h=0,lty=2,lwd=3)dev.off()

07

热图

hmExp=log10(heatmapData+0.001)library('gplots')hmMat=as.matrix(hmExp)pdf(file="热图.pdf",width=60,height=90)par(oma=c(10,3,3,7))heatmap.2(hmMat,col='greenred',trace="none")dev.off()

08

总结

     本次完成miRNA差异的分析和两个常见图的绘制，下期见。

作者介绍：医疗大数据统计分析师,擅长R语言。

欢迎各位在后台留言，恳请斧正！

更多阅读

长按二维码

易学统计