瀑布流网站,iis 浏览网站,泾川县建设局网站,上海网站建设公司联系方式新手遇到的问题都是类似的#xff0c;比如批量ID转换虽然我写过大量的教程#xff1a;ID转换大全 不过都需要R基础#xff0c;因为是大批量转换啊#xff01;但热心肠的植物生物信息学教学大佬还是友善的给出了解决方案我也狗尾续貂制作了一个网页工具教程#xff1a;简… 新手遇到的问题都是类似的比如批量ID转换虽然我写过大量的教程ID转换大全 不过都需要R基础因为是大批量转换啊但热心肠的植物生物信息学教学大佬还是友善的给出了解决方案我也狗尾续貂制作了一个网页工具教程简单的3个步骤不会代码也可以很容易把ID批量转换啦不过有趣的是我搜索电脑资料看到了2年前写的拟南芥教程。不过我为什么会花时间写拟南芥教程呢1 首先加载必要的包library(ggplot2)library(stringr)# source(https://bioconductor.org/biocLite.R)# biocLite(clusterProfiler)# biocLite(org.At.tair.db)library(clusterProfiler)library(org.At.tair.db)2 然后导入数据deg1read.table(https://raw.githubusercontent.com/jmzeng1314/my-R/master/DEG_scripts/tair/DESeq2_DEG.Day1-Day0.txt)deg1na.omit(deg1)head(deg1)## baseMean log2FoldChange lfcSE stat pvalue## AT3G01430 22.908929 18.989990 2.148261 8.839704 9.597263e-19## AT1G47405 20.709551 20.958806 2.434574 8.608820 7.381677e-18## AT2G33830 1938.159722 -2.560609 0.312663 -8.189678 2.619266e-16## AT5G28627 8.118376 -21.131318 2.875691 -7.348257 2.008078e-13## AT2G33750 9.789740 -19.989301 2.847513 -7.019915 2.220033e-12## AT3G54500 2238.314652 2.720430 0.386305 7.042182 1.892517e-12## padj## AT3G01430 1.858318e-14## AT1G47405 7.146571e-14## AT2G33830 1.690562e-12## AT5G28627 9.720602e-10## AT2G33750 7.164418e-09## AT3G54500 7.164418e-09prefixDay1-Day03 然后判断显著差异基因很明显这个时候用padj来挑选差异基因即可不需要看foldchange了。table(deg1$padj0.05)## ## FALSE TRUE ## 19166 197table(deg1$padj0.01)## ## FALSE TRUE ## 19303 60diff_geneList rownames(deg1[deg1$padj0.05,])up_geneList rownames(deg1[deg1$padj0.05 deg1$log2FoldChange 0,])down_geneList rownames(deg1[deg1$padj0.05 deg1$log2FoldChange 0,])length(diff_geneList)## [1] 197length(up_geneList)## [1] 89length(down_geneList)## [1] 1083.1 画个火山图看看挑选的差异基因合理与否colnames(deg1)## [1] baseMean log2FoldChange lfcSE stat ## [5] pvalue padjlog2FoldChange_Cutof 0## 这里我不准备用log2FoldChange来挑选差异基因仅仅是用padj即可deg1$change as.factor(ifelse(deg1$padj 0.05 abs(deg1$log2FoldChange) log2FoldChange_Cutof, ifelse(deg1$log2FoldChange log2FoldChange_Cutof ,UP,DOWN),NOT))this_tile Cutoff for log2FoldChange is ,round(log2FoldChange_Cutof,3), \nThe number of up gene is ,nrow(deg1[deg1$change UP,]) , \nThe number of down gene is ,nrow(deg1[deg1$change DOWN,]))g_volcano ggplot(datadeg1, aes(xlog2FoldChange, y-log10(padj), colorchange)) geom_point(alpha0.4, size1.75) theme_set(theme_set(theme_bw(base_size20))) xlab(log2 fold change) ylab(-log10 p-value) ggtitle( this_tile ) theme(plot.title element_text(size15,hjust 0.5)) scale_colour_manual(values c(blue,black,red)) ## corresponding to the levels(res$change)print(g_volcano)4 GO/KEGG注释4.1 首先进行KEGG注释diff.kk - enrichKEGG(gene diff_geneList,organism ath,pvalueCutoff 0.99,qvalueCutoff0.99)kegg_diff_dt - as.data.frame(setReadable(diff.kk,org.At.tair.db,keytype TAIR))up.kk - enrichKEGG(gene up_geneList,organism ath,pvalueCutoff 0.99,qvalueCutoff0.99)kegg_up_dt - as.data.frame(setReadable(up.kk,org.At.tair.db,keytype TAIR))down.kk - enrichKEGG(gene down_geneList,organism ath,pvalueCutoff 0.99,qvalueCutoff0.99)kegg_down_dt - as.data.frame(setReadable(down.kk,org.At.tair.db,keytype TAIR))4.2 可视化看看KEGG注释结果## KEGG patheay visulization: down_kegg$pvalue0.05,];down_kegg$group-1 up_kegg$pvalue0.05,];up_kegg$group1 datrbind(up_kegg,down_kegg) colnames(dat)## [1] ID Description GeneRatio BgRatio pvalue ## [6] p.adjust qvalue geneID Count group dat$pvalue -log10(dat$pvalue) dat$pvaluedat$pvalue*dat$group datdat[order(dat$pvalue,decreasing F),] g_kegg geom_bar(statidentity) scale_fill_gradient(lowblue,highred,guide FALSE) scale_x_discrete(name Pathway names) scale_y_continuous(name -log10P-value) coord_flip() ggtitle(Pathway Enrichment) print(g_kegg)4.3 接着进行GO注释for (onto in c(CC,BP,MF)){ ego OrgDb org.At.tair.db, keyType TAIR, ont onto , pAdjustMethod BH, pvalueCutoff 0.2, qvalueCutoff 0.9) ego TAIR) write.csv(as.data.frame(ego),paste0(prefix,_,onto,.csv)) #ego2 ego2ego pdf(paste0(prefix,_,onto,_barplot.pdf)) pbarplot(ego2, showCategory12)scale_x_discrete(labelsfunction(x) str_wrap(x,width20))print(p)dev.off() }ggsave(filename paste0(prefix,_volcano_plot.pdf),g_volcano)## Saving 7 x 5 in imageggsave(filename paste0(prefix,_kegg_plot.pdf),g_kegg)## Saving 7 x 5 in imagewrite.csv(x kegg_diff_dt,file paste0(prefix,_kegg_diff.csv))write.csv(x kegg_up_dt,file paste0(prefix,_kegg_up.csv))write.csv(x kegg_down_dt,file paste0(prefix,_kegg_down.csv))5 基因ID注释library(org.At.tair.db)ls(package:org.At.tair.db)## [1] org.At.tair org.At.tair.db ## [3] org.At.tairARACYC org.At.tairARACYCENZYME## [5] org.At.tairCHR org.At.tairCHRLENGTHS ## [7] org.At.tairCHRLOC org.At.tairCHRLOCEND ## [9] org.At.tairENTREZID org.At.tairENZYME ## [11] org.At.tairENZYME2TAIR org.At.tairGENENAME ## [13] org.At.tairGO org.At.tairGO2ALLTAIRS ## [15] org.At.tairGO2TAIR org.At.tairMAPCOUNTS ## [17] org.At.tairORGANISM org.At.tairPATH ## [19] org.At.tairPATH2TAIR org.At.tairPMID ## [21] org.At.tairPMID2TAIR org.At.tairREFSEQ ## [23] org.At.tairREFSEQ2TAIR org.At.tairSYMBOL ## [25] org.At.tair_dbInfo org.At.tair_dbconn ## [27] org.At.tair_dbfile org.At.tair_dbschema## Then draw GO/kegg figures:deg1$gene_idrownames(deg1)id2symboltoTable(org.At.tairSYMBOL) deg1merge(deg1,id2symbol,bygene_id)## 可以看到有一些ID是没有gene symbol的这样的基因意义可能不大也有可能是大部分人漏掉了head(deg1)## gene_id baseMean log2FoldChange lfcSE stat pvalue## 1 AT1G01010 58.25657 1.13105335 0.8000274 1.4137683 0.1574300## 2 AT1G01010 58.25657 1.13105335 0.8000274 1.4137683 0.1574300## 3 AT1G01020 542.64394 -0.05745554 0.3721650 -0.1543819 0.8773086## 4 AT1G01030 48.77442 -1.09845263 1.2875862 -0.8531100 0.3935983## 5 AT1G01040 1708.68949 0.32421734 0.2777530 1.1672865 0.2430947## 6 AT1G01040 1708.68949 0.32421734 0.2777530 1.1672865 0.2430947## padj change symbol## 1 0.6008903 NOT ANAC001## 2 0.6008903 NOT NAC001## 3 0.9805661 NOT ARV1## 4 0.8144858 NOT NGA3## 5 0.6992279 NOT DCL1## 6 0.6992279 NOT EMB60可以看到基因的ID和symbol的对应关系就出来了根使用网页工具是类似的感兴趣的朋友可以试试看网页工具和R代码的ID批量转换差别有多大。■ ■ ■ 全国巡讲约你第1-10站北上广深杭西安郑州 吉林武汉和成都(全部结束)七月份我们不外出只专注单细胞系统学习单细胞分析报名生信技能树的线下培训手慢无。一年一度的生信技能树单细胞线下培训班火热招生全国巡讲第11站-港珠澳专场(生信技能树爆款入门课)
