STARのSAMファイルパース時エラーの修正 (修正版)

[athos]

#16 の修正版のプルリクエストになります。

基本的には #16 と同様に、insertionを持つ融合遺伝子を検出した場合に発生するエラーを修正する変更です。
前回の修正では、read.seq という未定義変数の参照を避け、 F[9] (SAMファイルのSEQカラムに当たる)からinsertionの配列を抽出するようにしてエラーを修正しようしていました。

しかし、前回の修正ではアライメントの並びによっては正しいアライメントからinsertion配列を抽出できない場合がありました。具体的には、同一QNAMEを持つプライマリ、サプリメンタリ、およびそれらのペアからなるアライメントの3つ組が、入力SAMファイル中で、

プライマリ
サプリメンタリ
ペア

もしくは、

サプリメンタリ
プライマリ
ペア

の順に並んでいる場合、ペアのアライメントから誤ってinsertion配列を抽出してしまう問題がありました。

このプルリクエストでは、insertion配列を必ずプライマリのアライメントから抽出するように変更することで、このアライメントの順序依存性による問題を解決します。

[friend1ws]

ありがとうございます。ただ、このエラーはこれまで生じなかったので、、例えばこのエラーを起こすような共有できるテストデータなどあれば嬉しいですがいかがでしょうか？

[athos]

エラーを再現するためのサンプルデータを用意しました。

• pr17_repro.zip

こちらのzipファイルに以下のファイルを同梱しています：

• seq1.fq.gz, seq2.fq.gz
• realigned.fixed.sam

seq1.fq.gzおよびseq2.fq.gzは以下の融合遺伝子を検出するために人工的に作ったhg38用の配列です：

• CIC-DUX4
• ブレークポイント: 5'側 chr19:42295064 (strand: +), 3'側 chr4:190174745 (strand: +)
• ブレークポイント間のinsertion配列: GGGTGGAGGCAGTG

realigned.fixed.samは上記seq1.fq.gzおよびseq2.fq.gzをbwaでアライメントしたうえで、出力されたSAMをfusionfusionが期待するSTARのSAM形式に調整した結果のファイルです。手許の環境では、このSAMファイルを使って以下のようにエラーが再現することを確認しました：

$ fusionfusion --reference_genome <reference_file> --star realigned.fixed.sam --out <out_dir> --genome_id hg38 --no_blat --grc
Traceback (most recent call last):
  File "/usr/local/bin/fusionfusion", line 11, in <module>
    load_entry_point('fusionfusion==0.5.0', 'console_scripts', 'fusionfusion')()
  File "build/bdist.linux-x86_64/egg/fusionfusion/__init__.py", line 10, in main
  File "build/bdist.linux-x86_64/egg/fusionfusion/run.py", line 149, in fusionfusion_main
  File "build/bdist.linux-x86_64/egg/fusionfusion/parseJunctionInfo.py", line 596, in parseJuncInfo_STAR
  File "build/bdist.linux-x86_64/egg/fusionfusion/parseJunctionInfo.py", line 494, in getFusInfo_STAR
NameError: global name 'read' is not defined
$

[athos]

修正を最新のdevelブランチにリベースしました。

[aokad]

エラーが再現したのはbwaでアライメントしたうえで、出力されたSAMをfusionfusionが期待するSTARのSAM形式に調整した結果のファイル、とのことですが、fastqを直接STARでアライメントした結果のファイルをfusionfusionに入力していただけないでしょうか。
よろしくお願いいたします。

[athos]

コメントありがとうございます。

こちらの件について、#17 (comment) で共有致しましたサンプルデータをこちらでSTARでマッピングしてみましたが、Chimeric.out.samに期待するようなCICとDUX4に跨がるアライメントの組が出力されない(そのため、insertionも期待したように検出されない)ことが分かりました。

STARのオプションの調整で意図したアライメントの組が得られるようになるか、また別のfusionでよりエラーを再現しやすいものがあるか等、もう少しこちらで調べてみたいと思います。

[athos]

こちらの件について、改めて60件ほどinsertionを含むfusionの配列を準備して検証してみましたが、STARによって直接アライメントすることでエラーを再現できるサンプルは見つけられませんでした。

しかし、エラーが発生する箇所はfusionfusionがinsertionを検出した場合には必ず実行される処理で、サンプルが異なることで、あるいは将来的なSTARやfusionfusionの更新によってinsertionが検出可能になった場合、必ず未定義変数エラーを引き起こすことになるため、修正しておくのが望ましいではないかと考えています。

なお、このエラーは以下のようにPythonの静的解析ツール(pylint等)でも検出することができます:

$ pip install pylint
$ pylint fusionfusion | grep -E 'E[[:digit:]]+'
...
fusionfusion/parseJunctionInfo.py:494:30: E0602: Undefined variable 'read' (undefined-variable)
fusionfusion/parseJunctionInfo.py:566:30: E0602: Undefined variable 'read' (undefined-variable)
...
$