wPGSA解析

#wPGSA
遺伝子発現から転写因子の活性化を予測するエンリッチメント解析です。
詳しくはプレスリリースを参照下さい。
Webツールとしての提供もしています。
現在はマウスのChIP-Seqデータのみ公開しています。

#準備
必要なデータとスクリプトをダウンロードします。

shell

git clone https://github.com/eiryo-kawakami/wPGSA

wPGSA.pyの実行にはnympy,scipy,rpy2が必要なのでインストールします。

shell

pip install numpy
pip install scipy
pip install rpy2

MacでMacの標準コンパイラがgccではなくclangらしく、-fopenmpオプションが使えないとエラーが出ました。Macの方はこのようなエラーが出るかもしれませんが、技巧をMac環境で使うを参考にしてインストールしてみてください。

#wPGSA実行
サンプルデータは/wPGSA/sample_logFC_file/にもありますが、
今回は、DESeq2を使った正規化で得られた出力データを使って実行してみます。

shell

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で変化の大きい転写因子を抽出してみます。
複雑なので、実行スクリプトだけ載せておきますので、参考にして下さい。

R

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 (i in c(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でヒートマップを作図しました。
複雑なので、実行スクリプトだけ載せておきますので、参考にして下さい。

R

threshold <- 5
cutoff <- function(num){
  if(num > threshold){
    return(threshold)
  }else if(num < -threshold){
    return(-threshold)
  }else return(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までの分化のデータを解析して見ましたが、薬剤処理の有無や遺伝子のノックアウト等のデータでも使用可能です。遺伝子レベルだけではなく、転写因子レベルの活性化を見比べることで、これまで解釈が難しかった変化を理解することができるかもしれません。