共同研究でRNA-seqデーター解析をしており、性懲りもなくデーターをもらってきて勉強がてら解析してみることにした(*)。
Pair-endのデーターなので
> bowtie -S -p 2 mm9 -1 X_1.fastq -2 X_2.fastq X.sam (**)
でいいはず。
しかし、ほとんどがunmappedになってしまう。
ん？とりあえず、トラブルシュートを検討します。
とりあえず片側だけでやった方が、unmappedがすくなそうなので(***)、
ためしてみる。
> bowtie -S -p2 mm9 X_1.fastq X_1.sam
としunmapped sequenceを抽出してみることとする。
★Sam->Bam変換
> samtools view -bS X_1.sam > X_1.bam
★unmappedシーケンスの抽出
> samtools viwe -f4 -b X_1.bam | samtools sort -X_1-unmp (****)
> samtools index X_1-unmp.bam
と遅くなってきたので、今日はここまで。。
PCが遅い。。(涙）
(2014年9月3日追記）
★多色蛍光マウスなので、CAG－GFPの配列にあててみることとする。
詳しくは、STAP細胞関係のゲノムデータを解析してみた7619にある通りなのだが、
１）CAG-GFPの配列をとってき、fast形式のの配列をつくる（CAG.fa）（ここはSTAP細胞関係のゲノムデータを解析してみた7619を参照のこと）
２）インデックスをつくる。
>bwa index -a bwtsw cag.fa
３）次にアラインメントするが、ここはSTAP細胞関係のゲノムデータを解析してみた7619の通りだとなぜかうまくいかなかった。
そこでお手本のように引数を渡す手法でなく、分割しておこなった。
本当はbowtieが使えるとよいのだけれど、bamファイルでの入力ができない。
まずアラインメント
>bwa aln -b CAG.fa X_1_unmp.bam -> X_1-cag.sai
次にsai -> sam変換
>bwa samse CAG.fa X_1-cag.sai X_1_unmp.bam > X_1-cag.sam
sam -> bam変換
>samtools view -bS X_1-cag.sam > X_1-cag.bam
ソート
>samtools sort X_1-cag.bam X_1-cag_sort
インデックス作成
>samtools index X_1-cag_sort.bam
で無事作業は終わっていると思う。
最後にアラインメントされたところのみをみると、
>samtools view -F4 X_1-cag_sort.bam
のように出力できる。
ということで今日はここまで。
（2014年9月18日＆10月18日）
出力は見にくいので、
>samtools view -F4 X_1-cag_sort.bam > X_1-cag.txt
みたいにファイル出力して、あとでエクセルで見た方がよいかも。
また単にあたっているかどうか、見るだけなら、
>samtools idxstats X_1-cag_sort.bam > X_1-cag.txt
でもよいかもしれない。
列番号 説明 例
1 リファレンスの配列名 CAG
2 リファレンスの配列長 6118
3 マップリード数 311
4 非マップリード数 0
のように何本あたっているかが表示される。詳細はsamtoolsの使い方参照。
(*)いつぞやのSTAPデータでは、メモリ不足でこりました。今回は多少Ubuntuを高速化したりしてリトライです。何とかネットブックでも遊べそうな勢いです。何でSTAPデータはあんなに巨大だったのか。。
(**)Xはサンプル名です。
(***)何かやり方を間違っているのだろう。それでも30％くらいある。
(****)詳しくは以下のサイトが参考になる。
１）STAP細胞関係のゲノムデータを解析してみた7619
２）sam形式ファイルsam format
ちなみにsamtoolsのオプションの説明はsam形式ファイルsam formatを参考に。
-S はinput fileがsam形式であることを指定している。
-f は許可する数字，この場合4なのでマップされていないリードを許可します。
-Fは許可しない数字，この場合194（128+64+2）なのでpaied endでマップされているリードは許可しない
その他のオプションの説明はsamtoolsの使い方
より
Options: -b output BAM
-h print header for the SAM output
-H print header only (no alignments)
-S input is SAM
-u uncompressed BAM output (force -b)
-x output FLAG in HEX (samtools-C specific)
-X output FLAG in string (samtools-C specific)
-c print only the count of matching records
-t FILE list of reference names and lengths (force -S) [null]
-T FILE reference sequence file (force -S) [null]
-o FILE output file name [stdout]
-R FILE list of read groups to be outputted [null]
-f INT required flag, 0 for unset [0]
-F INT filtering flag, 0 for unset [0]
-q INT minimum mapping quality [0]
-l STR only output reads in library STR [null]
-r STR only output reads in read group STR [null]
-? longer help
Pair-endのデーターなので
> bowtie -S -p 2 mm9 -1 X_1.fastq -2 X_2.fastq X.sam (**)
でいいはず。
しかし、ほとんどがunmappedになってしまう。
ん？とりあえず、トラブルシュートを検討します。
とりあえず片側だけでやった方が、unmappedがすくなそうなので(***)、
ためしてみる。
> bowtie -S -p2 mm9 X_1.fastq X_1.sam
としunmapped sequenceを抽出してみることとする。
★Sam->Bam変換
> samtools view -bS X_1.sam > X_1.bam
★unmappedシーケンスの抽出
> samtools viwe -f4 -b X_1.bam | samtools sort -X_1-unmp (****)
> samtools index X_1-unmp.bam
と遅くなってきたので、今日はここまで。。
PCが遅い。。(涙）
(2014年9月3日追記）
★多色蛍光マウスなので、CAG－GFPの配列にあててみることとする。
詳しくは、STAP細胞関係のゲノムデータを解析してみた7619にある通りなのだが、
１）CAG-GFPの配列をとってき、fast形式のの配列をつくる（CAG.fa）（ここはSTAP細胞関係のゲノムデータを解析してみた7619を参照のこと）
２）インデックスをつくる。
>bwa index -a bwtsw cag.fa
３）次にアラインメントするが、ここはSTAP細胞関係のゲノムデータを解析してみた7619の通りだとなぜかうまくいかなかった。
そこでお手本のように引数を渡す手法でなく、分割しておこなった。
本当はbowtieが使えるとよいのだけれど、bamファイルでの入力ができない。
まずアラインメント
>bwa aln -b CAG.fa X_1_unmp.bam -> X_1-cag.sai
次にsai -> sam変換
>bwa samse CAG.fa X_1-cag.sai X_1_unmp.bam > X_1-cag.sam
sam -> bam変換
>samtools view -bS X_1-cag.sam > X_1-cag.bam
ソート
>samtools sort X_1-cag.bam X_1-cag_sort
インデックス作成
>samtools index X_1-cag_sort.bam
で無事作業は終わっていると思う。
最後にアラインメントされたところのみをみると、
>samtools view -F4 X_1-cag_sort.bam
のように出力できる。
ということで今日はここまで。
（2014年9月18日＆10月18日）
出力は見にくいので、
>samtools view -F4 X_1-cag_sort.bam > X_1-cag.txt
みたいにファイル出力して、あとでエクセルで見た方がよいかも。
また単にあたっているかどうか、見るだけなら、
>samtools idxstats X_1-cag_sort.bam > X_1-cag.txt
でもよいかもしれない。
列番号 説明 例
1 リファレンスの配列名 CAG
2 リファレンスの配列長 6118
3 マップリード数 311
4 非マップリード数 0
のように何本あたっているかが表示される。詳細はsamtoolsの使い方参照。
(*)いつぞやのSTAPデータでは、メモリ不足でこりました。今回は多少Ubuntuを高速化したりしてリトライです。何とかネットブックでも遊べそうな勢いです。何でSTAPデータはあんなに巨大だったのか。。
(**)Xはサンプル名です。
(***)何かやり方を間違っているのだろう。それでも30％くらいある。
(****)詳しくは以下のサイトが参考になる。
１）STAP細胞関係のゲノムデータを解析してみた7619
２）sam形式ファイルsam format
ちなみにsamtoolsのオプションの説明はsam形式ファイルsam formatを参考に。
-S はinput fileがsam形式であることを指定している。
-f は許可する数字，この場合4なのでマップされていないリードを許可します。
-Fは許可しない数字，この場合194（128+64+2）なのでpaied endでマップされているリードは許可しない
その他のオプションの説明はsamtoolsの使い方
より
Options: -b output BAM
-h print header for the SAM output
-H print header only (no alignments)
-S input is SAM
-u uncompressed BAM output (force -b)
-x output FLAG in HEX (samtools-C specific)
-X output FLAG in string (samtools-C specific)
-c print only the count of matching records
-t FILE list of reference names and lengths (force -S) [null]
-T FILE reference sequence file (force -S) [null]
-o FILE output file name [stdout]
-R FILE list of read groups to be outputted [null]
-f INT required flag, 0 for unset [0]
-F INT filtering flag, 0 for unset [0]
-q INT minimum mapping quality [0]
-l STR only output reads in library STR [null]
-r STR only output reads in read group STR [null]
-? longer help