2015年2月12日木曜日

Rで変数を使いたくてハマったのでメモ

そこそこ数のあるファイルを読み込んだりする際に、read.tableするだけでも冗長になる
ES1 <- read.table("./ES1/quant_bias_corrected.sf")
ES2 <- read.table("./ES2/quant_bias_corrected.sf")
ES3 <- read.table("./ES3/quant_bias_corrected.sf")
PE1 <- read.table("./PE1/quant_bias_corrected.sf")
PE2 <- read.table("./PE2/quant_bias_corrected.sf")
PE3 <- read.table("./PE3/quant_bias_corrected.sf")
$とかつけて変数として扱ってfor文でまわせたら楽なのにとかおもって
sample <- c("ES1", "ES2", "ES3", "PE1", "PE2", "PE3")
for (i in sample) {
        i <- read.table("./(i)/quant_bias_corrected.sf")
}
とかしても当然ダメですハイ。
それぞれのオブジェクトを別々に作ってもいいけど、どうせその後の処理が一緒ならapply関数使えるようにlistに格納する方が無難。一度pathをリストに格納して読み込んだ表を別のリストに順に格納する。こうすればevalを使わなくてもなんとかなる。
num = 3
filename <- as.list(rep(NA, num))
cell <- as.list(rep(NA, num))

for(i in seq(num)){
    filename[[i]] <- paste("./ES", i, "/quant_bias_corrected.sf", sep="")
    cell[[i]] <- read.table(filename[[i]], sep="\t")
}

for(i in seq(num)){
    filename[[i+num]] <- paste("./PE", i, "/quant_bias_corrected.sf", sep="")
    cell[[i+num]] <- read.table(filename[[i+num]], sep="\t")
}
names(cell) <- c("ES1", "ES2", "ES3", "PE1", "PE2", "PE3")

2014年11月3日月曜日

BEDファイルのソート

BEDファイルがでたらめにソートしていると、
chr1
chr10
chr11
chr12
…
chr19
chr2
chr20
chr21
chr22
chrX
chrY
という感じになっている。つまり、chr~ごと文字として認識しているから。
sort -k 1,1 -k2,2n hoge.bed 

#とすれば、以下のように理想的な形にソートされる。
chr1
chr2
chr3
…
chr21
chr22
chrX
chrY

2014年9月2日火曜日

MacでホストOSとゲストOSでファイルを共有する方法 on VirtualBox

だいたいググって出てくる結果と変わらない。

  1. 正常にMacにVirtualBoxがインストールされていれば、 ~/VirtualBox VMsができてる。
  2. ここに適当にVM-Shareという適当なディレクトリを作成する。
  3. VirtualBoxを起動して、設定の共有フォルダで右側の+をクリック
  4. 先ほど作ったディレクトリのパスを入力して、適当にshareと名前を付ける。
  5. 一度VirtualBoxを落として再起動、適宜イメージを起動。
  6. ターミナルで、sudo mkdir /mac
  7. パス聞かれるから入力。入力しても文字でないけど。気にしない。
  8.  sudo mount -t vboxsf share macとすれば、ゲスト側のmacというディレクトリとホスト側のVM-Shareが共有される。

2014年8月11日月曜日

liftover

既に解析済みのbedやgif, positionデータはあるものの、現在自分が解析しているゲノムのバージョンとは異なる場合があるが、UCSCにあるliftoverでそろえる事が可能。

http://hgdownload.cse.ucsc.edu/admin/exe/で環境に合わせて
私の場合はlinux.x86_64/以下にあるliftoverをダウンロードchmod +x liftoverでlocal/binなどパスが通っている場所に置いておく。さらに、例えば、mm8からmm10に変換したい場合は、http://hgdownload.soe.ucsc.edu/goldenPath/mm8/liftOver/mm8ToMm10.over.chain.gzをダウンロードして解凍。適当に共有ディレクトリ等に置いておく。

