Rなどを見よう見まねで勉強しています。備忘録と共有のためにブログに記録することとしました。@kazuigarashi
by kazuhiko_igarashi
カレンダー
S |
M |
T |
W |
T |
F |
S |
|
|
|
|
|
1
|
2
|
3
|
4
|
5
|
6
|
7
|
8
|
9
|
10
|
11
|
12
|
13
|
14
|
15
|
16
|
17
|
18
|
19
|
20
|
21
|
22
|
23
|
24
|
25
|
26
|
27
|
28
|
29
|
30
|
31
|
カテゴリ
以前の記事
最新の記事
記事ランキング
タグ
フォロー中のブログ
画像一覧
メモ帳
最新のトラックバック
ライフログ
検索
ブログパーツ
外部リンク
ファン
ブログジャンル
その他のジャンル
|
2021年 09月 05日
主成分分析の原理をシミュレーションしてみます。 nヶの変数からなるデータについて、固有ベクトルを使って線形変換して、nヶの新しい変数(分散の大きいものから小さいものまで順番に並んだ、すなわち、固有値の大きいものから小さいものに並んだ)を作り直し、それを使ってプロットします。 この方法の前に、まずはRの主成分分析の関数prcomp()を使ってみます。 用いる二次元正規分布のデータの作成法は こちらの記事を参考下さい。 作成したデータを使って、二つの方法で主成分分析を行います。一つ目は、Rのprcomp()関数を使う方法、二つ目は、分散共分散行列をcor()関数で作成し、その固有ベクトルをeigen()関数で求め、そのベクトルを使って線形変換で主成分分析を行う方法です。(ほぼ)同じ結果が得られます。 #相関係数xを与える対角行列の値を出力する関数(詳しくは こちらの記事) com_norm <- function(x){ #xは設定したい相関係数とする b = 1/x + ((1/x)**2-1)**0.5#ここで右項の+は本当は+/-の二通りを取り得ることに注意。 return(b) }#相関係数がおおよそxとなる二次元乱数を100ヶ作成する関数(詳しくは こちらの記事) #わかりやすくするためデータの各行毎に線形変換していますが、この記事の後半にあるように、行列の積で一度に求めることもできます。 Synth <- function(x){X1 <- rnorm(100, mean = 0, sd = 1) #平均0、標準偏差1の正規分布から乱数100ヶ X2 <- rnorm(100, mean = 0, sd = 1)#平均0、標準偏差1の正規分布から乱数100ヶ XX <- cbind(X1, X2) #二つの値から成るデータセットにする A <- matrix(c(1, x, x, 1), nrow = 2, ncol = 2) #入力値xを使って線形変換に使う対角行列を作る com <- c() for (i in 1:100){ a <- A%*%XX[i,] #XXの各行毎に二値を線形変換します。 com <- cbind(com, a) } return(com) }#相関係数がおおよそ0.9になる線形変換の対角行列の数値を出します DiagMatVal <- com_norm(0.9)#その数値を使って二次元正規分布を作ります data <- Synth(DiagMatVal)data2 <- t(data) #転置しておきます head(data2)#作成したデータをそれぞれ平均0、標準偏差1に標準化しておく data3 <- scale(data2) #scale()関数を使用 mean(data3[,1])sd(data3[,1])mean(data3[,2])sd(data3[,2])#得たデータの分布を確認しておく plot(density(data3[,1]), xlim = c(-6,6), ylim = c(0, 0.5)) #data3[,1]について密度分布をプロットする curve(dnorm(x, 0, 1), add = T, lty = 2, col = 2) #平均0、標準偏差1の正規分布を赤点線で上書き #Synth()でサンプル数を10000ヶにすればほとんど正規分布となること確認できます。 #データをプロットして確認 plot(data3[,1], data3[,2])cor(data3[,1], data3[,2]) #設定した相関係数に近い分布が得られている #これでシミュレーションに使うデータdata3ができました #例1、主成分分析の関数prcomp()を使って軸を回転させる(分散が最大の軸、次に最大の軸を作ってプロットし直す) pca <- prcomp(data3) #scale = TRUEとすると標準偏差で割って標準化する。今回は例2の固有ベクトルを使った例と比較するために標準化しない。 print (pca) #縦方向に固有ベクトルが出てくる PCAvalues <- pca$x #主成分得点を出力する head(PCAvalues)plot(PCAvalues, xlim = c(-5,5), ylim = c(-5,5), xlab = 'PC1', ylab = 'PC2') #主成分得点を使ってプロットする、x軸y軸のスケールを同じにすること #例2、主成分分析を線形変換で実施する一連の過程をコードしてみます。分散共分散行列から固有ベクトルを計算し、それを使って二つの軸を設定します(分散が最大の軸、次に最大の軸、を決定します)。 varcovar <- cor(data3, method = 'pearson') #分散共分散行列を作成します、methodで分散共分散の計算方法を選ぶことができます varcovar#固有値と固有ベクトルを得ます。 ei <- eigen(varcovar) #固有ベクトルは縦方向が一つのセットになる ei#固有ベクトルは行方向に格納されている(列方向に並んでいる) #各要素へのアクセス ei$values[1] #第一の固有値 ei$values[2] #第二の固有値 ei$vectors[,1] #行方向に格納されている固有値1に対応するベクトル ei$vectors[,2] #行方向に格納されている固有値2に対応するベクトル #固有ベクトルを使ってdata3を線形変換する rotation2 <- t(ei$vectors) #転置しないとPC1が縦軸になってしまう PCAvalues2 <- rotation2 %*% t(data3)PCAvalues2 <- t(PCAvalues2) #転置して行方向にデータをならべる head(PCAvalues2 )plot(PCAvalues2, xlim = c(-5,5), ylim = c(-5,5), xlab = 'PC1', ylab = 'PC2') #図としては(ほとんど)同じ結果が得られます。 #第一の方法と第二の方法で基底ベクトルが逆方向になることがあり、その場合は、PC1やPC2の展開が対称になります。 #(ほとんど)同じ結果であることは以下の数値からわかります。(場合によっては完全に同じにはならないようです。これは統計処理の詳細が異なるためと思われます) summary(PCAvalues) #第一の方法の結果 summary(PCAvalues2) #第二の方法の結果 sd(PCAvalues[,1]) #第一の方法の結果 sd(PCAvalues[,2]) #第一の方法の結果 sd(PCAvalues2[,1]) #第二の方法の結果 sd(PCAvalues2[,2]) #第二の方法の結果 二変数で記述されていたサンプルをPC1の変数一つで記述できそうです。多次元データをPC1~PC2の二変数(あるいはPC3までの三変数)でプロットするなど、よく使われる方法です。 分散共分散行列の固有値と固有ベクトルは、分散共分散行列で表現される関係性を作った力について、その大きさと方向を示していると捉えることもできます。例えば遺伝子発現の場合、無秩序な状態(遺伝子間に相関がない状態)に多数の様々な「制御の力」が作用して相関が生じると考えてみます。そうすると、その状態を作り出した制御の力、それぞれについて、その方向を固有ベクトルは示している、力の大きさを対応する固有値が示している、と捉えることもできます。その力の大きなものから(寄与度の高いものから)順に使ってデータを説明する、とも言えそうです。
#
by kazuhiko_igarashi
| 2021-09-05 22:08
| 統計
2020年 03月 21日
「事始め」なんでRを始めたのか? https://datator.exblog.jp/23084502/Rを勉強してよかったことhttps://datator.exblog.jp/23322591/ コードアカデミーhttps://datator.exblog.jp/23089512/プログラムのダウンロードhttps://datator.exblog.jp/23095223/R関連パッケージのダウンロードhttps://datator.exblog.jp/23103814/「基本操作」オブジェクトとその操作、ベクトルvectorとマトリクスmatrix、データフレームdata framehttps://datator.exblog.jp/23112006/表データの読み込み (read.table()) https://datator.exblog.jp/23118193/基本統計量と検定https://datator.exblog.jp/23142769/1次元散布図(beeswarm(), RColorBrewer) https://datator.exblog.jp/23151925/棒グラフhttps://datator.exblog.jp/23156497/棒グラフ(横) https://datator.exblog.jp/23193328/集合解析 その1 (VennDiagram, venn.diagram()など) https://datator.exblog.jp/23210931/集合解析 その2 (union(), intersect(), setdeff(), write.table(), など) https://datator.exblog.jp/23244165/集合解析 その3 (色指定) https://datator.exblog.jp/23252423/3群の積集合をリスト化する (sets, set_intersection()) https://datator.exblog.jp/24214521/樹形図の作成 (hclust()) https://datator.exblog.jp/23303990/樹形図の作成 その2、クラスター毎のデータの取り出し方 (rect.hclust() ) https://datator.exblog.jp/23308181/画像サイズの調整(画像の保存) (png(),dev.off() ) https://datator.exblog.jp/23417549/ヒートマップの作成、遺伝子発現パターンの可視化 (gplots, heatmap.2(), hclust() ) https://datator.exblog.jp/23736047/ヒートマップの色指定の仕方(RColorBrewer)(display.brewer.all(), brewer.pal()) https://datator.exblog.jp/23739335/single cell PCRのヒートマップの例、遺伝子発現パターンの可視化 (gplots, hclust(),colorRampPalette()) https://datator.exblog.jp/23940279/single cell PCRのヒートマップの例、その2 (hclust()) https://datator.exblog.jp/24281698/ヒートマップ、色変化の方向を変える (RColorBrewer, brewer.pal(), rev(), png(), barplot()) https://datator.exblog.jp/29970800/散布図(scatter plot)でデータポイントを見る、ラベル入れる (identify()) https://datator.exblog.jp/24397313/散布図(scatter plot)でデータポイントの列名(変数名)や補助線を入れる (maptools, pointLabel(), par()) https://datator.exblog.jp/24464336/分子間相互作用のネットワーク図をigraphで描いてみる (igraph) https://datator.exblog.jp/24523122/分子間相互作用ネットワークをタンパク質遺伝子機能で色分けしてみる (igraph, RColorBrewer, ) https://datator.exblog.jp/24552429/「データ操作」merge()関数を用いたデータの結合、列の入れ換え、属性(項目名)つけ (names(), merge(), ) https://datator.exblog.jp/24580887/reshape2、dplyrパッケージを使ったデータ処理 (reshape2, dplyr, melt(), summarize(), ) https://datator.exblog.jp/24814124/「ggplot2を使った可視化」ggplot2を使ったグラフ作成(棒グラフ)(ggplot2, apply(), cbind(), ggplot(), geom_bar(), geom_errorbar() ) https://datator.exblog.jp/24796622/ggplot2を使ったグラフ(棒グラフ)の体裁の調整 (ggplot(), geom_bar(), geom_errorbar(), coord_cartesian(), theme()) https://datator.exblog.jp/24800112/ggplot2を使ったグラフ(棒グラフ)(実験グループ毎の表示、エラーバーのズレを調整)(reshape2, dplyr, ggplot(), geom_bar(), geom_errorbar, position_dodge(), ) https://datator.exblog.jp/24820291/ggplot2で棒グラフと散布図の組み合わせをつくってみる (aggregate(), ggplot2, ggplot(), geom_boxplot(), geom_bar(), geom_col(), stat_summary(), geom_point()) https://datator.exblog.jp/27008169/ggplot2を使ったグラフ作成(折れ線、時系列)(reshape2, melt(), ggplot(), geom_line(), geom_point(), ) https://datator.exblog.jp/24805889/ggplot2を使ったグラフ作成(折れ線、時系列、エラーバー)(dplyr, ggplot(), geom_line(), geom_point(), geom_errorbar(), scale_y_continuous(), theme_classic(), coord_cartesian()) https://datator.exblog.jp/24809125/ggplot2で線形変換の写像をプロットするプログラム (ggplot2, matrix(), %*%, rbind(), geom_point(), geom_line()) https://datator.exblog.jp/25972346/ggplot2で線形変換の写像をプロットするプログラム〜その2〜 (ggplot2, matrix(), %*%, rbind(), geom_point(), geom_line()) https://datator.exblog.jp/25978889/「モンテカルロシミュレーション」モンテカルロシミュレーションに向けた準備 (for(), runif()) https://datator.exblog.jp/25376795/モンテカルロ法で面積割合いと円周率を求めてみる (for(), runif(), sqrt()) https://datator.exblog.jp/25381470/ggplot2でモンテカルロ法シミュレーション結果を図にしてみる (ggplot2, ggplot(), stat_binhex()) https://datator.exblog.jp/25385771/モンテカルロ法シミュレーションで誕生日パラドックスを検討してみる、その1(ChIPseqのピークの分析へ??)) (duplicated(), as.integer(), runif()) https://datator.exblog.jp/25431928/モンテカルロ法シミュレーションで誕生日パラドックスを検討してみる、その2(ChIPseqのピークの分析へ??)) (duplicated(), as.integer(), runif()) https://datator.exblog.jp/25438493/モンテカルロ法シミュレーションでChIPシークエンスのピークのオーバーラップを考える (intersect(a), runif()) https://datator.exblog.jp/25831140/モンテカルロ法シミュレーションでChIPシークエンスのピークのオーバーラップを考える2 〜確率分布の推定〜 (ggplot2, ggplot(), geom_point(), intersect(a), runif()) https://datator.exblog.jp/25834627/正規分布に従う乱数の線形変換で2次元正規分布モデルを得るプログラム (rnorm(), %*%) https://datator.exblog.jp/25986972/正規分布に従う乱数の線形変換で2次元正規分布モデルを得るプログラム〜その2、関数の統合とggplot2で可視化など〜 (rnorm(), cbind(), %*%, ggplot2, ggplot(), geom_point(), cor()) https://datator.exblog.jp/25993516/二項分布とボルツマン分布 (rbinom(), hist(), runif(), as.integer(), ggplot2, ggplot(), geom_histogram, geom_vline()) https://datator.exblog.jp/26330166/カイ二乗分布のシミュレーション (runif(), rnorm(), hist(), ggplot2, ggplot(), geom_histogram()) https://datator.exblog.jp/27124636/カイ二乗分布をモンテカルロシミュレーションプログラムで遊ぶ (rnorm(), sum(), for(), ) https://datator.exblog.jp/27127208/RとImageMagickでアニメーション作成、カイ二乗分布シミュレーションをアニメに (animation, ggplot2, ggplot(), geom_histogram(), ggtitle(), paste(), gifplot(), ani.replay(), ani.options()) https://datator.exblog.jp/27130635/カイ二乗分布シミュレーションのアニメーションを改善 (animation, ggplot2, ggplot(), geom_histogram(), ggtitle(), paste(), gifplot(), ani.replay(), ani.options()) https://datator.exblog.jp/27152947/カイ二乗検定のシミュレーション (seq(), for(), hist()) https://datator.exblog.jp/27233623/カイ二乗検定のシミュレーション結果をggplot2で描画 (ggplot2, data.frame(), ggplot, geom_histogram(), stat_function(), geom_vline(), order()) https://datator.exblog.jp/27257026/ブートストラップ法を試す (matrix(), for(), rnorm(), replicate(), arrows(), sample()) https://datator.exblog.jp/27984979/p値の問題を考える (rnorm(), hist(), sample(), t.test()) https://datator.exblog.jp/29490861/「文字列操作の基本とパスワード作成プログラム」文字列の操作、パスワード自動生成プログラムの準備 (sample(), paste(), ) https://datator.exblog.jp/26842675/文字列の操作、パスワード自動生成プログラムの準備〜その2 (sample(), append(), paste()) https://datator.exblog.jp/26851715/パスワード自動生成プログラムを完成させます (for(), runif(), sample(), append(), paste(), sprintf()) https://datator.exblog.jp/26874013/組み合わせの爆発をシミュレーションでモデリングして確率を推定するプログラム、パスワードを例に (subset(), which(), sample(), paste(), length()) https://datator.exblog.jp/26951345/「反応シミュレーション」deSolveパッケージを使ったシミュレーション計算 (deSolve) https://datator.exblog.jp/24736820/酵素反応のシミュレーションをコードしてdeSolveパッケージで解析する(deSolve) https://datator.exblog.jp/24867771/酵素反応のシミュレーションをプログラム的なものに改変してみる (deSolve) https://datator.exblog.jp/24873940/酵素反応シミュレーションプログラム、基質濃度と速度定数を振るプログラム (deSolve, for()) https://datator.exblog.jp/25172679/酵素反応シミュレーションプログラム、基質濃度と複数の速度定数を振るプログラム (deSolve, for()) https://datator.exblog.jp/25187521/「その他」NGS解析用ゲノムブラウザGenomeJackにRefSeqデータを取り込むhttps://datator.exblog.jp/28405046/
#
by kazuhiko_igarashi
| 2020-03-21 19:08
| 目次
2020年 03月 21日
ヒートマップの色をRColorBrewerで作る際、色変化の方向を変更したいことがあります。そこで色コードのベクトルを逆転させて使ってみます。 library(RColorBrewer)mypalette <- brewer.pal(11, "RdBu")mypalette #[1] "#67001F" "#B2182B" "#D6604D" "#F4A582" "#FDDBC7" "#F7F7F7" "#D1E5F0" "#92C5DE" "#4393C3" "#2166AC" "#053061" このまま使うと、数値小さい方が赤、大きい方が青になります。そこで、逆転させたベクトルを作ります。 mypalette2 <- rev(mypalette)mypalette2# [1] "#053061" "#2166AC" "#4393C3" "#92C5DE" "#D1E5F0" "#F7F7F7" "#FDDBC7" "#F4A582" "#D6604D" "#B2182B" "#67001F"パレットを画像にして保存しておきます。 png("color_pallete.png", 400, 50)par(mar = c(0,0,0,0))barplot(rep(1, 11), col=mypalette2, axes=F)dev.off()heatmap.2で描画するさいにcol= mypalette2 と指定します。例えば heatmap.2(as.matrix(y), Rowv = TRUE, Colv= NA, distfun = dist, hclustfun = hclust, dendrogram = c("row"), col=mypalette2, trace="none")ヒートマップについてはこちらをご覧下さい。
#
by kazuhiko_igarashi
| 2020-03-21 08:41
| ヒートマップ
2019年 06月 25日
帰無仮説有意差検定におけるp値の誤用やp値ハッキングがいろいろ話題になっています。p値の使用を止めた心理学系ジャーナルもあるようです。実験科学者としてはその利点と盲点を理解した上で活用したいものです。Rを使うと、p値分布のシミュレーションが簡単にできるので、より理解を深めることができます。 まず、帰無仮説H 0が真である場合、すなわち、同じ母集団に由来する標本が二つある場合、そのt検定で得られるp値は一様分布になることをシミュレーションで確認します。 正規分布母集団をつくり、そこから2回独立にサンプリングし、そのt検定を行います。これを繰り返してp値をシミュレーションします。 population <- rnorm(1000, mean = 100, sd = 30) #正規分布する母集団、要素数1000、平均100、標準偏差30をつくる hist(population, col = "#ff00ff40", border = "#ff00ff", breaks=20)p.val <- c() #p値を収納する空ベクトル for (i in 1:5000){ #5000回繰り返します s1 <- sample(population, 50, replace=TRUE) #n=50の実験とします s2 <- sample(population, 50, replace=TRUE) res <- t.test(s1, s2, var.equal=TRUE) p <- res$p.value #テスト結果からp値を取り出します。 p.val <- append(p.val, p) #p値をベクトルに収納します }print(p.val)plot(p.val, main = "p value when Ho is true, 5000 trials", col = "#ff00ff40")hist(p.val, breaks = 20, freq = FALSE, col = "#ff00ff40", border = "#ff00ff", main = "distribution of p value when Ho is true, 5000 trials")帰無仮説H 0が真であっても、p値が有意水準(多くの場合はp = 0.05)を超えて小さくなる場合があることがわかります。真の帰無仮説を棄却し、二つの標本に差があると誤って述べる疑陽性、第一種過誤です。有意水準は、許容される第一種過誤の最大限の確率αを示します。 次に帰無仮説H 0が偽である場合(標本が異なる母集団に由来する)、そのt検定のp値の分布をシミュレーションします。この例は母集団平均に1.1倍の違いがある場合です(生化学的にはかなり微妙な差でしょうか、○○学ではよくあるのかもしれません)。サンプルサイズを変えて比較してみます。 pop1 <- rnorm(1000, mean = 100, sd = 30) #正規分布の母集団1をつくる pop2 <- rnorm(1000, mean = 110, sd = 30) #正規分布の母集団2をつくる 、平均は1.1倍 hist(pop1, col = "#ff00ff40", border = "#ff00ff", breaks=20, main = "two populations with 1.1-fold difference in the means")hist(pop2, col = "#0000ff40", border = "#0000ff", breaks=20, add = T)p.val <- c()for (i in 1:5000){s1 <- sample(pop1, 50, replace=TRUE) #n=50の実験とします s2 <- sample(pop2, 50, replace=TRUE) res <- t.test(s1, s2, var.equal=TRUE) p <- res$p.value p.val <- append(p.val, p) }plot(p.val, main = "p value when Ho is false, 5000 trials", col = "#0000ff40")hist(p.val, breaks = 20, freq = FALSE, col = "#0000ff40", main = "distribution of p value when Ho is false, sample size = 50, 5000 trials")帰無仮説が偽であってもp値は確率的に分布することが改めて確認できます。p値が0.05よりも大きくなる場合(偽の帰無仮説を間違って採択し、サンプル間に差がないと判断してしまう)が意外に多いことがわかります。これが偽陰性、第二種過誤です。βはその確率、1-βは検出力です。 サンプル数が小さいと、βが増え検出力が低下することを確認します。直上のコードでサンプリング数を5にしてみます。 s1 <- sample(pop1, 5 replace=TRUE) #n=5の実験とします s2 <- sample(pop2, 5, replace=TRUE) たまたまある実験で有意差が出たとしても再現実験では否定される確率の方が高い、そういう状況になったり論文にしてしまわないよう、気をつける必要があります。 母集団の違いの大きさ(効果)や標準偏差、サンプリングサイズをいろいろ振ってみると面白いですし、危険性も実感できます。 p値は、「帰無仮説が真であった場合」、観察されたデータとそれ以上に偏ったデータが得られる確率とされています。帰無仮説が真である確率を示すものではありませんし、(1-p)が標本間に差がある確率を示すわけでもありません。
#
by kazuhiko_igarashi
| 2019-06-25 19:09
| 統計
2018年 06月 25日
DRY解析教本でトレーニングを受けたとき、NGSデータを可視化するソフトとして GenomeJackを使いました。教本の230ページからインストール法、そして遺伝子アノテーションとしてRefSeq Genesを読み込み表示する方法が紹介されています。RefSeqGenesの表示を最近試したところ、うまくデータをダウンロードできませんでした。そこで別の方法を工夫したのでまとめておきます。 UCSCのマウスゲノム(mm9)データダウンロードサイトからターミナルを使ってtxtファイルをダウンロードしてゲノムブラウザ用に加工します。以下はマウスの例ですが、 ヒトゲノム(hg19)データダウンロードサイトなどから同様にできるはずです。 そもそもの問題は、ゲノムブラウザで使うファイル形式、gene transfer fromat (GTF)をどうやって作るのか?でした。 以下のサイトの解説読むと、 Table BrowserからGTF形式でダウンロードしても遺伝子名は失われていて、このままではGenomeJackに読み込んでも遺伝子名は表示できませんでした。バグも含まれる可能性ありとも書いてあります。 http://genomewiki.ucsc.edu/index.php/Genes_in_gtf_or_gff_format“GenePred is a table format commonly used for gene prediction tracks in the UCSC Genome Browser. The genePredToGtf command-line utility can be used to convert genePred to GTF. While the Table Browser does contain an option to output query results in GTF, the output is limited, and in some cases, may contain bugs. The best method to convert genePred to GTF is the genePredToGtf command-line utility.” やったこと コマンドラインから作業用ディレクトリに入って以下のコマンドでgenePredToGtf command-line utilityをダウンロードします。 $ rsync -aP rsync://hgdownload.cse.ucsc.edu/genome/admin/exe/macOSX.x86_64/genePredToGtf ./(注 UCSCサーバは大学ネットからは接続できずmobile wifiで接続してダウンロードしました。最初はmacOSXではなくlinuxを指定したのでプログラムが動かなかった!!) $ ./genePredToGtf でユーティリティーの説明が出てくることを確認します。 ちなみにパスの通っていないコマンドの実行にはパスを通す必要があります。カレントディレクトリにおかれているgenePredToGtf であれば ./genePredToGtf のように。 genePredToGtfの使い方はこちらにあります。 http://genomewiki.ucsc.edu/index.php/Genes_in_gtf_or_gff_format次にターミナルからmm9のrefGene.txtをダウンロードします。 $ wget --timestamping 'ftp://hgdownload.cse.ucsc.edu/goldenPath/mm9/database/refGene.txt.gz' -O refGene.txt.gzカレントディレクトリにgzファイルがダウンロードされるのでクリックして解凍。ターミナルコマンドでもいいです。 次にターミナルから加工します。 まず最初の列(ビン)を取り除き $ cut -f 2- refGene.txt > refGene.inputついでgenePredToGtfを用いてgtfに変換します。 $ ./genePredToGtf file refGene.input mm9refGene.gtfつぎに染色体、座標でソートします。ソートしないとGenomeJackに読み込めない!!パイプでコマンドをつなぎます。いわゆるパイプライン。 $ cat mm9refGene.gtf | sort -k1,1 -k4,4n > mm9refGene_sorted.gtfこのファイルをgenomejackのlocal storage > import > gtf で読み込みます。 この時、途中set parametersのウインドウでgene_nameをラベルに選びdata typeはstringにします。これで遺伝子名が表示されるようになります。
#
by kazuhiko_igarashi
| 2018-06-25 12:15
| NGS等、その他
|