東北大学 生物多様性進化分野 河田研究室
河田研の人々
PukiWikiの使い方
おすすめのソフトウェア

詳しい人がどんどん更新していってください。知識は共有しましょう。

関連Software覚書

  • トラッキング
    • idTracker
  • 自然選択検出
    • dN/dS
      • PAML(codeml):二つのモデル(ある枝でdN/dS>1になるモデルと全ての枝が同じdN/dSの値をとるモデル)を比較(尤度計算)するというのが基本的な使い方(branch-siteモデルがよく使われる)。選択を受けた可能性のあるアミノ酸も検出してくれる。branchモデルでは各枝のdN/dSを違う値にすることもできる。概要はこちら。詳しい説明は井上潤さんのページ参照。オプションの設定は以下がオススメ。
        fix_blength = 0 #系統樹ファイルはトポロジーだけ使い、枝長は無視。再計算して求める。 "Branch length is defined as number of nucleotide substitutions per codon (not per neucleotide site)."
        cleandata=1 #Nなど曖昧な塩基を取り除く。
        今後まとめてきます
  • HyPhy:BSRELでPAMLのbranch-siteと同じようなことができるらしい
    PAMLと違って、正の選択が起きた枝の数に制限が無い?
  • Fst, Tajima's D, HKA testとか
    • Arlequin:Windows向けのGUIが基本っぽい。アルレッキーノのフランス名らしい。
    • VCFtools:CUIで使えて便利。vcfファイルから計算できる。Fstはwindow解析もできる。使い方が簡単でとにかく楽。Tajima's Dはwindow stepの指定ができない。のでスクリプト作った。最新版(v0.1.14)はうまくconfigureできないので注意。configure.acをconfigureに変換する必要あり。詳しくはメリーさんまで。Watterson's θは-SNPdensityコマンドでsegregating site数を数えた後、n(=配列数)-1の調和数(n-1番目までの自然数の逆数の和)で割ればいい。
    • DnaSP:windows、GUI。論文でよく見る割にあまり使い方の情報がなくて困っている。誰か助けて。
    • PopGenome?:Rのパッケージ。
      アライメントの非同義置換、同義置換の「数」を数えてくれて、MK testに使える。
      自分で数えるのは、大変だし色々考慮すべき問題があるので、やめましょう。
          GENOME.class <- readData("~/Desktop/Collagen/Sequences/for_Popgenome")

      #指定したディレクトリ以下のファイル全部に対してやってくれる

          GENOME.class <- MKT(GENOME.class,list(1:1,2:2))

      #どっからどこまでが集団かを指定(たぶん)

          get.MKT(GENOME.class)

      #これで、同義・非同義置換の数とαの値が出てくる。一気に。楽勝。

    • PLINK:GWASやQTL解析なんかによく使われるソフトだが、pairwise連鎖不平衡解析やら近交係数やらHW平衡からの逸脱やらなんでもできる優れもの。
    • libsequence:C++の知識を必要とする。
    • MLHKA:maximum likelihood based HKA testを行える。C++のファイルをコンパイルすれば使える。説明も簡単。複数座位を使うので集団動態の影響は受けづらいらしいが、それでも集団の拡大やボトルネックでp値が有意になることがあるらしい。
    • HKA:こちらは集団動態のモデルも組み込むので、MLHKAよりもさらに集団動態の影響を受けづらい。
  • 有害変異検出
    • Provean:配列の保存度から推定。値が-2.5以下が有害変異ということだが、変異前の配列と変異後の配列を逆に入力すると、絶対値がほぼ同じで+-が逆の値が出てくる。どちらの変化も同じと考えれば+2.5以上も有害変異(?)。CUIでできるはずだが牧野先生しか成功していない。だれか助けて。やっとできたー。
    • Polyphen2:たぶんProveanと同じ原理。ヒトの多型データが対象?EnsemblのVEP(?)ファイルには既に結果が盛り込まれてた気がする。自分でデータベースを構築することで、様々な種に使えるらしい。
    • SIFT:Polyphenと共に、この手の解析で最もよく見るソフトウェア。スコアは0~1の値をとり、0.05以下が有害。sequence homologyを基に計算すると書いてあるが、polyphenやproveanとどこが違うのか、よう分からん。
    • MutationTaster?:"disease causing automatic" or "disease causing"を示したら有害変異。
    • GERP++:値が5以上だと有害変異。
    • PhyloP:値が3以上だと有害変異。
    • CADD:ソフトウェアというより計算済みのスコア。配列の保存性だけでなく、病原性やら調節効果も包括しているとか。ヒトだけ。
  • ハプロタイプ解析
    • Haploview:LDを逆三角形の図で可視化できる。しかしあまりにサイト数が多いとsvgで保存できなくなる。haplotypeも見れる。いずれにせよ、画像がpdfで保存できないため、その後の図の編集がしづらい。
    • PopART:最近注目されている新しいソフト。綺麗な図が描ける。こんな感じのインプットファイルを用意する。
      #NEXUS
      Begin data;
      Dimensions ntax=10 nchar=5;
      Format datatype=dna missing=? gap=-;
      Matrix
      
      E       00000
      H       11100
      A       01001
      I       01011
      D       11001
      G       11111
      C       11011
      F       00011
      B       11101
      J       01000
      ;
      End;
      
      BEGIN TRAITS;
      Dimensions NTRAITS=4;
      Format labels=yes missing=? separator=Comma;
      TraitLabels Africa Europe China Japan;
      Matrix
      
      E       13,50,54,59
      H       0,0,1,0
      A       115,36,78,58
      I       0,0,0,2
      D       14,2,0,0
      G       0,0,14,40
      C       22,72,36,35
      F       0,1,0,0
      B       52,37,23,13
      J       0,0,0,1
      ;
      End;
      haplotypes.png
    • NETWORK:windows専用。よく使われている気がする。インプットがrdfというやつでめんどくさい。お絵描き専用?のPopARTとは違い、コアレセントによる解析なんかもできるらしい。
    • TCS:簡素。インプットはNEXUSかPHYLIP。図がしょぼい。
  • コアレセント
    • 集団動態推定
      • dadi:集団の動態を推定する。
      • ms:推定された集団動態のパラメーターを代入して、中立を仮定した場合の多型情報を作り出す。
      • msms:自然選択もパラメーターとして入れられる。
    • 起源年代推定
      • DMLE+:
      • GENETREE:ハプロタイプの起源年代(TMRCA)を推定できる。
  • マッピング
    • tophat2:RNA-seqのreadをマッピングする際、選択的スプライシングやイントロンによって張り付かないreadがあるが、その問題を解決してくれる。readを分けて貼り付けてくれる。
    • bowtie2:ゲノムにゲノムのreadをマッピングするならこれ?
    • samtools, bcftools:2016年9月現在、スパコンのsamtoolsとbcftoolsのバージョンが合致していない。bcftoolsのバージョンが古いままなので、最新版をダウンロードする必要がある。
    • BWA:bowtie2よりいいらしいが、時間がかかる。
  • 系統樹作成
    • RAxML:MLで作成してくれる。のに速い。15種、250万アミノ酸配列を1時間半くらいで仕上げる(ML serachの場合)。boot strapを計算するとけっこう時間かかるかも。
      RAxML PTHREADS:マルチスレッド版。ML探索に有効。
      RAxML MPI:マルチプロセス版。bootstrapを並列できる。
      RAxML HYBRID:上2つを混合。
    • ExaML:RAxMLの大規模解析版みたいな。大量データをRAxMLよりも速く計算してくれる。ただbootstrapを出せないので微妙?
    • PhyML:メールすればソースコードをもらえる。正確だってよ。
    • FastTree:とにかく速いとのこと。
    • FigTree?:作った系統樹ファイルを画像で表示してくれるソフト。一番上の部分が欠けるのはどうにかならないのか。
    • iTOL:ブラウザ上から系統樹の見た目を編集するツール。様々なオプションがあり、結構きれいに作れます。無料。
  • アライメント
    Ziheng Yangら(PAMLの作者)によると、PRANK > MUSCLE & MAFFT > ClustalW の順に正確(dN/dSなどを計算する時のみ?)。
    • ClustalW2,ClustalX:WはCUI、XはGUI。
    • PRANK:codonモデルがオススメ。実際、Clustalよりも正しいと思われる結果が返ってくる。こちらからソースコードをダウンロード。乱数を指定しないと毎回結果が微妙に変わる。-seedオプションで乱数を指定すべし。
    • MUSCLE:
    • MAFFT:個人的にはPRANKの次点扱いだが、とにかく速い。mafft-linsiが正確。--auto オプションを使うと配列数に対して最適なオプションを勝手に選んでくれる(情報科学研究科山田さんより)
    • SATe:反復性があるとか。
    • PASTA:SATeを改良したやつらしい。
    • T-Coffee:
    • trimAl:ギャップのあるサイトを取り除く。NとかXは残る。
  • 相同性検索
    • blast+:オプションはblastとは異なるので注意。"-num_threads 数字"でマルチスレッドでのCPUを増やせる?(def_slotを忘れずに)
  • アセンブル
    • Trinity:RNA-Seqのde novoアセンブル
  • 変異のコール(RADとかMIGとか)
    • Stacks
    • pyRAD:Stacksと違ってインデルを考慮してくれる。アウトプットは.nex/.phy/.str/.vcfなど。pythonで書かれているためインストールがちょっとめんどい。de novo専用。
    • GATK:Bayesian genotype likelihood modelというのを使ってSNPsやインデルのコールをしてくれる&そのために必要なインデル領域を探したり、リアライメントしたり、クロリティスコアの再計算をしたり、いろんなことをしてくれる。ゲノムサイズとデータサイズ次第ではあるが、時間はかかる。
  • 遺伝子予測
    • AUGUSTUS : 真核生物用の遺伝子予測ソフト。Linuxbrewで入れましょう。
  • RNAエディティングサイトの検出
    • RES-Scanner
  • 分子動力学計算
    • gromacs:分子動力学法によって分子の挙動をシミュレーションするもの。アミノ酸の違いによって、タンパク質の働きがどう変わるかを見ることができる(はず)。
  • その他便利なもの
    • ファイル変換とか
      • PGDSpider:Arlequin,SAM,BAM,FASTA,STRUCTURE,VCF(もっとある)など色々なファイルの形式を相互に変換してくれる。便利。
      • EMBOSS:色々なパッケージが入ってる。あまり把握してない。とりあえずヌクレオチド→アミノ酸に変換してくれるtranseqが便利。不完全なコドン(塩基が2つとかNが入ってるとか)でも、アミノ酸が確定する場合は何事も無く翻訳されるので注意。
    • newick形式の編集
    • ダウンロード
      • lpm:本来は所有者権限でないとインストールできないものをインストールする。遺伝研のスパコンに新たなソフトを入れたい時とかに。local権限でインストールできるようになったので、もういらないっぽい。
    • ターミナルを使いやすく
      • cdto:Finderのタブに追加しておけば、好きな所でターミナルを開ける。めっちゃ便利。こちらに解説。

