Julia での VolcanoPlot を描くプログラムが見つからなかったので,以下を参考にして,プログラムを書いた。
1) は非常に参考になった。ただ近接するデータ点のアノテーションは Julia では難しいので,描画結果を時後に部分選択・拡大を行うことで代替することとした。
参考サイト
1) 19.11 Volcano plots
https://biocorecrg.github.io/CRG_RIntroduction/volcano-plots.html
2) Rでvolcano plotを描いてみる
https://qiita.com/tomo_metal/items/619e08f6dc18100884ee
3) RNA-seq解析 - R言語でvolcano plotを描く
https://qiita.com/shibanosuke/items/b4aa286993fd23a5c2ea
4) Pythonで作るVolcano plot
https://qiita.com/insilicomab/items/c382f10208e8bd2482b5
プログラム
function volcano_plot( # 最初の3つの引数さえ与えればよい。描画色の指定をするなら4番目の引数まで指定する。あとは,お好み次第。
log2FoldChange, # log2(変化倍率)
pvalue, # p 値
name; # 名前
color=[:black, :blue, :red], # 描画点の色, "NO", "DOWN", "UP" の順。':' に続いて一般的な色の名称を指定
markersize=3, # 描画点の大きさ
alpha = 0.5, # 描画点のα値
grid=false, # グリッドを書くか? false/true
tickfontsize=6, # 目盛りの文字サイズ
xlabel="log2FoldChange", # 横軸の名前
ylabel="-log10(pvalue)", # 縦軸の名前
guidefontsize=8, # 軸名の文字サイズ
vlinecolor=:red, # 追加の縦線の色の指定
vlinewidth=1, # 追加の縦線の線の幅の指定
vlinestyle=:dot, # :solid, :dash, :dot, :dashdot, :dashdotdot # 線の種類
hlinecolor=:red, # vline* に準ずる
hlinewidth=1,
hlinestyle=:dot, # :solid, :dash, :dot, :dashdot, :dashdotdot
annotatefontsize=8, # 名前のフォントサイズ
width=500, height=600, # 出力画像ファイルの横幅と高さをピクセル単位で指定
threshold_log2FoldChange = 0.6, # log2FoldChange の閾値 ± threshold_log2FoldChange
threshold_pvalue = 0.05, # 水平線の描画位置
leftright=15, # 名前 name を右揃えにするか揃えにするかの閾値
ypos = 999 # "DOWN"/"UP" を描く y 座標(デフォルトでは自動算出)
)
assigncolor(x) = x == "NO" ? color[1] : x == "DOWN" ? color[2] : color[3]
pyplot(grid=grid, size=(width, height), label="")
n = length(log2FoldChange)
diffexpressed = fill("NO", n);
diffexpressed[(log2FoldChange .> threshold_log2FoldChange) .& (pvalue .< threshold_pvalue)] .= "UP";
diffexpressed[(log2FoldChange .< -threshold_log2FoldChange) .& (pvalue .< threshold_pvalue)] .= "DOWN";
plt = scatter(log2FoldChange, -log10.(pvalue),
color=assigncolor.(diffexpressed),
markersize=markersize, markerstrokewidth=0, alpha=alpha,
tick_direction=:out, tickfontsize=tickfontsize,
xlabel=xlabel, ylabel=ylabel,
guidefontsize=guidefontsize)
vline!([-0.6, 0.6], color=vlinecolor, linewidth=vlinewidth, linestyle=vlinestyle)
hline!([-log10(0.05)], color=hlinecolor, linewidth=hlinewidth, linestyle=hlinestyle)
ypos = ypos == 999 ? (maximum(-log10.(de.pvalue))-log10.(0.05)) / 1.5 : ypos
annotate!(-threshold_log2FoldChange, ypos, text(" DOWN", annotatefontsize, color[2], :left))
annotate!( threshold_log2FoldChange, ypos, text("UP ", annotatefontsize, color[3], :right))
for i = 1:size(de, 1)
x = log2FoldChange[i]
y = -log10(pvalue[i])
d = diffexpressed[i]
if d != "NO"
annotate!(x, y,
text(" $(name[i]) ", 6, assigncolor(d),
y > leftright ? (d == "DOWN" ? :left : :right) : (d == "DOWN" ? :right : :left)))
end
end
plt
end
使用例
using RData, DataFrames, Plots
データ読み込み
# Download the data we will use for plotting
# download.file("https://raw.githubusercontent.com/biocorecrg/CRG_RIntroduction/master/de_df_for_volcano.rds", "de_df_for_volcano.rds", method="curl"
tmp = load("de_df_for_volcano.rds");de = dropmissing(tmp);
引数に与えるベクトルデータの設定
log2FoldChange = de[!, 2];
pvalue = de[!, 3];
name = de[!, 1];
関数呼び出し
plt = volcano_plot(log2FoldChange, pvalue, name)
savefig("fig4.png")
データ点が密になる箇所がある場合は,関数の戻り値を変数に代入し(例では plt),x 範囲,y 範囲を指定して再描画する。
この操作(範囲の再設定・再描画)は何回でも行うことができるので,描画結果を見ながら最適な描画範囲を設定できる。
plot(plt, xlims=(-0.9, -0.5), ylims=(6, 18), size=(400, 400))
savefig("fig5.png")
※コメント投稿者のブログIDはブログ作成者のみに通知されます