追記
1、クオリティチェック
変なリードが混じっていると、本来シーケンスされた領域以外の遺伝子にマッピングされたりして解析結果を歪める恐れがある。そのため、シーケンスが終わりBCLからfastqデータが出力されたら、一度クオリティをチェックする。
クオリティがイマイチな場合、リード低クオリティ領域のトリミングや、リードの除去を行う必要がある。
トリミングはphread quality scoreに基づいて行う。クオリティ分布をグラフ化したり、指定した閾値にしたがって5'側や3'側をトリムするソフトはたくさん存在するし、有償の次世代データ解析ソフトにもほぼ間違いなく実装されている。
オススメはfastx-toolkit
使い方だが、
https://bi.biopapyrus.net/rnaseq/qc/fastx-toolkit.html
の東大の方？の運営されているHPにわかりやすくまとめられている。
（いつも勉強させてもらっています）
インストールはhomebrewを使って
brew install homebrew/science/fastx_toolkit
使い方は先ほどのページ通りです。
Rならqrqcが色々見れて便利。
ターミナルを起動し、シーケンスfastqがあるディレクトリに移動してRを起動(osx 10.11以上でRを標準搭載）
$R #R起動のコマンド
>install.packages("qrqc", dependencies = TRUE)
library(qrqc) #ライブラリの読み込み
fq <- readSeqFile("DRR016695_1_fill.fastq") #データの読み込み
makeReport(fq)
これで実行ディレクトリにフォルダができているはず。その中に下のデータが全て保存されている。
勉強を兼ねて個別に実行してみる。
クオリティ分布を見るには
qualPlot(fq)
塩基の出現頻度は
basePlot(fq)
リード長の分布
seqlenPlot(fq)
K-mer 頻度を見るには（アセンブリの参考になる）
kmerKLPlot(fq)
複数ウィンドウを残しながら進めるなら
dev.new() #これで新しい白紙ウィンドウができる。
qualPlot(fq)
dev.new()
basePlot(fq)
dev.new()
seqlenPlot(fq)
dev.new()
kmerKLPlot(fq)
結果は以下のようになる。
あとから追加編集したくなった場合、例えば2番ウィンドウをアクティブにするには dev.set(2)
アクティブなデバイスを閉じる
dev.off()
すべてのデバイスを閉じる
graphics.off()
実際のトリミングはastx-toolkitで行えばよい。処理スピードも早いし、柔軟なトリミングが可能。問題はどこまで貴重なデータを捨てるか？ということになる。真核生物のRNA seqでは最低2000万リードは必要と言われたりする。これは細胞内で100万分の１分子数だけ発現している（=低発現のmRNA）を２０リードシーケンスでき、ある程度の統計的処理に対応できる数ということをどこかで聞いた。
時々低発現のmRNAはstochasticな発現なのかよく分からんし高発現のものだけ捉えられたらいいというような発現を聞くことがあるが、ome解析なのだから低発現のRNAも含めてシーケンスされているのが理想である。場合によっては隠れた登場人物を見過ごして結論を出してしまう恐れもあるのだから。
ということで、トリミングを欲張ればデータがなくなるし、トリムを全くしないと、エラーリッチなリード異なるmRNAにマッピングされる恐れがある。phread quality scoreが20の場合、経験則で1%の確率でエラーであることがわかっている(画像データから塩基を決定する時のアナログ->デジタル変換時の間違いなど？）。
1%というと低いようだが、100bpシーケンスすれば、１リードにつき１塩基は間違っている可能性が高いことになる。Q=10なら エラー率10%となる。100bpシーケンスすれば、10塩基くらい間違っている。リピートに対して十分な精度があると言えるだろうか？
ということで、Q=20くらいでトリムする人が多いのではないだろうか？
ただし個人が開発したプログラムによってはペアが同一数あるのを前提に作られているソフトもある。トリムすると、ペアの個数が異なりその後の処理が行えない場合もある。その時はhomemadeなスクリプトを書いてペアードエンドのリード数を少ない方に合わせる必要がある。私は海外の次世代フォーラムみつけたperlのスクリプトを使ってます。実行するには
perl To_Remove_Unpaired_Reads.pl trimmed_R1.fastq trimmed_R2.fastq
To_Remove_Unpaired_Reads.plは自分がコピペして保存した適当な名前。
引数としてトリムして数が異なってしまったR1とR2のfastqを指定する。
ただし、ペアの個数がより現実的な問題になるには、何らかの理由によってR2だけクオリティが悪くてR2のリードがR1より随分減ったときなどである。ペアリードとしてマッピングするなら、数が少ない方にリード数を揃える必要があり、ますますリードが減ってしまう。ぶっちゃけそんなデータを解析に使うな、となるが、高額の研究費を使って初めてRNA seqやる人ならそこで諦める人はいないだろう。
ということでトリミングなしでマッピングしてまた問題になったりする。書いているときりがないが、クオリティ以前の問題で、非モデル生物からRNAを無理やりとってRIN値の低いサンプルを使ってRNA seqすることもあるかもしれない。rRNA/mRNA値が変わると特に低発現のmRNAの分布のばらつきが大きくなることがわかっている。
RNA seqをしっかりやるには簡単ではない。
脱線してしまったが、クオリティ解析が終わればマッピング解析を行う。それは次回に
