パールのモジュールをローカル権限のLinuxにインストールする方法
-
Perl -MCPAN -e shell
-
ユーザーのホームディレクトリーに ~/.cpan/CPAN/MyConfig.pm というのができてる。
-
その中の
makepl_arg の部分を以下のように変更します。
'makepl_arg' => q[LIB=~/perl/lib
-
ログインシェルの環境変数にさきほど指定した ~/perl/lib の場所を PERL5LIB として設定しておく必要があります。
.bashrcの中身を
export PERL5LIB=$HOME/perl/lib:$PERL5LIB
とする。
-
Perl -MCPAN -e shell とすると、install Hogeとかいろいろなものをインストールできる。
Rのpackageをローカル権限のLinuxにインストールする方法
-
mkdir $HOME/R/library
-
$HOME/.Renvironを作り、その中に「 R_LIBS_USER="~/R/library"」と書く
-
その後は、R CMD INSTALL ダウンロードしたプログラム --library=./R/library/
-
または、install.packages("name-of-your-package",lib="~/R/library")
Generic Genome Browser, v2のインストール
- 徹底的に、ここに書いてあるモジュール等をインストール。
- http://gmod.svn.sourceforge.net/gmod/ をマニュアルでインストール、いろいろしたが、マニュアルじゃないエラーでインストールできなかった。
- 完全に、Defaultでインストール。いろいろ変えるといろいろ問題が起こった。つまり、perl Makefile.PL、./Build test、./Build
install
- 重要なのは、./Build apache_confで、でてくるものを、”http.conf”に加えること。以下のも、”http.conf”に加える。常に、apachectlでStop、Startを繰り返して様子を見る。あと、GBrowse's
temporary data である[/var/tmp/gbrowse2]をchmod a+wとして書き込み可能にする。
ScriptAlias /cgi-bin/ "/var/www/cgi-bin/"
<Location "/cgi-bin">
AllowOverride None
Options None
Order allow,deny
Allow from all
</Location>
HMMERの使い方
- phmmer： Search a sequence against a sequence database. (BLASTP-like)
- jackhmmer： Iteratively search a sequence against a sequence database. (PSIBLAST-like)
- hmmbuild： Build a profile HMM from an input multiple alignment.
- hmmsearch： Search a profile HMM against a sequence database.
- hmmscan： Search a sequence against a profile HMM database.
- hmmalign： Make a multiple alignment of many sequences to a common profile HMM.
- hmmconvert： Convert profile formats to/from HMMER3 format.
- hmmemit： Generate (sample) sequences from a profile HMM.
- hmmfetch ： Get a profile HMM by name or accession from an HMM database.
- hmmpress ： Format an HMM database into a binary format for hmmscan.
- hmmstat： Show summary statistics for each profile in an HMM database.
改行コードをパールで変更
- Mac -> UNIX：perl -i -pe 's/\r/\n/g' file
- Windows -> UNIX：perl -i -pe 's/\r\n/\n/g' file
ログアウトした後もバックグランド実行を続ける方法
-
nohup perl test.pl test > out &
-
ログアウトしても仕事をしてくれる
シロイヌナズナのAFFYのマイクロアレイデータの集め方
-
NCBIのTAXOMYの部屋から、Arabidopisi
thalianaのGEOのデータに行く。
-
Platformsに行き、Affymetrix GeneChip Arabidopsis ATH1 Genome Arrayを選び、そのGSEから始まるデータセットをみる。
-
モジュールNet::FTPを使い、データをダウンロードする。
例）
$ftp=Net::FTP->new("ftp.ncbi.nih.gov", Debug=>0);
$ftp->login("anonymous", "-anonymous@");
$ftp->cwd("/pub/geo/DATA/supplementary/series/ディレクトリ名");
または、$ftp->cwd("/pub/geo/DATA/SeriesMatrix/ディレクトリ名");
$ftp->binary;
$ftp->get(ディレクトリ名_RAW.tar);
@files=$ftp->ls() #これをやるとカレントディレクトリ内のファイルを取得でき、これを基に、mget、mputみたいなこともできる。
$ftp->get("filelist.txt");
$ftp->quit;
Rを使ってAFFYのマイクロアレイデータのCELファイルからの解析
-
同じディレクトリにcelファイルを置いとく。
library("affy")
Data<-ReadAffy()
list.celfiles() #確認、読み込んだCel Fileを見せてくれる。
-
一つ一つのチップデータを独立に正規化する。
library(affy)
data <- ReadAffy()
eset_mas <- mas5(data)
write.exprs(eset_mas, file="data_mas.txt")
-
TAIRで、プローブとAGIコードの関係を示すファイルをこのFTP(ftp.arabidopsis.org/home/tair/Microarrays/Affymetrix/)から得る。一つの遺伝子に対して、プローブが複数あるときは、平均をとる。
参考のPERLのスクリプト
ファイル名を"data_mas.gene”として保存
-
SAMで、PvalueとLIMMAでP.valueを推定する統計テスト
例）
library("affy")
library("limma")
library("samr")
t<-read.table("data_mas_gene", header=T, row.names=1);
data.cl1 <- c(rep(1, 4), rep(2, 2))
data1 = list(x=t, y=data.cl1, geneid=rownames(t), genenames=rownames(t),
logged2=TRUE)
samr.obj <- samr(data1, resp.type="Two class unpaired",
nperms=10)
p_value_sam <- samr.pvalues.from.perms(samr.obj$tt, samr.obj$ttstar)
p_value_sam <- p_value_sam[sort(names(p_value_sam))]
data.cl2 <- c(rep(0,
4), rep(1, 2))
design <- model.matrix(~data.cl2)
fit <- lmFit(t, design)
fit <- eBayes(fit)
tmp.out<-topTable(fit, coef=2, number=(length(row.names(t))), adjust="fdr")
EBall<- tmp.out[order(tmp.out[,1]),]
out<-cbind(EBall$ID, p_value_sam, EBall$P.Value, EBall$adj.P.Val)
write.table(out,file="data_mas_out",quote=F,eol='\n',row.names=F,col.names=F,sep='\t')
PerlからRを呼び出す方法。
-
スクリプトを書き、それを呼び出すという方法ができる。
例）
open(SCRIPT, "> tmp.R");
print SCRIPT "library(\"affy\")\n";
print SCRIPT "data<-ReadAffy()\n";
print SCRIPT "eset_mas <- mas5(data)\n";
print SCRIPT "write.exprs(eset_mas, file=\"data_mas.txt\")\n";
close SCRIPT;
system("R --vanilla --slave < tmp.R");
MrBayesの使い方。
-
NEXUS FORMATのアライメントされたファイルを用意
-
"mb"でプログラムスタート、"execute ファイル名”で配列データを読む
-
lset ：根本となるモデル、 prset：事前確率を構築するモデル を設定する。”showmodel”でモデルを確認できる。
-
Nucmodelは、特に情報がなかったら4by4で、コドン配列だとしたら、codonに”lset Nucmodel=codon Rates=Gamma”をタイプして変える。そうでないとしたら、”lset Nstl=2 Rates=Gamma”でいいと思う。Nstはおそらくパラメータの数、個人的には１（JC？)か２（Twoパラ？）でいいと思う。でも、６をTutorialでは薦めているっぽい。
-
Prsetには、6個のタイプのモデルがある（トポロジー, ブランチの長さ、各塩基の頻度、塩基置換のパターン、置換しないサイトの割合、塩基置換のガンマパラメーター）。正直どういう効果がでるかなど不明であるため、Defaultでいいのではないか？
-
mcmcで解析が始まる。"mcmc ngen=10000" 10000回行い、やめるかどうかを決める。理想的には、0.01を切るまでやったほうがいいらしい。
-
sumt burnin=250
Rで図を4分割にし、ヒストグラムを書く方法。
以下のようなデータフレームがあるとする。
V1 V2
1 0.02886309 V
2 0.88217513 V
3 0.48842642 S
4 0.39763334 R
4 0.39763334 O
...............
.........
....
..
1. par(mfrow=c(2,2))で図を４つに分ける。（ちみみにpar(new=T) は 1 つ目の図に上書きする）
2. hist(t[t$V2=="S",]$V1) ＃Sのみのヒストグラムを作る。
3.
hist(t[t$V2=="R",]$V1) ＃Rのみのヒストグラムを作る。
4.
hist(t[t$V2=="V",]$V1) ＃Vのみのヒストグラムを作る。
5. hist(t[t$V2=="O",]$V1) ＃Oのみのヒストグラムを作る。
6. brlist<-seq(0,1000,by=100)
7. hist(t[t$V2=="O",]$V1, breaks=brlist)でX軸を100ずつづらしたものに変えることができる。
6. ちなみに、boxplot(V1~V2,data=t) で、boxplotのグラフを作ることができる。
7. boxplot(V1~V2,data=t,ylim=c(0,2))で、Y軸の範囲を指定できる
図の設定と保存
1. pdf("test.pdf"1, width=7, height=9.9) #A4サイズのPDFファイル名を指定
2. par(ps = 12) #文字サイズの変更
3. par(oma=c(2,2,2,2)) #ページ全体の余白を設定
4. par(mfrow=c(2,2))で図を行、列に各二個づつ。
5. par(mar=c(3,3,2,0))で、各図の余白を指定、下、左、上、右の順番
6.
par(mgp=c(2,1,0))で、各図の説明、ラベル、軸の位置
7.
plot(x, y, main="上部にタイトル",xlab="x軸のタイトル",ylab="y軸のタイトル")
8. dev.off()
False positive, False negative, Sensitivity and Specificity
|
検査がTrue |
検査がFalse |
真にTrue |
A |
C |
真にFalse |
B |
D |
False positive rate= C/(A+C)
False negative rate=B/(B+D)
Sensitivity=A/(A+C)
Specificity=D/(B+D)
つまり、
１−False positive rate = Sensitivity
１−False negative rate = Specificity