tophat mapping statistics와 관련한 수수께끼

tophat-cufflinks pipeline을 이용하여 RNA-seq data에 대한 분석을 시도하였다. 총 여섯 개의 Illumina data file을 매핑하였는데, 이 중에서 3개 결과가 좀 이상하다. 예를 들어 3400만개의 read가 투입되었는데 mapped read의 수는 이보다 많은 5400만개라는 결과가 나온 것이다.

상식적으로 이해가 되지 않는다. fastq file 안에 들어 있는 read의 수를 아무리 헤아려 봐도 34,165,827개이다. tophat 결과 디렉토리에 들어 있는 정보 파일도 거짓말을 하지 않는다.

$ cat prep_reads.info
min_read_len=40
max_read_len=40
reads_in =34165827
reads_out=34129283

그런데,

$ samtools flagstat accepted_hits.bam
54565324 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
54565324 + 0 mapped (100.00%:-nan%)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (-nan%:-nan%)
0 + 0 with itself and mate mapped
0 + 0 singletons (-nan%:-nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

이건 도무지 말이 되지 않는다. spliced alignment를 했으므로, 혹시 하나의 read가 여러 곳에 나뉘어 정렬한 경우 이것이 SAM file 안에서 별도의 line으로 취급되나 싶어서 다시 한 번 점검을 해 보았다. read name을 추출하여 만일 중복이 있다면 이를 제거하고 그 수를 센다.

$ samtools view -F 0x004 accepted_hits.bam | awk '{print $1}' | uniq | wc -l
54565305

으헉! 미치겠군.

2019년 8월 6일 업데이트

당시에는 multiple mapping에 대해서 별로 진지하게 생각해 보지 않았었다. 부끄러운 일이다!