R语言通路基因ID提取表达table
<h1>背景</h1><p>用基因ID提取基因表达量，输出csv文件，作为下一步分析的输入文件。</p>
<h1>应用场景</h1>
<p>场景一：GO富集分析获得了GO term里的多个基因，需要获得这些基因的表达矩阵。</p>
<p>场景二：手里有已分类的基因列表，想提取相应的表达矩阵。</p>
<p>下一步，根据实际需要，用下面这些代码进行可视化：</p>
<ul>
<li>绘制heatmap，为单个GO term对应的基因进行聚类分析。</li>
<li>已分类的heatmap绘制多个GO term对应基因的heatmap</li>
<li>批量绘制单个基因的box plot</li>
<li>批量绘制表达谱曲线</li>
</ul>
<h1>输入文件</h1>
<p>需要两种文件：</p>
<ul>
<li>基因列表文件，<code>not_easy_input.txt</code>或 <code>GO*.txt</code></li>
<li>表达矩阵文件，<code>not_easy_input_expr.txt</code>或 <code>easy_input_expr.txt</code>。</li>
</ul>
<h1>基因表达矩阵</h1>
<p>每行一个基因，每列一个sample</p>
<pre><code class="language-{r}">exprSet<-read.table("not_easy_input_expr.txt",as.is = T)
exprSet
</code></pre>
<p>如果你不需要排序，可以跳过下面这段，直接进入“基因列表”。</p>
<p>下面按照表达量由高到低给基因排序：</p>
<pre><code class="language-{r,message=FALSE}">library(dplyr)
library(tidyr)
library(tibble)
exprSet <- exprSet %>%
#把行名变成一列，命名为symbol
rownames_to_column(var = "symbol")%>%
#新建一列rowMeans，每一行求平均值并填充进去
dplyr::mutate(rowMeans =rowMeans(.)) %>%
#按照求得的平均值把其他列从高到低排序
dplyr::arrange(desc(rowMeans)) %>%
#把symbol这一列变成行名
tibble::column_to_rownames(var = "symbol")%>%
#去掉平均值那一列
dplyr::select(-rowMeans)
write.table(exprSet[,1:10],"easy_input_expr.txt",quote = F)
</code></pre>
<p><code>easy_input_expr.txt</code>是按基因表达量从高到低排好序的表达矩阵。</p>
<h1>基因列表</h1>
<p>如果你的基因列表已经整理成 <code>GO*.txt</code>那样，就可以跳过这步，直接进入“开始提取”。</p>
<p>此处的输入文件是GO富集分析的输出文件：<code>not_easy_input_GO.txt</code>，每行一个GO term，GO terms对应的基因位于第8列，并且以“/”号分隔。</p>
<p>以下代码同样适用于其他富集分析结果，如果文件里基因名以“, ”分隔，只需把下面代码里的 <code>/</code>改为 <code>, </code>。</p>
<pre><code class="language-{r}">ego_BP_df<- read.table("not_easy_input_GO.txt",sep = "\t")
ego_BP_df
#输出排名靠前的GO term里的基因
#此处输出前两个
for (i in 1:2){
#此处用"/"分隔基因列表，这取决于你的基因ID之间用的是哪个分隔符
genelist <- unlist(strsplit(as.character(ego_BP_df$geneID),"/"))
write.table(genelist,paste0(rownames(ego_BP_df),".txt"),row.names = F,col.names = F,quote = F)
}
#或者提取特定的几个GO term
GO_id <- c("GO:0060333","GO:0050900") #把你想要提取的GO term放在这里
for (i in 1:length(GO_id)){
index <- grep(GO_id,ego_BP_df$ID) #在总GO term列表中的位置
genelist <- unlist(strsplit(as.character(ego_BP_df$geneID),"/"))
write.table(genelist,paste0(GO_id,".txt"),row.names = F,col.names = F,quote = F)
}
</code></pre>
<p>到这里，基因名就被整理成 <code>GO*.txt</code>文件的格式：</p>
<p>基因列表保存在多个文件里，每个文件包含一个GO term的基因ID。</p>
<p>基因名呈一列，每行一个基因名。</p>
<h1>开始提取</h1>
<h2>读入基因表达矩阵</h2>
<pre><code class="language-{r}">exprSet<-read.table("easy_input_expr.txt",as.is = T)
exprSet
</code></pre>
<h2>从基因列表文件逐一提取表达矩阵，保存到文件</h2>
<pre><code class="language-{r,warning=FALSE}">#按照文件名的规律，读取基因列表文件
#即使你只有一个基因ID列表文件，也可以这样操作
fnames<-Sys.glob("GO*.txt")
for (i in 1:length(fnames)){
genelist<-read.table(fnames)
#按表达矩阵中基因的顺序排列
genelist <- rownames(exprSet)
genelist_expr <- exprSet
write.csv(genelist_expr,paste0(unlist(strsplit(fnames,".txt")),".csv"),quote = F)
}
</code></pre>
<p>这里输出的 <code>GO*.csv</code>文件里就是基因表达矩阵。</p>
<p>通常情况下，到这里就结束了。除非：</p>
<p>你要用一条pheatmap命令画出多个GO term里的基因，就要继续运行下面的代码：</p>
<h2>重复出现的基因名处理</h2>
<p>例如你要画已分类heatmap，一步画出所有GO term里的基因表达谱heatmap。</p>
<p>会遇到报错提示“基因名不唯一”。</p>
<p>那是因为很多基因在多个GO term间重复出现，因此，用下面代码在基因名后面加上数字，以区分来源于不同GO term的同一基因。</p>
<pre><code class="language-{r}">fnames<-Sys.glob("GO*.txt")
for (i in 1:length(fnames)){
genelist<-read.table(fnames)
#按照表达矩阵中的位置排序
genelist <- rownames(exprSet)
genelist_expr <- exprSet
#在基因名后面加上“.数字”
rownames(genelist_expr)<-paste0(rownames(genelist_expr),".",i)
write.csv(genelist_expr,paste0(unlist(strsplit(fnames,".txt")),".csv"),quote = F)
}
</code></pre>
<p>文中所用数据可以关注公众号“生信喵实验柴”<br />
<img src="https://ipfs.io/ipfs/QmWMrgpjHGN5raAJdpML1fJL721mQNub87FsrS7bphNWFp" alt="1.png" /><br />
发送关键词“20250413”获取</p>
页:
[1]