インストール関連覚書

インストールする際に困ったこととか。

ProveanのCUI化(2016/8/19更新)

ProveanはSIFTやPolyphenと同じく、アミノ酸置換のマトリックスと配列の保存度を用いて、有害変異を推定するソフト。
今までCUI化に挑戦しては失敗し続けてきたが、今度こそ成功した気がする。
まずソースコードをダウンロード。http://provean.jcvi.org/downloads.phpから。
次にCDHITというソフトをインストール。解凍してmakeする。
nrデータベース(SIFTのuniref90みたいな全種の全配列データ,2016年8月10日時点でnr00~53まである)をダウンロード、解凍する。量が多くてちょっとめんどくさい。
proveanをダウンロード後、解凍したらインストールしていく。ここでCDHITとnrデータベースのパスを指定する。
nrのディレクトリ/nr←nrという名前まで指定しておくことが重要。
Blastが既にインストール済みでパスを通している場合は、PSIBLASTとBLASTDBCMDは勝手にパス指定される(?)。

./configure CDHIT=**/cdhit-master/cd-hit BLAST_DB=**/nr/nr --prefix=${HOME}/local
make
make install #色々エラーっぽいの出てきても気にしない。パス指定せいと言われたら従う。

/scriptsディレクトリに行って、実行ファイルであるprovean.shを少し編集。
スクリプトの下の方にある以下の部分を変える。

COMMAND="$SCRIPT_DIR/provean #$SCRIPT_DIRのままだとうまくない。
COMMAND=“**/provean-1.1.5/src/provean #フルパス指定してやる。

あとはアミノ酸の配列ファイル(fasta形式,input1)と置換リストのファイル(SIFTと同じ,input3)を用意して、

provean.sh -q input1 -v input3

で実行。しばし待たれい。実行時間はアミノ酸配列の長さに比例するらしい。

save_supporting_setのオプションでアウトプットファイル名を指定する。

provean.sh -q input1 -v input3 --save_supporting_set output

ただProveanの場合、結果はprintされるので、

provean.sh -q input1 -v input3 > output

としたほうがいい。

SIFT(ver.6.0.0)のCUI化(2016/8/19更新)

SIFTはヒト以外の生き物にも使うことが出来る汎用性の高いソフト。
ちょっと頑張ったらスパコン内で動かせるようになったので報告。
ちなみにSIFTは様々なバージョンがあり、開発者が異なっていたりする。
今のところsift4.0.3b/jcvi-sift-1.0.3/sift6.0.0と3つあるが、ここではsift6.0.0について扱う。
まず、http://sift.bii.a-star.edu.sg/から、SIFTの最新版(sift6.0.0)ソースコードをダウンロードする。
(8/19日追記:sift6.0.0にエラーがあったので報告したら6.0.1を作ってくれた。そちらを使うべし。)

Blastの最新版(2.4.0+)をftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/からインストール。

