sshfs ユーザー名@サーバー名(hoge.co.jp):/home/ユーザー名/PATH/TO/指定ディレクトリ /Users/ユーザー名/Desktop/hoge指定ディレクトリ以下にディレクトリがあってもアクセスできる。
2014年7月25日金曜日
サーバーの特定のディレクトリをマウントする
サーバー側で作成したpdfやpngなどのファイルをいちいちscpするのはめんどくさいのでマウントして手元のPCから直接覗けるようにする。
基本的にhttp://blog.offline-net.com/2014/03/02/mavericks-osxfuse-sshfs-local2remoteをみればわかる
私も上記解説と同じ環境である。
デスクトップに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文字目を示す事に注意する。
2014年7月24日木曜日
マップしなかったリード・マップしたリード・ユニークなリードの数え方
bowtieが終わったあとに表示される以下の様なリード数が表示されるがこれをあとでコマンドから確認する方法
<<<<<追記>>>>>>
XS:による方法は非常に限られた状況でのみ成立するっぽいです.
bowtie2でSEでやった時は少なくともbowtieを走らせたあとの結果と一致しました。
ちなみにPEでやったものを同様に行うと、想定よりも大きい値が出たため、PEにはこの方法は適応できません。
XS:i:<N>
Alignment score for the best-scoring alignment found other than the alignment reported. Can be negative. Can be greater than 0 in --local mode (but not in --end-to-end mode). Only present if the SAM record is for an aligned read and more than one alignment was found for the read. Note that, when the read is part of a concordantly-aligned pair, this score could be greater than AS:i. ってbowtie2のmanualに記載があるけれどよくわかりません( TДT)
SE/PEどちらの場合も、よくいわれるNH:i:のフラグはbowtie2では立ちませんでした。MAPQの値で行う事も出来ませんでした。
確実にmulti alignを取り除く方法知っている方いたら教えてくださいm(_ _ )m
135582741 reads; of these: 135582741 (100.00%) were unpaired; of these: 22383842 (16.51%) aligned 0 times 61518065 (45.37%) aligned exactly 1 time 51680834 (38.12%) aligned >1 times 83.49% overall alignment rate
基本的にフラグが立っているのでそれを利用する。 #reads; of these:これはfastQCした結果のhtmlを見るかbowtieにかける前のfastqをwc -lして4で割れば良い。 #aligned 0 times:マップされなかったリード数の表示は samtools view -cf 4 filename.bam #マップされたリード数の表示は samtools view -cF 4 filename.bam #aligned >1 times:マルチマップされたリード数の表示は grep 'XS:' ./filename.sam | wc -l #aligned exactly 1 time:ユニークにマップされたリード数の表示は マップされたリード数からマルチマップされたリード数を引くか sed '/XS:/d' Enh1.sam | wc -lとして、全リード数からマルチマップのリードをsamから取り除いた行カウントを行う。その値からsamのheaderのライン数とマップされなかったリード数を引く。ちなみに、pairをきちんと組んでいるものを抽出するには この組み合わせ
samtools view -S -f 99 -S filename.sam > ./filename.99.sam samtools view -S -f 147 -S filename.sam > ./filename.147.samこの組み合わせ
samtools view -S -f 83 -S filename.sam > ./filename.83.sam samtools view -S -f 163 -S filename.sam > ./filename.163.samhttp://ppotato.wordpress.com/2010/08/25/samtool-bitwise-flag-paired-reads/がわかりやすい。
<<<<<追記>>>>>>
XS:による方法は非常に限られた状況でのみ成立するっぽいです.
bowtie2でSEでやった時は少なくともbowtieを走らせたあとの結果と一致しました。
ちなみにPEでやったものを同様に行うと、想定よりも大きい値が出たため、PEにはこの方法は適応できません。
XS:i:<N>
Alignment score for the best-scoring alignment found other than the alignment reported. Can be negative. Can be greater than 0 in --local mode (but not in --end-to-end mode). Only present if the SAM record is for an aligned read and more than one alignment was found for the read. Note that, when the read is part of a concordantly-aligned pair, this score could be greater than AS:i. ってbowtie2のmanualに記載があるけれどよくわかりません( TДT)
SE/PEどちらの場合も、よくいわれるNH:i:のフラグはbowtie2では立ちませんでした。MAPQの値で行う事も出来ませんでした。
確実にmulti alignを取り除く方法知っている方いたら教えてくださいm(_ _ )m
Trimmomatic
だいぶ慣れてきたので、インストールの仕方よりもツールの使い方をメインに記載するかな。
アダプタートリマーツールであり、クリーニングも出来るTrimmomatic (http://www.usadellab.org/cms/?page=trimmomatic)
使い方としては、デフォルトにadoptersにいくつか入っているが、適当に自作してよい。例えば、Truseq4-SE.faとして
このツールはシーケンスするフラグメントが非常に短い時にインサートを突き抜けてアダプターもシーケンスしてしまった時の除去に非常に有効なツール
>TruSeq_v1v2_Index12Adapter GATCGGAAGAGCACACGTCTGAACTCCAGTCACCTTGTAATCTCGTATGCCGTCTTCTGCTTG TruSeq_v1v2_UniversalAdapter AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCTというファイルを作って、以下のコマンドを叩く。
java -jar Trimmomatic-0.32/trimmomatic-0.32.jar SE \ -threads 8 -phred33 \ -trimlog ./DNseqE14.log.txt \ ./DNseqE14.fastq \ ./trimmomatic/DNseqE14.trim.fastq \ ILLUMINACLIP:../../../../../local/Trimmomatic-0.32/adapters/TruSeq4-SE.fa:2:40:15 \ LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:20ILLUMINACLIP自体がアダプタートリム。その他はクオリティ値でトリム。 ペアエンドの時は少し書き方が異なるがhttp://www.usadellab.org/cms/?page=trimmomaticに参考例の記載がある。
このツールはシーケンスするフラグメントが非常に短い時にインサートを突き抜けてアダプターもシーケンスしてしまった時の除去に非常に有効なツール
登録:
投稿 (Atom)