#変換したいファイルがbedなら(chr1 1000099 1000200)
liftOver Oct4_mm8.bed mm8ToMm10.over.chain Oct4_mm10.bed unmapped.txt
Oct4_mm10.bed>>変換されたbed 
unmapped.txt>>変換できなかったもの
#変換したいファイルがpositionなら(chr1:1000100-1000200)
liftOver -positions Oct4_mm8.txt mm8ToMm10.over.chain Oct4_mm10.txt unmapped.txt
Oct4_mm10.txt>>変換されたposition
unmapped.txt>>変換できなかったのもの
liftOver_XXXXXX.bed>>変換前のbedファイル
liftOver_XXXXXX.bedmapped>>変換後のbedファイル
positionでもliftover出来て、そのbedを得られるのは便利

GEOのRAWデータをターミナルからダウンロード

あるアクセッション番号のRAWデータをダウンロードしたい時に、webのリンクからだと直接ダウンロードできないので、以下の様に自分のダウンロードしたい番号を入れる。親ディレクトリは下三桁をnnnとして、それより上の桁をダウンロードしたい番号と同じにする。
以下は例
wget ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE11nnn/GSE11431/suppl/GSE11431_RAW.tar

2014年7月25日金曜日

サーバーの特定のディレクトリをマウントする

サーバー側で作成したpdfやpngなどのファイルをいちいちscpするのはめんどくさいのでマウントして手元のPCから直接覗けるようにする。 基本的にhttp://blog.offline-net.com/2014/03/02/mavericks-osxfuse-sshfs-local2remoteをみればわかる 私も上記解説と同じ環境である。 デスクトップにhogeという名前のフォルダを先に作成する。 以下の様にターミナルで叩けばマウントできる。
sshfs ユーザー名@サーバー名(hoge.co.jp):/home/ユーザー名/PATH/TO/指定ディレクトリ /Users/ユーザー名/Desktop/hoge
指定ディレクトリ以下にディレクトリがあってもアクセスできる。

fastaファイルの特定の位置をマスクする方法

bedtoolsの中のmaskFastaFromBedというツールを使う。 使い方は
maskFastaFromBed -fi [マスクかけたいfasta ファイル] -bed  [位置を指定した bedファイル] -fo [ 出力fa staファイル名] -mc N(Nでマスクする)
マスクかけたいfasta ファイルはmulti fastaでOK bedも複数行でOK 例として
test.fa
>chr1
TCTCATGCTAATCAGGATCTGGGGCACAGTACAAACCTCAAATGCTGTAGTTGAG
>chr2
GCTATTTCAGCCCTAGTGGAAGGGGGTGGGGCGGCACTAGGGCAGCAGCCCGTGCCCGCTGCACTTCCCTTATCGGATTACTGCTTAACTCCCACCCCCTCCTAGGTCCAAAAAGTATACTAATTACATGGA
>chr3
GCAGTACATTCCACAGCATCCTTAGGTGTTGCTCAGAGTCTACCTAACACCTAAGTACAAAAATGAAGCCCTAGGGCAGCCCTAAATGGTTTGTTTCCGCAAAATAGTGATCAACTCAAGTAAACAAGGGTT
>chr4
CATAAGACACATTGAAGACACCGCCACAGAACTTTACTGCCCTCAAGAGTCTTCTCTCATCTCTGGGATGCCTGTGAGGC

test.bed
chr1 5 10
chr2 10 20
chr2 21 23
chr5 10 20

maskFastaFromBed -fi test.fa -bed  test.bed -fo out.fa -mc N

out.fa
>chr1
TCTCANNNNNATCAGGATCTGGGGCACAGTACAAACCTCAAATGCTGTAGTTGAG
>chr2
GCTATTTCAGNNNNNNNNNNANNGGGTGGGGCGGCACTAGGGCAGCAGCCCGTGCCCGCTGCACTTCCCTTATCGGATTACTGCTTAACTCCCACCCCCTCCTAGGTCCAAAAAGTATACTAATTACATGGA
>chr3
GCAGTACATTCCACAGCATCCTTAGGTGTTGCTCAGAGTCTACCTAACACCTAAGTACAAAAATGAAGCCCTAGGGCAGCCCTAAATGGTTTGTTTCCGCAAAATAGTGATCAACTCAAGTAAACAAGGGTT
>chr4
CATAAGACACATTGAAGACACCGCCACAGAACTTTACTGCCCTCAAGAGTCTTCTCTCATCTCTGGGATGCCTGTGAGGC
となるはず。 bedは定義がstartの数字を含まないので5と書いてあれば6文字目を示しendは10と書いてあれば10文字目を示す事に注意する。