ハーバード留学研究2年目編010：ちょいしょぼNGS解析覚書2

そんなに重要な解析ではないが、QCの一環として、それぞれ別々の蛍光タンパク（Tdimer2とCerulean）を発現する細胞から、とってきたRNA-seqデーター(74.bam, 75.bam)から本当に前情報と同じ蛍光タンパクを発現しているかどうかチェックしてみる。

前回行った解析手法を利用してみる。Tdimer2とCeruleanのfastaqファイル及びindexファイル作成は前回に準じる。

１）BAM fileからunmapped sq抽出
$ samtools viwe -f4 -b 74.bam | samtools sort -74-unmp
$ samtools viwe -f4 -b 74.bam | samtools sort -74-unmp

$ samtools index 74-unmp.bam
$ samtools index 75-unmp.bam

2)Tdimer アラインメント

$ bwa aln -b Td.fa 74-unmp.bam -> 74_T.sai
$ bwa aln -b Td.fa 75-unmp.bam -> 75_T.sai

3)SAI -> SAM変換

$ bwa samse Td.fa 74_T.sai 74-unmp.bam > 74_T.sam
$ bwa samse Td.fa 75_T.sai 74-unmp.bam > 75_T.sam

4)SAM->BAM変換

$ samtools view -bS 74_T.sam > 74_T.bam
$ samtools view -bS 75_T.sam > 75_T.bam

5)ソート
$ samtools sort 74_T.bam 74_T_sort
$ samtools sort 75_T.bam 75_T_sort

6)インデックス作成

$ samtools index 74_T_sort.bam
$ samtools index 75_T_sort.bam

7)結果表示

$ samtools idxstats 74_T_sort.bam
$ samtools idxstats 75_T_sort.bam

てな具合である。
主力結果は

74が
Tdimer 1395 1117 17
* 0 0 15976919

75が
Tdimer 1395 0 0
* 0 0 17205327

のようにでてくる。
このうち2つ目のカラムの上が、マップされたリード数だから（詳しくはsamtools使い方参照）、74のunmapped sequence のうち、1117リードがTdimerの配列にマップされた、一方で75のunmapped sequence のうちTdimerの配列にマップされたものはなかったということになる。

Ceruleanへのアラインメント及び結果表示も同様。

74が
Ceru 825 7 0
* 0 0 17205454

75が
Ceru 825 4451 114
* 0 0 15972354
となっている。74のほうは4リードが、75のほうは4451リードがunmapped sequenceからCeruleanの配列にマップされたことになる。

定性的ではあるが、これで74がTdimer, 75がCeruleanを発現していることがわかり、前情報とサンプルの中に発現している蛍光タンパクは同じであると考えられる。ふー。。