3-1. MiSeq i100からのデータが
うまく処理できない？
内容の要約
直面した問題：
Prokaryotic 16SやITSの次世代シーケンス結果の解析において、RのDADA2でエラー学習ができない
(リバースリードのみ)
原因：
300bpに満たないショートリードが取り除けていない場合に生じる無色・高品質Gがエラー学習を阻害？
2色法に基づくシーケンスで生成したfastqあれば生じ得る？
解決：
各FASTQファイルのリードのcomplexityを確認(dada2::seqComplexity + dada2::plotComplexity)
⇒通常であれば12~16ぐらいに集中、1や4あたりにややヒストグラムの棒が見える
(polyGとかが入っている可能性のある)Low complexityがある可能性
dada2::filterAndTrimでrm.lowcomplex=10を配置
コメント：
完全には除けないので、dada2後に一部キメラっぽいものは残っている
ただし、微量であり、また簡単にフィルターできる
詳細
fastqファイルの解析に先立ってquality scoreを確認すると異常にQ値が高い…
まぁ最近MiSeqからMiSeq i100に機械変わったしなぁ…と思い、そんなこともあるのかなと解析を継続。
(左2枚が問題のQ値、右2枚が普段のQ値)
Primerをtrimしてもあまり変わらず、やや気持ち悪さを抱えたまま進行。
(なお、primer trimmingはcutadapt v4.8を用いておりますが、--discard-untrimmedフラッグを立てております。
これで行けると思ってたんだけどなぁ…)
その後はlength filtering等無事に終了、いざdada2のerror学習するぞ！という段で以下のエラーとなりました。
Error rates could not be estimated (this is usually because of very few reads).
Error in getErrors(err, enforce =TRUE) : Error matrix is NULL
ここで、MiSeq i100になってからIllumina Adapter配列らしきものがあって
それ以降Gとなるreadが存在していてQiime2で難儀したことを思い出す。
adapter配列を頼りに探してみるとQualityがGのpolyGを今回も発見。
QualityがGなので、Q値38ですね。 こりゃぁだめだ。
慣れ親しんだMiSeqは4色法SBSという４色の傾向でATGCの塩基を判断するシステムだった一方、
新しいIlluminaのシステムは2色法SBSという２色でATGCを判別するシステムとなったことが原因のようです。
2色法SBSだと、色の付くC、TおよびCTの組み合わせで表現されるAは良いのですが…
short readで塩基がない部分が色のないGと判断される≒Qualityの高いGが生じてしまうようです。
PCRがうまく行っていてPrimer dimer等の短いリードが適切に除けている場合は気にするものではないのですが、
低濃度サンプルで頑張ってPCRして頑張ってLibraryを作ったうえに、SizeSelectionでも若干のコンタミが生じたときに
生じ得る問題でありそう？
記憶が正しければ通常は短いリードから読まれていくので、
Sequencing depthが十二分に確保できるような、Mixed libraryへの混合サンプル数が少ない状況では生じづらいかもしれません。
そのため、この問題は今後貧乏研究室の指標となるかもしれません笑
僕はずっと向き合っていくんだろうな…。
なお、--discard-untrimmedを立てても改善しなかった理由は、恐らくshort readもprimer配列を持っているからです。
primer dimerあるいはprimer自身のコンタミが究極の原因っぽいです。
問題の探索と解決
PolyG問題、どうするかというところですが、公式のIssueがやはりこういう時は参考になります。
https://github.com/benjjneb/dada2/issues/2147
実際同じ問題を抱えている人がおり、このissueに従いsequenceのcomplexityを描写すると
(計算の詳しい中身
・Biostrings::oligonucleotideFrequencyを使ってk-mer (defaultでは2mer)のカウント
oligonucleotideFrequencyは各kmerを列、各readを行にして出力
・各サンプルに対して kmerの出現頻度の割合を計算 (x/sum(x))
・これに基づき、shannon多様性を計算 sum(-y*ln(y))
・自然対数eのshannon乗したものをEffective Oligonucleotide Numberとして集計
・定義からshannon多様性は一番低い場合は１、一番高い場合はln (4^k)となる
≒つまり、Effective Oligonucleotide Numberは2 merだと1~16, 3 merだと1~64までの範囲を取る)
これで1(≒すべて同じ塩基)や4(かなり偏りがある)のところにバーが見えると、塩基に偏りが生じている
≒通常ではありえない
としてfilterAndTrim(…, rm.lowcomplex=X)で削除
X=10ぐらいとしておくのがよさそうらしい。
Effective Oligonucleotide Numberが10ということは、GG以外の2merが均等に出現していればGGの割合が40%ぐらいになる計算
(残りの15個の2kmerの出現頻度が均一ということはない(AG/TG/CG/GA/GT/GCが高くなる？)が)
なお、完全にランダムであれば16個のkmerは6.25%ぐらいの頻度で出現する
ここまですれば無事動く
が、完全には抜けきれないようで、例えば16Sだとrepresentative sequenceに対するlength filtering等で取り除くとよさそう
(現在中山研で使用している515F/806Rでは253bpが通常)
詰まったときにこれで解決できるようになることを祈っております。
3-2. R作図関連関数メモ
残念ながら(?)最近いろんな解説がggplot2をベースになされているようです。
でもなんかあれ、慣れないよなぁ…というロートル向けの関数メモのページです。
ggplot2を使わなくても自由に作図できるぞ！となりたいですね。
・基本的な設定
par()
作図に関するパラメータはここで設定します。
par(mar=c(a,b,c,d))
作図領域内で、図の余白を指定します。デフォルトはc(4,4,4,4)です。
右上の余白はいらないけど、左と下は軸ラベルを書くからある程度余白が欲しいなどなどの場合にこれをいじります。
a, b, c, dはそれぞれ下、左、上、右に対応しており、
par(mar=c(6,6,1,1))とすると左と下は余裕がある一方、右上の余白は狭くなります。
(0にするのはあまりお勧めしません。)
par(oma=c(a,b,c,d))
作図領域の外の余白です。１枚だけ図を描く場合はほぼ気にしなくて良いかと思います。
a, b, c, dの対応は同じです。
どういうときに使うのか、というと例えば図を2枚とか4枚とか同じに載せて、fig. 1 aとかで語りたい場合です。
４枚図を載せる場合、mar=c(6,6,1,1)と指定すると右上の図の左と下に余分な余白ができて見栄えが良くないです。
そこで、
par(oma=c(4,4,0,0)); par(mar=c(1,1,1,1))
など指定します。
この場合、各小さい作図領域のplot外の余白は1で指定されていて小さいですが、
全体として見た場合に左と下にいい感じの余白ができます。
mar, omaに関しては生物統計学さんのHPがわかりやすいです。
par(mfcol=c(a,b))
par(mfrow=c(a,b))
図を複数書くときに指定します。
c(a,b)の中のaとbは行・列に対応しており、例えばc(3,2)だと縦3枚横2枚の計6枚分の描画スペースが作られます。
mfcolは列方向、mfrowは行方向に描き進めていくもので、
c(3,2)の場合、2つ目のプロットは左列の2行目(mfcol)または右列の1行目(mfrow)となります。
par(fig=c(a,b,c,d))
ええいもっと自由に図を描きたいんだ！というときの設定です。
a,b,c,dはそれぞれ図のx方向の始点、x方向の終点、y方向の始点、y方向の終点を示しています。
a~dは0~1の範囲で示し、0が左端または下端、1が右端または上端です。
例えばc(0, 0.5, 0.5, 1)だと左上に全体の1/4サイズの図を描写することになります。
par(new=TRUE)
一度作った図に上書きしていくためのコマンドです。
例えば
par(mfcol=c(1,2))として
plot(~) #A
plot(~) #Bとした場合、左と右に1枚ずつ、計2枚の図が描かれます。
ただし、
par(mfcol=c(1,2))として
plot(~) #A
par(new=TRUE)
plot(~) #Bとした場合、左1枚、AとBが重なった図が描かれます。
これを活用することで、y軸を2軸にしたり棒グラフを折れ線グラフに重ねたり…とすることができます。
なお、par(fig=c(a,b,c,d))で描写位置を指定する場合で
他の場所に図を追加する場合もpar(new=TRUE)を記述する必要があります。
par(new=TRUE)はplot(~, add=TRUE)でも可能です。
他にもたくさんparの設定がありますが、主に使うのはこれぐらいかなという気がします。
他のものが気になった場合はRドキュメントを見るのが良いでしょう。