Paired-endのデータが不揃いの場合に起こるエラーについて

mailing listに投稿された以下のエラー

Message posted: The data seems to be corrupt or originates from different sources since the input files contain different numbers of reads!

これは、paired-endの2つのfastqファイルのread数が合わない場合にその後の解析でたまにみるエラー。fastqファイルを事前にQCフィルタリングするなどすると、２つのファイルが不揃いになる。では、fastqファイルをQCせずにmappingするのか？というと、少し気持ち悪く思う人もいるし、タグカウントするのみであればmapping qualityでいいんじゃないか？と思う人もいると思う。また、件の場合には、シークエンシングサービスがQC後のデータを提供してくれているなどで、元データの提供がない場合もあるだろう。やはり、「一から全部自分でやるのが一番！」と思う瞬間でもある。

さてせっかくなので、この不揃いのfastqデータが本当にmappingできないのかみてみたい。使ったのはSTARとbowtie2である。

• STAR: STAR_2.4.0j
• bowtie: version 2.2.5

まずはデータを準備。たまたま利用しているpublic dataを使いたい。元のデータはpaired-endのデータを合体させたfastqであったので、わざわざ２つのfastqに分けなくてはならないが省略。あとは、5'を15basesトリムしているデータを用意した。これを適当にquality filteringしてみよう。

qc.q25.p80.sh

fastq_quality_filter -Q33 -q 25 -p 80 -i SRR1158259.f.trm.fastq -o SRR1158259.f.trm.q25p80.fastq
fastq_quality_filter -Q33 -q 25 -p 80 -i SRR1158259.r.trm.fastq -o SRR1158259.r.trm.q25p80.fastq

filtering後の行数を確認してみよう。

ls.sh

wc -l  SRR1158259.f.trm.q25p80.fastq
220663356 SRR1158259.f.trm.q25p80.fastq

wc -l SRR1158259.r.trm.q25p80.fastq
210424212 SRR1158259.r.trm.q25p80.fastq

すでに不揃い。しかしこのままmappingするのは時間がかかるので、テストデータを作ろう。

head -n 100000 SRR1158259.f.trm.q25p80.fastq >test.f.trm.q24p80.fastq
tail -n 100000 SRR1158259.r.trm.q25p80.fastq >test.r.trm.q24p80.fastq

完成した2つのデータは、行数（read数）は同じだが、readのpairは一致しない。さてこのデータをmappingしてみよう。STARから結果を示す。

star.sh

FQ1=test.f.trm.q24p80.fastq
FQ2=test.r.trm.q24p80.fastq

#STARの実行
STAR --genomeDir /home/suimye/genome/hg19  --readFilesIn $FQ1 $FQ2 --runThreadsN 3 

###実行結果
STAR --genomeDir /home/suimye/genome/hg19  --readFilesIn $FQ1 $FQ2
May 30 18:29:31 ..... Started STAR run
May 30 18:32:20 ..... Finished successfully

正常に終了したようだ。同様にbowtile2を使ってみよう。

bt2.sh

INDEX_NAME=/home/suimye/genome/hg19/hg19

#bowtie2の実行
bowtie2 -x $INDEX_NAME -1 test.f.trm.q24p80.fastq -2 test.r.trm.q24p80.fastq -S test.bt.sam

#### 実行結果
25000 reads; of these:
  25000 (100.00%) were paired; of these:
    25000 (100.00%) aligned concordantly 0 times
    0 (0.00%) aligned concordantly exactly 1 time
    0 (0.00%) aligned concordantly >1 times
    ----
    25000 pairs aligned concordantly 0 times; of these:
      9278 (37.11%) aligned discordantly 1 time
    ----
    15722 pairs aligned 0 times concordantly or discordantly; of these:
      31444 mates make up the pairs; of these:
        6266 (19.93%) aligned 0 times
        11935 (37.96%) aligned exactly 1 time
        13243 (42.12%) aligned >1 times
87.47% overall alignment rate

ここまでの結論としては、2つのreadのpairは合っていなくてもread数さえ一致すれば、どちらのalignerも使える。さて、ではread数が異なる場合はどうなるか。正解は、

• STAR: Warningもなく実行できる。
• bowtie2: エラーメッセージとともにmappingできずに終了する。

という結果だった。STARの結果は載せませんが、bowtie2のエラーメッセージについては紹介しよう。

異なるread数のpaired

head -n 100000 SRR1158259.f.trm.q25p80.fastq >test.f.trm.q24p80.fastq
head -n 100004 SRR1158259.r.trm.q25p80.fastq >test.r.trm.q24p80.fastq

bowtie2のエラー

#bt2の実行
bowtie2 -x $INDEX_NAME -1 test.f.trm.q24p80.fastq -2 test.r.trm.q24p80.fastq -S test.bt.sam

#エラーメッセージ  
Error, fewer reads in file specified with -1 than in file specified with -2
terminate called after throwing an instance of 'int'
(ERR): bowtie2-align died with signal 6 (ABRT) 
-bash-4.1$