More than 5 years have passed since last update.
wPGSA
遺伝子発現から転写因子の活性化を予測するエンリッチメント解析です。
詳しくはプレスリリースを参照下さい。
Webツールとしての提供もしています。
現在はマウスのChIP-Seqデータのみ公開しています。
準備
必要なデータとスクリプトをダウンロードします。
git clone https://github.com/eiryo-kawakami/wPGSA
wPGSA.pyの実行にはnympy,scipy,rpy2が必要なのでインストールします。
pip install numpy
pip install scipy
pip install rpy2
MacでMacの標準コンパイラがgccではなくclangらしく、-fopenmpオプションが使えないとエラーが出ました。Macの方はこのようなエラーが出るかもしれませんが、技巧をMac環境で使うを参考にしてインストールしてみてください。
wPGSA実行
サンプルデータは/wPGSA/sample_logFC_file/にもありますが、
今回は、DESeq2を使った正規化で得られた出力データを使って実行してみます。
python PATH-TO-wPGSA/wPGSA/wPGSA.py --logfc-file PATH-TO-DESeq2_count_file --network-file PATH-TO-wPGSA/wPGSA/network/merged_mouse_v7_160212.network
・TF_wPGSA_p_value.txt
・TF_wPGSA_q_value.txt
・TF_wPGSA_t_score.txt
が出力されていると思います。
研究の目的に合わせて、適宜使い分けて下さい。
wPGSAの結果の可視化例
wPGSAで得られた結果を可視化する例を挙げておきます。
今回は得られた出力データからHSCとProBで変化の大きい転写因子を抽出してみます。
複雑なので、実行スクリプトだけ載せておきますので、参考にして下さい。
library(car)sample_info<-read.table("sample_info.txt",header=TRUE,row.name=1,sep='\t')#sample_info.txtはDESeq2の項目に記載してありますdf<-read.table("normalized_count_TF_wPGSA_t_score.txt",header=TRUE,row.name=1,sep='\t')pval<-c()for(iinc(1:nrow(df))){data<-na.omit(data.frame(t(df[i,]),sample_info[rownames(t(df[i,])),]))data<-data[data$sample_name=="Adult_BM_HSC"|data$sample_name=="Adult_BM_proB",]data$sample_name<-factor(data$sample_name,levels=c("Adult_BM_HSC","Adult_BM_proB"))colnames(data)<-c("value",colnames(sample_info))e<-try(est<-lm(value~sample_name,data=data),silent=FALSE)pval<-c(pval,coef(summary(est))[2,4])}df_ext<-data.frame(df,pval=pval,FDR=p.adjust(pval,method="BH"))write.table(df_ext,file="TF_wPGSA_t_score_ex.txt",sep="\t")作図
得られた出力データのFDRやp値を使って変化の大きい転写因子を抽出して作図します。
今回は、HSCとProBを比較をして得られたpvalで0.005以下の転写因子を抽出し、
hclustでクラスタリングをして並び替え、ggplot2でヒートマップを作図しました。
複雑なので、実行スクリプトだけ載せておきますので、参考にして下さい。
threshold<-5cutoff<-function(num){if(num>threshold){return(threshold)}elseif(num<-threshold){return(-threshold)}elsereturn(num)}library(ggplot2)library(amap)library(RColorBrewer)library(reshape2)library(ggthemes)myPalette<-colorRampPalette(rev(brewer.pal(name="PiYG",n=11)),space="Lab")df=read.table('TF_wPGSA_t_score_ex.txt',header=T,sep='\t',row.names=1)TF<-rownames(df[df$pval<0.005,])df<-df[df$pval<=0.005,]df<-df[,3:(ncol(df)-2)]df<-as.matrix(df)df<-na.omit(df)hc<-hclust(Dist(df,method="correlation"),method="ward.D2")df<-apply(df,c(1,2),cutoff)df<-df[hc$order,]# hc.t <- hclust(Dist(t(df),method="euclidean"),method="ward.D2")# df <- t(t(df)[hc.t$order,])# write.table(df,"VHvsK_KEGG_summary.ordered.txt",quote=FALSE,sep='\t')pdf="normalized_count_TF_wPGSA_t_score.pdf"data<-melt(df)(p<-ggplot(data,aes(as.factor(Var2),as.factor(Var1)))+geom_tile(aes(fill=value),colour="white")+scale_fill_gradientn(colours=myPalette(10),limit=c(-threshold,threshold))+theme(axis.text.x=element_text(angle=90,hjust=1))+ylab("TF (pval<0.005)")+xlab(""))ggsave(file=pdf,plot=p,width=4,height=8)
PAX5はBcellの代表的な転写因子ですが、ヒートマップの結果からproBで強く活性化していることが見られました。
また、BcellでMEIS1の不活性化が起こることも知られていますが、それについても今回の解析で確認することができました。
分化過程であるLMPPとCLPはほぼHSCとproBの中間のような転写制御パターンを示しましたが、replicate間でばらつきが見られました。これは分化過程において、一過的に不安定な状態を取っていると解釈できるのではないでしょうか。
まとめ
今回は例として、HSCからproBまでの分化のデータを解析して見ましたが、薬剤処理の有無や遺伝子のノックアウト等のデータでも使用可能です。遺伝子レベルだけではなく、転写因子レベルの活性化を見比べることで、これまで解釈が難しかった変化を理解することができるかもしれません。
Register as a new user and use Qiita more conveniently
- You get articles that match your needs
- You can efficiently read back useful information
- You can use dark theme