curl -O ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.4.0/ncbi-blast-2.4.0+-x64-linux.tar.gz

Uniprotから全種の全アミノ酸配列をダウンロードする。

curl -O ftp://ftp.uniprot.org/pub/databases/uniprot/uniref/uniref90/uniref90.fasta.gz

その後fasta化。

zcat uniref90.gz | perl -pe 's/^>UR090:UniRef90_/>/'  > uniref90.fa

makeblastdbでデータベース化する。30分くらいかかる。

makeblastdb -in uniref90.fa -dbtype prot

sift6.0.0/binにあるSIFT_for_submitting_fasta_seq.cshとSIFT_for_submitting_NCBI_gi_id.cshを編集。

setenv NCBI **/ncbi-blast-2.4.0+/bin
setenv SIFT_DIR **/sift6.0.0
set tmpdir **/** #アウトプットを置くところ。

後は調べたい配列と置換リストのファイルを作って、sift6.0.0/bin/に移動し、以下のコマンドを打つ。

csh SIFT_for_submitting_fasta_seq.csh input1 uniref90.fa input3

input1は見たいアミノ酸配列(fasta形式)。input3はアミノ酸の置換リスト。↓みたいな感じ。

S76T
V189I

これで、さっき指定したディレクトリに結果が出るはず。ただ20分くらいかかるかもしれない。
ちなみにinput3の部分を-とすると、全てのアミノ酸変化を全パターン調べて結果を出してくれるらしい。

たまにblastの時間がかかりすぎる(30分以上)と”CPU時間制限オーバー”などと出力され、強制終了することがある。
スパコンの問題かと思い、遺伝研に問い合わせたところ、そうではないっぽい。
blastの計算時間に30分という時間制限が課されていたらしい。
/sift6.0.0/bin/seqs_chosen_via_median_info.cshの以下の部分を変更すればいい。

limit cputime 30m #60行目。30mは30分の意味。これを10h(10時間)とか、unlimited(無制限)とかにすればOK。

遺伝研スパコンのRにパッケージを追加する

遺伝研スパコンには2つのバージョンのRがインストールされている(2.14.1と3.1.1)。(遺伝研スパコンプログラミング環境
デフォルトでは古い方にパスが通っているが、パスを書き換えれば新しい方(といっても2014年だけど)を使える。

$ export PATH=/usr/local/pkg/R/3.1.1_icc_mkl_shared/bin:$PATH


遺伝研スパコンのRにinstall.packages()でパッケージを追加しようとすると、既存のディレクトリに書き込めないと言われるが、勧められるがままにホームディレクトリ直下にディレクトリをつくりミラーサイトを選択すれば、そこにパッケージをインストールできる。

'lib = "/usr/local/pkg/R/3.1.1_icc_mkl_shared/lib64/R/library"' は書き込み可能ではありません
Would you like to use a personal library instead?  (y/n) y
Would you like to create a personal library
~/R/x86_64-unknown-linux-gnu-library/3.1 to install packages into?  (y/n) y
--- このセッションで使うために、CRAN のミラーサイトを選んでください --- 

もちろん元々入っていたパッケージと新しく入れたパッケージとで場所が異なることになるが、Rの.libPaths()関数でパッケージライブラリのパスを調べると

.libPaths()
[1] "/lustre3/home/username/R/x86_64-unknown-linux-gnu-library/3.1"
[2] "/usr/local/pkg/R/3.1.1_icc_mkl_shared/lib64/R/library" 

とちゃんと2つ出る。library()で読み込める。

おまけ(https://qiita.com/uri/items/8fa7e194013977899f86
パッケージのインストール先(ディレクトリ)を知りたいときにはfind.package, path.packageの2つの関数を使うと良い。

find.package("dplyr")
path.package("ggplot2")

この2つの関数の違いは、すでにパッケージが読み込まれているかどうかで挙動が変わってくる。目的のパッケージを読み込んだ後であればpath.packageで良いが、読み込む前ならばfind.packageを使用しないとパッケージは見つからない。

遺伝研のスパコン内で用いるPerlスクリプトにRを組み込みたい

今まで、スパコン内ではPerlスクリプト内のRを使えないと思っていたが、ついに成功したのでここに報告。
Statistics::RというCPANのモジュールを用いてPerlとRをつなげる。
もしかしたら必要ない部分あるかも。各自で判断してけろ。

まずスパコンのホームディレクトリでCPANをインストール。

yum install perl-CPAN.x86_64

CPANを開く。質問が出てきたら[yes]→[manual]→[yes]を選択。

cpan

次にStatistics::Rをダウンロードする。
http://search.cpan.org/~fangly/Statistics-R-0.34/lib/Statistics/R.pmから。
→.tarファイルをスパコンに転送(スパコン内でcurlでダウンロードしたら行方不明になった)。

送ったtarファイルを解凍。

tar xvf /home/daikisato177/Downloads/Statistics-R-0.34.tar

Statistics::Rのあるディレクトリに移動。

cd /home/daikisato177/Downloads/Statistics-R-0.34/

makeしていく。ルート権限のない環境でやるためには下のPREFIX=~の部分が肝心らしい。

perl Makefile.PL PREFIX=~/lib/perl5
make
make test
make install

あとはパスを通す(?)作業。下のいずれかを書き加える。.bash_profileをいじる方が楽なのでオススメ。
1) .bash_profileに以下を加える。

export PERL5LIB=/home/daikisato177/Downloads/Statistics-R-0.34/lib

2) Perlスクリプトの先頭に以下を加える。

use lib '/home/daikisato177/Downloads/Statistics-R-0.34/lib’; #各自のディレクトリ名で。

最後に、使いたいPerlスクリプトの先頭に下の一文を加えるのを忘れずに。

use Statistics::R; 


これで使えるようになっているはず。
以下、Perlスクリプト内でRを使う際のコード例。

#!/usr/bin/perl

use strict;
use warnings;
use lib '/home/daikisato177/Downloads/Statistics-R-0.34/lib';
use Statistics::R;

my $R = Statistics::R->new();
$R->startR;
$R->send(qq` #ここから下がR内の処理。Perlの記号も通用する(\nとか、作った変数とか)。
a <- 30\n
a
`) ; #ここまで。
my $res = $R->read;
$R->stopR();
print "$res\n”;

すると

[1] 30

と出力される。

pyRADを遺伝研のスパコンにインストールする

ちょっと厄介だったのでまとめておく(文責:稲田)。
みっちーさん方式に則ってpythonの仮想環境を立ちあげて使います。
遺伝研にログインしてから、以下のコマンドを順番にやればうまくいくはず、たぶん。
ホームディレクトリに以下の3つのディレクトリをつくることになるので、気に食わない人は適宜改変してください。
・もろもろのインストール先の「local」ディレクトリ
・python2.7の仮想環境用の「py27」ディレクトリ
・ダウンロードしたものを入れておく「software」ディレクトリ

下ごしらえ。

cd ~/
mkdir local
mkdir software

仮想環境を整える。

virtualenv -p python2.7 --no-site-packages py27
source ~/py27/bin/activate

(py27)と表示されるようになるはず。
これが仮想環境が立ち上がっている証拠らしい。

pip install numpy
pip install scipy

これで完了。一旦、仮想環境を閉じても良い。

deactivate

次に、pyRADに必要なもの(計4つ)をどんどんインストールしていく。
まず、1つ目。

cd ~/software
curl -O http://ftp.gnu.org/gnu/autoconf/autoconf-latest.tar.gz
tar xvfz autoconf-latest.tar.gz
rm xvfz autoconf-latest.tar.gz
cd autoconf-2.69/                        ##versionは2016/08現在のもの
./configure --prefix=/home/○○/local/   ##○○はユーザー名
make
make install

2つ目に行く前にパスを通す

nano ~/.bash_profile

以下の一文を追加

PATH="/home/○○/local/bin":$PATH

ちなみに「PATH=$PATH:"/home/◯◯/local/bin"」としてしまうとうまくいかないので注意。
  ※実はautoconfはすでにスパコンに入っている。ただversionが古くて使えない。
   今インストールしたautoconfに、この古いversionのautoconfよりも優先度の高いパスを通さないといけない。
   よって、"/home/◯◯/local/bin"を前に追加しないといけない。
次に、2つ目。

cd ~/software
curl -O http://ftp.gnu.org/gnu/automake/automake-1.15.tar.gz       ##versionは2016/08現在のもの
tar xvfz automake-1.15.tar.gz
rm automake-1.15.tar.gz
cd automake-1.15/
./configure --prefix=/home/○○/local/
make
make install

次に、3つ目。

cd ~/software
git clone https://github.com/torognes/vsearch
cd vsearch
./autogen.sh
./configure --prefix=/home/○○/local/
make
make install

最後に、4つ目。

cd ~/software
curl -O http://www.drive5.com/muscle/downloads3.8.31/muscle3.8.31_i86linux64.tar.gz     ##versionは2016/08現在のもの
tar xvfz muscle3.8.31_i86linux64.tar.gz
rm muscle3.8.31_i86linux64.tar.gz
mv muscle3.8.31_i86linux64 muscle
mv muscle ~/local/bin/

準備完了。
あとはpyRAD本体をインストールして終了。

cd ~/software
git clone https://github.com/dereneaton/pyrad.git

使用する時は常に仮想環境を立ち上げる。

source ~/py27/bin/activate                   ##仮想環境の立ちあげ
python ~/software/pyrad/pyrad/pyRAD.py -h    ##pyRADのヘルプを表示
deactivate                                   ##仮想環境の終了

qsubでジョブを投げて使うときはjobの中で立ちあげるべし。

コマンド&スクリプト覚書

極めて個人的なメモ。

General command

  • root権限とインストール
    su #linuxでroot権限でインストールする時。
    sudo su #macの場合。
    perl Makefile.PL PREFIX= #スパコン内では、パス指定によりlocal権限でインストールする。
    ./configure --prefix=/~ #こんな感じもあり。いずれにせよmakeする前にやる。
  • tar系解凍コマンド
    tar zxvf ~.tar.gz / ~.tgz
    tar xvf ~.tar
    gunzip ~.gz
    tar jxf ~.bz2
    tar -czvf dir.tar.gz dir #dirという名前のディレクトリをtar.gz形式に圧縮
  • cd
    cd - #一個前に作業していたディレクトリに戻る。
  • cp
    cp -a ~ #パーミッション状態なんかも含めてコピー。
    cp -r ~ #ディレクトリのコピー。
  • rm
    rm -r ~ #ディレクトリの削除。
    rm -f ~ #ファイルの強制削除。
    rm -rf ~ #ディレクトリの強制削除。
  • less
    Space #次のページに飛ぶ。
    b #前のページに飛ぶ。
    g #最初の行に飛ぶ。
    G #最後の行に飛ぶ。
    (行数)< #指定行に飛ぶ。
    /(文字) #ファイルの後方に向かって検索。
    ?(文字)#ファイルの先頭方向に向かって検索する。
    n #再度検索。
    N #先頭方向に再度検索。
    less -N #左側に行数表示。
  • du
    現在使用している容量を表示してくれる。
    du -m #メガバイト表示。
    du -s #ディレクトリの指定。
  • ls,wc
    ls #現在のディレクトリにあるファイルを全て表示。
    ls -1 #一列で表示。
    wc ~ #指定したファイルの行数、文字数を表示。
    ls -1 | wc -l #現在のディレクトリにあるファイルの個数を表示。
  • head,tail,sed
    head -n 100 input > output  #上から100行だけ取り出す。
    tail -n 100 input > output  #下から100行を取り出す。
    tail -n +2 input > output       #2行目以降を取り出す。ヘッダーを除きたい時に便利。
    head -n 100 input | tail -n 1 > output  #100行目だけを取り出す。
    sed -e '1d' input > output  #1行目だけ削除したファイルを作る。
    sed -e '$d' input > output  #最終行だけ削除したファイルを作る。
  • grep
    grep aaa input > output     #aaaという文字列を含む行だけ抜き出す
    grep -v bbb input > output    #bbbという行を含まない行を抜き出す
    grep aaa input | grep -v bbb > output #aaaを含み、かつbbbを含まない行を抜き出す
  • while
    条件が成り立つ限りdoからdoneまでの動作を繰り返すコマンド。perlと同じように、ファイルを一行ずつ読み込ませることができる。
    while read name #nameは変数名でなんでもいい。ここでは$はいらない。読み込んだ行は$nameに格納される。
    do #処理の開始。
    ~ ~ ~ $name #適当にコマンドを。
    ~ ~ ~ $name
    done < file.txt #処理の終了。file.txtは読み込むファイル。
  • chmod
    ファイルやディレクトリの権限を変更するコマンド。
    ソフトの実行ファイルはあるのに起動しない、みたいなときには以下のようにする。
    chmod u+x ファイル名

peco

  • pecoってなに?
    インクリメンタルサーチツール。日本語で言うと、逐語検索ツール。
    逐語検索とは、検索したい単語をすべて入力した上で検索するのではなく、入力のたびに即座に候補を表示させること。
    なんだ、ただの検索ツールか、と思うことなかれ。
    こいつのすごいところは、検索結果の一覧から矢印キーで選択できるところにある。
    pecoの機能自体はシンプル故にいろいろ応用も効く。
    常日頃、cdやファイル選択などにストレスを抱えているあなたも、そうでもないあなたも、一度設定したら使わない日はない、そんなツールです。
    類似品にfzfとかいうのもあるらしい。
  • インストール
    Homebrewを導入してる人は
    brew install peco
    でいけるはず。
    してない人はこのへん↓みてください。
    https://github.com/peco/peco
  • 基本動作
    pecoは単独では何の役にも立たない。
    他のコマンドとパイプでつなぐことで、はじめて効力を発揮する、そんなやつです。
    基本動作は、
    パイプで渡されたリストを表示して、ユーザーに選ばせて、選ばれた項目を出力する
    といった感じのコマンドです。
    echo -e "aaa\nbbb\nccc" | peco
    ls -F | peco
    例えばファイル一覧をリストにして、出力結果を変数に入れてlessすれば、選択したファイルの内容を見るコマンドができる。
    less "$(ls -F | grep -v / | peco )"
    ### ls の -F はディレクトリに / を付けるオプション ###
    ### grep -v でリストから / の付く文字列(=ディレクトリ)を除いている ###
    さらにこれにショートカットを作ったりすれば…
    ### .bashrcや.bash_profileに以下の文を記述 ###
    alias less-peco='less "$(ls -F | grep -v / | peco )"'
    
    ### もしくは、以下のようにしてもよい ###
    function less-peco(){
     less "$(ls -F | grep -v / | peco )
    }
    
    ### 変更後は再起動するか、sourceで読み込む必要がある ###
    source .bashrc         # .bashrcに書いた人
    source .bash_profile   # .bash_profileに書いた人
    less-pecoを使えるようになる。
    less-peco
    こんな感じで、いろいろ応用していくことが可能です。
    if条件分岐とかでこちょこちょやれば、好きな仕様のコマンドができます。
  • 応用例
    いくつか応用例を載せておきます。bashならスパコンでも動くはずです。
    拙いところは目をつぶって下さい。
    • ディレクトリを選択して移動:sd
      ### .bashrcに記述 ###
      function peco-cd {
       local sw="1"
       while [ "$sw" != "0" ]
        do
      		if [ "$sw" = "1" ];then
      			local list=$(echo -e "---$PWD\n../\n$( ls -F | grep / )\n---Show hidden directory\n---Show files, $(echo $(ls -F | grep -v / ))\n---HOME DIRECTORY")
      		elif [ "$sw" = "2" ];then
      			local list=$(echo -e "---$PWD\n$( ls -a -F | grep / | sed 1d )\n---Hide hidden directory\n---Show files, $(echo $(ls -F | grep -v / ))\n---HOME DIRECTORY")
      		else
      			local list=$(echo -e "---BACK\n$( ls -F | grep -v / )")
      		fi
      		
      		local slct=$(echo -e "$list" | peco )
      		
      		if [ "$slct" = "---$PWD" ];then
      			local sw="0"
      		elif [ "$slct" = "---Hide hidden directory" ];then
      			local sw="1"
      		elif [ "$slct" = "---Show hidden directory" ];then
      			local sw="2"
      		elif [ "$slct" = "---Show files, $(echo $(ls -F | grep -v / ))" ];then
      			local sw=$(($sw+2))
      		elif [ "$slct" = "---HOME DIRECTORY" ];then
      			cd "$HOME"
      		elif [[ "$slct" =~ / ]];then
      			cd "$slct"
      		elif [ "$slct" = "" ];then
      			:
      		else
      			local sw=$(($sw-2))
      		fi
        done
      }
      alias sd='peco-cd'
      カレントディレクトリの中のディレクトリをリスト → pecoで選択 →
      そのディレクトリに移動してさらにディレクトリをリスト → pecoで選択 →
      …ループ… → 一番上を選択したら終了

      という感じで動作するはず。
      一番上にカレントディレクトリのフルパスを表示したり、
      一番下にホームディレクトリに移動する項目とか、
      オプションで隠しディレクトリも含めて表示する項目とか、
      さらにオプションでカレントディレクトリのファイルを表示する項目とかを実装した。
      ちょっと重いかも。もっと軽いのとかつくれないかな。
      あと、オプションを選択肢の中に放り込んじゃったから、検索にも引っかかってしまう欠点がある。
      検索避けみたいなことができたらいいんだけど。
      実行は以下の通り、.bashrc編集後、sourceで読み込むのをお忘れなく。
      sd
      文字で書いても、イメージしにくいと思うので、とにかく一回お試しあれ。
      我ながらなかなか良いものを作ったと自負してます。
      不具合報告などは稲田まで。
  • ファイル選択
    ### .bashrcに記述 ###
    function peco-〇〇 {
       local dir="$PWD"
       local sw="1"
       while [ "$sw" != "0" ]
    	 do
    	     if [ "$sw" = "1" ];then
    		     local list=$(echo -e "---CANCEL\n$( ls -a -F | grep -v / )\n---HOME DIRECTORY\n---CHANGE DIRECTORY")
    	     elif [ "$sw" = "2" ];then
    		     local list=$(echo -e "---$PWD\n---HOME DIRECTORY\n---Show hidden directory\n../\n$( ls -F | grep / )\n\n$( ls -a -F | grep -v / )")
    	     elif [ "$sw" = "3" ];then
    		     local list=$(echo -e "---$PWD\n---HOME DIRECTORY\n---Hide hidden directory\n$( ls -a -F | grep / | sed 1d )\n\n$( ls -a -F | grep -v / )")
    	     fi
    	 
    	     local slct=$(echo -e "$list" | peco )
    	     
    	     if [ "$slct" = "---CANCEL" ];then
    		     local sw="0"
    	     elif [ "$slct" = "---$PWD" ]||[ "$slct" = "" ];then
    		     local sw="1"
    	     elif [ "$slct" = "---CHANGE DIRECTORY" ]||[ "$slct" = "---Hide hidden directory" ];then
    		     local sw="2"
    	     elif [ "$slct" = "---Show hidden directory" ];then
    		     local sw="3"
    	     elif [ "$slct" = "---HOME DIRECTORY" ];then
    		     cd "$HOME"
    	     elif [[ "$slct" =~ / ]];then
    		     cd "$slct"
    	     else
    		     local sw="0"
    		     ●● "$slct"
    	     fi
    	 done
    	 cd "$dir"
    }
    alias □□='peco-〇〇'
    カレントディレクトリのファイルを一覧にして、pecoで選択。
    必要に応じて、別のディレクトリにもアクセスすることも可能、というコマンド。
    ●●は引数に1つのファイル(ディレクトリは不可)を要求するコマンドなら何でもよい。たぶん。
    less、cat、head、tail、nanoなどなど。
    head -n 100 "$slct"
    みたいにオプションをつけても問題ない。
    ○○は、特殊文字を含まなければなんでも良い、lessで作るなら○○もlessとかにしとけばいい。
    □□にはお好きなショートカットを。
    lessなら、slctlessとか、sleとか、pecoleとか、pleとか。
    短くてわかりやすいのがオススメ。

    これもオプションとして色々つけてみた。
    まず、一番上にCANCELという項目を設けて、実行をキャンセルするオプションをつけた。
    あと、一番下にディレクトリ変更という項目を作った。
    通常ではファイルのみ表示されるが、この項目を選択するとディレクトリも表示されるようになる。
    そこから選択して、一時的に別のディレクトリにアクセスすることができる。
    ディレクトリ変更は上のsdと同様に、ループ処理が施されているため、何度でも可能。
    ただし、sdと違ってディレクトリは移動しない。
    つまり、今のディレクトリにいながら別のディレクトリのファイルを参照できるという仕様です。
    これまた我ながらなかなか便利。
    ただ、sdと同じくオプションも検索対象に含まれてたり、動作が重かったりするのはご愛嬌。
    お試しあれ。

スパコン関連

  • qsub:ジョブ投入
    qsub -l month -l gpu job1 #ジョブを投入するノードを指定。デフォルトはweekノード。
    qsub -l s_vmem=16G -l mem_req=16G  job1 #メモリを16GBに指定。デフォルトは4GB。
    qsub -pe def_slot 4 job1 #スロットを増やす。使用するソフトが対応していないと意味がない。デフォルトは1。
    qsub -N job1 job2 #job1が終わったらjob2を投入。
  • quota;自分が現在使用している容量を確認。
    lfs quota -u ユーザー名 ~/
  • qreport 終了したジョブの詳細を表示。
    qreport -j ジョブID
  • ジョブを一気に終了させたい
    qstat > temp
    perl -lane 'print "qdel $F[0]"' temp > temp2 #temp2を適宜編集してからsh temp2でOK。
  • monthでもweekでもどこでもいいから、空いてるところにジョブを投げたい
    qsub -soft -l month,gpu,short #month-phi/month-gpu/short/weekノードのどれかに投入してくれる
  • phase1とphase2間のやり取り
    phase2→phase1では/home_ph1/ユーザー名/~でファイルを見たり移動が出来るが、phase1→phase2は基本できない。
    ただphase1のmonth fatノードに投入するジョブに限っては/home_ph2/ユーザー名/~が使用可能。
    詳しくは遺伝研の『Phase2システム 計算機利用方法』参照。

  • 大容量ファイルの転送
    遺伝研では、大容量ファイルの転送にはAsperaの使用が推奨されてる。体感、速い気がしないでもない。
    参考URL
    インストール後にここに飛んでGUI操作するか、コマンドライン上で操作する。
    Asperaはphase2にしか繋がっていないので、phase1へのシンボリックリンクを作成すべし。

  • 遺伝研に最新版Rをインストール
    Linuxbrewで簡単に入れられる。
    brew install homebrew/science/r

Perlスクリプト覚書

  • forよりもforeachの方が回転が微妙に速い。
  • push @array,$line;の方が@array=(@array,$line);とするよりも格段に速い。これはマジ。
  • 配列を文字列にする
    join("",@array); #1文字ずつ入れた配列を、再度出力したいときとか。
  • 配列の最後の添字を得る
    $#array;
  • 配列の要素数を得る。以下のどれでも可。
    @array;
    $#array+1;
    scalar(@array);
  • perlスクリプトを終了させる。
    exit;
  • 繰り返し(for,foreach,whileなど)を終了する。
    last;
  • 10進数を2進数に変換。
    $dec = sprintf("%0"."$bin"."b",$a); #$binのところは桁数指定。$aが変換したい数字。
  • 小数点以下の桁を指定する
    $b=sprintf("%.1f", 0.1234) #$b=0.1
  • 1行perl
    perl -オプション 'スクリプト' インプット > アウトプット
    インプットはデフォルト変数に格納される。
    実は1行じゃなくてもよかったり(そん代わり、スクリプトの修正ができない)。
    • オプション
      色々と手間が省けて便利。
      基本、-laneってしとけば問題ない。
      • e:1行perlには必須。
      • n:インプットを一行ずつ読み込んでデフォルト変数に格納。
      • p:インプットを一行ずつ読み込んでデフォルト変数に格納。かつ、各行の最後に"print $_;"
      • l:print時の\nを自動でしてくれる。
      • a:デフォルト変数をスペース区切りで@Fに格納。区切る文字列は-Fで指定できる。
    • 1行perlでできること。
      • Excelファイルをtab区切りテキストにした時に改行がうまくいかないのを直す。
        perl -pe 's/\r/\n/g' input > output
  • 変数名に変数を用いる
    簡単にやるには、以下のように
    ${$name};
    @{$name};
    %{$name};
    しかし、これらはシンボリックリファレンスってのに触れる危険な行為なので、use strictを使うとできなくなる。

    安全にやるには、ハッシュのキーに名前を、値にリファレンスを入れる。
    他のやり方があったら追記してください。

Perlモジュール

  • CPANではなくCPANMを使うとlocal権限でインストール可能。他にも色々できるらしい。こちら参考。
  • モジュールがインストールされているかの確認
    find `perl -e 'print "@INC";'` -name "*.pm" -print | grep モジュール名
  • モジュールを読み込むパスの確認
    perl -e 'print join("\n",@INC)."\n";'
  • モジュールのパスを追加
     .bash_profileに以下を追加。
    export PERL5LIB=パス
     この方法だとcronなどで混乱が起きるとかなんとか。よーわからん。

  • おすすめモジュール
    • Kyoto Cabinet
      ハッシュのメモリ消費量を節約できる。ただし時間はかかるっぽい。参考ページ
      ほかにも、valueで検索できるとか。
      モジュールインストールの際には、環境変数LD_LIBRARY_PATHにKyoto Cabinetのlibディレクトリを追加する必要がある。

Perlの注意事項

  • use strict;
    これをつけるとmyで宣言していない変数は使えなくなるので、スペルミスが防げる。
    危険な方法も使用禁止になる(シンボリックリファレンスとか)。
    この1文があるせいで動作しないような場合もあるから、時には省くのもアリっちゃアリ。
    けど、本来はこれをクリアするようなスクリプトを書くべき。
  • foreachで配列の要素を格納した変数を変更すると、元の配列も変化する。
    my @array = (1,2,3);
    foreach my $item(@array){
    $item++; #配列の要素を格納した変数を変更
    }
    print join("\s",@array) . "\n";
    
    ↓
    
    2 3 4
  • 正規表現
    • 優先順位
      マッチする文字列が複数箇所ある場合、以下の優先順位に則って選ばれる。
      1. 文字列がマッチするもの。
      2. より左側にあるもの(マッチする文字列の左端で判断する)。
      3. より長くマッチするもの。

    • 置換で文字列を全て削除すると、undefではなく空文字列になる。
      空文字列は長さが0なので、条件式にはlengthが有効。(definedでは真が返ってくる)
      参考サイト
  • 小数点以下の計算を行うと、バグる時がある。いまいち理由がわからない。

Java

遺伝研に標準インストールされているJavaが以下のようにバグっていたので一からインストール

Error occurred during initialization of VM
Could not allocate metaspace: 1073741824 bytes

ここからLinux x86のをダウンロードし解凍するだけ。コンパイルは必要なし。
ホームで作業しすぎると上記のエラーが出るっぽい。ログインし直せば問題なし。

次に環境変数のPATHとCLASSPATHをいじる

export PATH=Javaの絶対パス/bin:$PATH
export CLASSPATH=.:Javaの絶対パス/lib:$CLASSPATH ###カレントディレクトリ'.'に注意

これで使えるようになるハズ

自作スクリプト

プライマー設計

注意点:EnsemblのcDNAデータに対しBlastする前提で書いてます。
(大)問題点:Blastは14bp以下のヒットを見つけられません。
自分が作ったプライマーが他にヒットする箇所があるかどうかを調べたい。
filePrimerSearch.pl

使い方↓

perl PrimerSearch.pl ~/Data/Ensembl_v88/cDNA/Homo_sapiens.GRCh38.cdna.all.fa_1line Xho1_SLC18A1_whole_Fw Not1_SLC18A1_whole_Rv

primer_Fとprimer_Rは配列を書いたファイルを用意。
以下のコマンドで調べたい種のfastaファイルをmakeblastdbでデータベース化してから使う。

makeblastdb -in ~/Data/Ensembl_v88/cDNA/Homo_sapiens.GRCh38.cdna.all.fa_1line -dbtype nucl -hash_index

アウトプットはこんな感じ↓

3'GC_content_PrimerF:   0.4
3'GC_content_PrimerR:   0.5
geneID	genename	 PrimerF_start	PrimerF_end	PrimerR_start	PrimerR_end	Length	GC_content_PrimerF	GC_content_PrimerR	Tm_PrimerF	Tm_PrimerR
ENST00000437980.3	SLC18A1	375	396	2033	2011	1659	0.545454545454545	0.521739130434783	59	59
ENST00000440926.3	SLC18A1	375	396	2167	2145	1793	0.545454545454545	0.521739130434783	59	59
ENST00000519026.5	SLC18A1	137	158	1833	1811	1697	0.545454545454545	0.521739130434783	59	59
ENST00000517776.5	SLC18A1	90	111	1984	1962	1895	0.545454545454545	0.521739130434783	59	59
ENST00000276373.9	SLC18A1	171	192	1963	1941	1793	0.545454545454545	0.521739130434783	59	59
ENST00000265808.11	SLC18A1	284	305	1980	1958	1697	0.545454545454545	0.521739130434783	59	59

Tm値とか3'の10塩基分のGC含量も計算しました。

#!/usr/bin/perl
#The way to calculate Tm value is obrained from http://catalog.takara-bio.co.jp/com/tech_info_detail.php?mode=1&masterid=M100005115 
use strict;
use warnings;


my $database=$ARGV[0];
my $primerFw=$ARGV[1];
my $primerRv=$ARGV[2];
my @primerF=();
my @primerR=();

open (IN4,$primerFw);
while(<IN4>){
chomp;
	@primerF=split(//,$_);
	my $GC3=0;
	for (my $a=-1;$a>=-10;$a--){
		$GC3++ if($primerF[$a] eq "G" or $primerF[$a] eq "C");
	}
	$GC3=$GC3/10;
	print "3\'GC_content_PrimerF:\t$GC3\n";
}

open (IN5,$primerRv);
while(<IN5>){
chomp;
	@primerR=split(//,$_);
	my $GC3=0;
	for (my $a=-1;$a>=-10;$a--){
		$GC3++ if($primerR[$a] eq "G" or $primerR[$a] eq "C");
	}
	$GC3=$GC3/10;
	print "3\'GC_content_PrimerR:\t$GC3\n";
}

my $blastF="blastn -outfmt 7 -evalue 1 -db $database -query $primerFw -out 
result_$primerFw -task blastn-short";
system($blastF);
my $blastR="blastn -outfmt 7 -evalue 1 -db $database -query $primerRv -out result_$primerRv -t ask blastn-short";
system($blastR);

open (IN,"./result_$primerFw");
open (IN2,"./result_$primerRv"); 
open (IN3,$database);

my %hash;
while(<IN3>){
chomp;
	if ($_=~/^\>(\S+)/){
		my $name=$1;
		if ($_=~/gene_symbol:(\S+)/){
			$hash{$name}=$1;
		}
	}
}

my %start;
my %end;
my %GC_F;
my %temp_F;
while(<IN>){
chomp; 
	if($_=~/^Query_1/){
		my @array=split(/\t/,$_);
		$start{$array[1]}=$array[8];
		$end{$array[1]}=$array[9];
		my ($count_AT,$count_GC)=(0,0);
		for (my $t=$array[6]-1;$t<$array[7];$t++){
			if ($primerF[$t] eq "A" or $primerF[$t] eq "T"){
				$count_AT++;
			}else{
				$count_GC++;
			}
		}
		$GC_F{$array[1]}=$count_GC/($count_AT+$count_GC);
		$temp_F{$array[1]}=4*$count_GC+2*$count_AT+35-2*($count_AT+$count_GC);
	}
}

print "geneID\tgenename\tPrimerF_start\tPrimerF_end\tPrimerR_start\tPrimerR_end\tLength\tGC_content_PrimerF\tGC_content_PrimerR\tTm_PrimerF\tTm_PrimerR\n";
while(<IN2>){
chomp;
	if($_=~/^Query_1/){
		my @array=split(/\t/,$_);
		if (defined $start{$array[1]} or defined $end{$array[1]}){
			my ($count_AT,$count_GC)=(0,0);
			for (my $t=$array[6]-1;$t<$array[7];$t++){
				if ($primerR[$t] eq "A" or $primerR[$t] eq "T"){
					$count_AT++;
				}else{
					$count_GC++;
				}
			}
			my $GC_R=$count_GC/($count_AT+$count_GC);
			my $temp_R=4*$count_GC+2*$count_AT+35-2*($count_AT+$count_GC);
			if($start{$array[1]} < $end{$array[1]} and $array[8] > $array[9]){
				my $length=$array[8]-$start{$array[1]}+1;
				print "$array[1]\t$hash{$array[1]}\t$start{$array[1]}\t$end{$array[1]}\t$array[8]\t$array[9]\t$length\t$GC_F{$array[1]}\t$GC_R\t$temp_F{$array[1]}\t$temp_R\n";
			}elsif($start{$array[1]} > $end{$array[1]} and $array[8] < $array[9]){
				my $length=$start{$array[1]}-$array[8]+1;
				print  "$array[1]\t$hash{$array[1]}\t$start{$array[1]}\t$end{$array[1]}\t$array[8]\t$array[9]\t$length\t$GC_F{$array[1]}\t$GC_R\t$temp_F{$array[1]}\t$temp_R\n";
			}
		}
	}
}

Tajima’s Dの計算(Window analysis)

#!/usr/bin/perl
use strict;
use warnings;

my $filename_in = $ARGV[0];	#input_vcf
my $chr_num = $ARGV[1];		#chr_num
my $filename_in2 = $ARGV[2];	#pop
my $filename_in3 = $ARGV[3];    #start_end_information_for_each_chr
my $window_first = $ARGV[4];	#window_first
my $window_last = $ARGV[5];	#window_last
my $window_size = $ARGV[6];	#window_size
my $window_step = $ARGV[7];	#window_step
my $filename_out = $ARGV[8];	#output

if($window_first!~/\-/ and $window_last!~/\-/ and $window_first>$window_last){
        print "window_firstはwindow_lastより小さい値を設定してください。\n" ;
        exit;
}

##各染色体のデータの始まりと終わりを取得
open(IN3,$filename_in3);
my %first_last;
while(<IN3>){
chomp;
next unless($_=~/(\S+)/);
	my @array=split(/\t/,$_);
	$first_last{$array[0]}{$array[1]}=$array[2];
}
close(IN3);
my $first=$first_last{$chr_num}{first};
my $last=$first_last{$chr_num}{last};
###############################

##解析に用いる集団のIDを取得
open(IN2,$filename_in2);
my %pop;
while(<IN2>){
chomp;
next unless($_=~/(\S+)/);
	$pop{$_}=1;
}
close(IN2);
my $n=(keys %pop)*2;
###############################

##Tajima's Dに必要な数字の計算
my $a1=0;
for (my $i=1;$i<$n;$i++){
	$a1=$a1+(1/$i);
}
my $a2=0;
for (my $i=1;$i<$n;$i++){
	$a2=$a2+1/$i/$i;
}
my $b1=($n+1)/(3*$n-3);
my $b2=2*($n*$n+$n+3)/(9*$n*$n-9*$n);
my $c1=$b1-1/$a1;
my $c2=$b2-($n+2)/($a1*$n)+$a2/$a1/$a1;
my $e1=$c1/$a1;
my $e2=$c2/($a1*$a1+$a2);
my $C=$n*($n-1)/2;
###############################

##windowの端の値を決める
$window_first=$first if($window_first=~/\-/);
$window_last=$last if($window_last=~/\-/);
print "SNPのある範囲にwindow_firstを設定してください。\n" if($first>$window_first or $last<$window_first);
print "SNPのある範囲にwindow_lastを設定してください。\n" if($last<$window_last or $first>$window_last);
print "window_sizeはwindow_stepの倍数に設定してください。\n" unless(($window_size % $window_step)==0);
my $temp3=(int($window_size/$window_step)-1)*$window_step;
my $first2=int($window_first/$window_step)*$window_step;
my $last2=int($window_last/$window_step)*$window_step;
my $n_temp=$n/2;
print "SNP_first=$first\twindow_first=$first2\nSNP_last=$last\twindow_last=$last2\nnumber_of_individual=$n_temp\n";
###############################

my $pos;
my %hash;
my %pop2;
my $present_window=$first2;
my ($pi,$D,$e,$e_t,$S,$num)=(0,0,0,0,0,0);
open(OUT,">$filename_out");
print "window_mid\twindow_first\twindow_last\tS\tPi\tTajimaD\n";
print OUT "window_mid\twindow_first\twindow_last\tS\tPi\tTajimaD\n";
open(IN,$filename_in);
while(<IN>){
chomp;
	if($_=~/^(\d+)\t/ or $_=~/^X\t/ and $_!~/^\#/){	#実際の多型データの行は染色体番号から始まる
		my @array=split(/\t/,$_);
		$pos=$array[1];
		next if($pos<$window_first);
		last if($pos>$window_last+$window_size);
		if($pos>$present_window+$window_size){
			my $window_end=$present_window+$window_size;
			my $window_mid=($present_window+$window_end)/2;
			($S,$e_t)=(0,0);
			foreach my $data (@{$hash{$present_window}}){
				$e=0;
				my @temp_array=split(/\t/,$data);
				$num=$n-1;
				for (my $c=0;$c<$num;$c++){
					for (my $d=$c+1;$d<=$num;$d++){
						$e++ if($temp_array[$c] ne $temp_array[$d]);
					}
				}
				$e_t=$e_t+$e;
				$S++ if($e>0);
			}
			$pi=$e_t/$C;
			if($S>0){
				$D=($pi-$S/$a1)/sqrt($e1*$S+$e2*$S*($S-1));
			}else{
				$D="nan";
			}
			print "$window_mid\t$present_window\t$window_end\t$S\t$pi\t$D\n";
			print OUT "$window_mid\t$present_window\t$window_end\t$S\t$pi\t$D\n";
			delete($hash{$present_window});
			$present_window=$present_window+$window_step;
		}
		my @array2=();
		foreach my $key (sort{$a<=>$b} keys %pop2){
			push @array2,$array[$key];
		}
		my $line=join("\t",@array2);
		$line=~tr/\|/\t/;
		my $temp=int($pos/$window_step)*$window_step;
		my $temp2=int($pos/$window_step+0.5)*$window_step;
		my $temp4=$temp-$temp3;
		for(my $a=$temp;$a>=$temp4;$a=$a-$window_step){
			if($pos>=$window_first and $pos<=$window_last){
				push @{$hash{$a}},$line;
			}
		}
	}elsif($_=~/^\#CHROM/){ #header行から見たい集団の個体の要素番号だけ得る
		my @array=split(/\t/,$_);
		for (my $a=0;$a<=$#array;$a++){
			my $name=$array[$a];
			$pop2{$a}=1 if(defined $pop{$name});
		}
	}
}
close(IN);
if($pos>$present_window+$window_size){
	my $window_end=$present_window+$window_size;
	my $window_mid=($present_window+$window_end)/2;
	($S,$e_t)=(0,0);
	foreach my $data (@{$hash{$present_window}}){
		$e=0;
		my @temp_array=split(/\t/,$data);
		$num=$n-1;
		for (my $c=0;$c<$num;$c++){
			for (my $d=$c+1;$d<=$num;$d++){
				$e++ if($temp_array[$c] ne $temp_array[$d]);
			}
		}
		$e_t=$e_t+$e;
		$S++ if($e>0);
	}
	$pi=$e_t/$C;
	if($S>0){
		$D=($pi-$S/$a1)/sqrt($e1*$S+$e2*$S*($S-1));
	}else{
		$D="nan";
	}
	print "$window_mid\t$present_window\t$window_end\t$S\t$pi\t$D\n";
	print OUT "$window_mid\t$present_window\t$window_end\t$S\t$pi\t$D\n";
	delete($hash{$present_window});
}
close(OUT);

inputの形式はVCFtoolsに倣った。 2016/7/24 更新。インプットが多くなったが、動作は少し速くなったはず。
今のところ、windowサイズをstepサイズの倍数にしないといけないのがネック。
染色体番号はそのまま"8"のように入力。集団の情報は、下みたいなファイルを作っておく。

HG02051
HG02325
HG02485
HG02308

また、以下のような各染色体のSNPデータの座位(最初と最後)情報をあらかじめ用意する。

1       first   10177
1       last    249240543
2       first   10179
2       last    243188367
3       first   60069

windowの開始位置と終了位置(というか含めたいSNPの座位)は数値で入力しても良いし、"-"と入力すれば最初(もしくは最後)の位置になる。
vcftoolsには勝てないが、2000万SNPsくらいあっても2日くらいでは終わるはず。



添付ファイル: filePrimerSearch.pl 44件 [詳細] filehaplotypes.png 200件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2017-11-09 (木) 10:38:08 (13d)