Genomon-Project / fusionfusion

script for detecting fusion genes from several transcript alignment tools
GNU General Public License v3.0
9 stars 8 forks source link

中間fastaファイルに混入されるリージョン情報の除去 #10

Closed athos closed 3 years ago

athos commented 4 years ago

現状のfusionfusionでは、中間ファイルとして生成されるfastaファイル(*.chimeric.clustered.splicing.contig.fa)にリージョンの情報のようなものが意図せず含まれてしまっているようです。たとえば、以下のようなfastaが生成されます:

>1:-47767000-1:+47770700_contig1
1:47767000-47767050AAGTACATCTTGAACTTTCTTATCTCGGTTACACTGCTTCTTCTGAAAGGC
>1:-47767000-1:+47770700_contig2
TGTGCTTTAATTTTACATTGTTTTTTCAAAAGGTTTCCTTCAAGCAGGATG00707774-05607774:1

2行目の先頭部および4行目の末部が該当箇所です。

これは、fastaのシーケンス文字列を構築する seq_utils.getSeq() とその内部で使用される pysam.faidx() で想定するデータ構造が食い違っていることに起因しているようです。seq_utils.getSeq()pysam.faidx() が結果を改行区切りのリストとして返すことを想定していますが、実際には pysam.faidx() は結果を改行では区切らずひとかたまりの文字列として返します。 その結果、seq_utils.getSeq() 内にあるリージョンの行を除く処理がうまく機能しなくなっているようです:

>>> import pysam
>>> import fusionfusion.seq_utils as seq_utils
>>> pysam.faidx('reference.fa', "1:47767000-47767050")
'>1:47767000-47767050\nAAGTACATCTTGAACTTTCTTATCTCGGTTACACTGCTTCTTCTGAAAGGC\n'
>>> seq_utils.getSeq('reference.fa', [("1", 47767000, 47767050)])
'1:47767000-47767050AAGTACATCTTGAACTTTCTTATCTCGGTTACACTGCTTCTTCTGAAAGGC'
>>> 

このプルリクエストは、pysamの split_lines を使って、結果を改行区切りにすることでシーケンス文字列部分のみを抽出するように修正します:

>>> pysam.faidx('reference.fa', "1:47767000-47767050", split_lines = True)
['>1:47767000-47767050', 'AAGTACATCTTGAACTTTCTTATCTCGGTTACACTGCTTCTTCTGAAAGGC']
>>> seq_utils.getSeq('reference.fa', [("1", 47767000, 47767050)])
'AAGTACATCTTGAACTTTCTTATCTCGGTTACACTGCTTCTTCTGAAAGGC'
>>> 

従前の方法により生成されるfastaファイルでも、それを受け取る後段のBLATではシーケンス文字列に現れるアルファベットとハイフン以外の文字を読み捨てるようになっているため、fusionfusion全体としての解析結果に与える影響は限定的なものと思われますが、シーケンス文字列中のハイフンはギャップと解釈されるため影響がまったくないとも言い切れず、修正をした方が懸念が少ないかと思います。