About TopHat accepted.bam

TopHat 의 align 결과 stats 를 뽑기 위해 samtools flagstat을 사용하려 한다. 관련된 내용을 정리한다.

http://www.biostars.org/p/12475/  참조
위 URL이 조금 설명이 잘못된 부분이 있어 아래와 같이 수정한다

<example>

아래 내용은 tophat에서 생성된 bam 파일을 samtools flagstat을 한 결과이다.

40578140 + 0 in total (QC-passed reads + QC-failed reads)  #bam파일에 recording 되어 있는 read 수

0 + 0 duplicates 

40578140 + 0 mapped (100.00%:nan%)  #tophat 자체는 mapping 된 read의 정보만 bam 파일을 만들기 때문에 bam 파일에 있는 read는 모두 mapping 된 것으로 100%가 나옴, 0x4가 없는 read

40578140 + 0 paired in sequencing  #paired-end sequencing, 0x1 을 갖는 read

20716091 + 0 read1  #bam 파일에서 read1의 갯수 #0x40을 갖은 read

19862049 + 0 read2  #bam 파일에서 read2의 갯수 #0x80을 갖은 read

14778858 + 0 properly paired (36.42%:nan%)  #read1과 read2가 같은 chromosome 에 적절한 방향성을 같고 align 이 제대로 된 경우, 0x2를 갖고 있는 read

31546712 + 0 with itself and mate mapped #pair이고 mate도 mapping이 됐지만 proper 하게 mapping이 안된경우, 0x1을 가지고 있으나 0x2는 없고 0x8도 없는 경우

9031428 + 0 singletons (22.26%:nan%) #read itself는 mapping이 됐지만 pair가 mapping이 안된 경우, 0x8을 가지고 있는 read

0 + 0 with mate mapped to a different chr

0 + 0 with mate mapped to a different chr (mapQ>=5)

결과적으로..
tophat의 align 결과인 accepted_hits.bam 파일은 multiple alignment option를 허용했을 때, 하나의 read에 대해 여러개의 record를 남긴다. samtools flagstat은 read의 flag 값만을 가지고 stat을 뽑는걸로 보이는데 그렇기 때문에 read가 중복 count 된다. 
이해가 안되는 것은 마직막 stats, "with mate mapped to a different chr" 은 왜  count 하지 않는냐는 것.

<SAM format flag>

16진수의 3자리의 수의 10진법으로 변환한 것이 flag이므로 3번째 자리의 숫자가 있는지 확인하기 위해서는 16^2 = 256 의 몫을 2번째 자리의 숫자를 확인하기 위해서는 16의 몫을 그리고 나머지가 1번째 자리의 수.

예)

163 = 16*10 + 3

에서 10 은 0x80 + 0x20 을 의미하고 3 은 0x1 + 0x2 를 의미

또한 각 bit 의 의미는 아래와 같음.

0x1  read paired

0x2  read mapped in proper pair

0x4  read unmapped

0x8  mate unmapped

0x10  read reverse strand

0x20  mate reverse strand

0x40  first in pair

0x80  second in pair

0x100  not primary alignment

0x200  read fails platform/vendor quality checks

0x400  read is PCR or optical duplicate

<Illumina sequence identifiers>

http://en.wikipedia.org/wiki/FASTQ_format 참조

한 샘플에 대해 여러 lane을 걸게 되면 하나의 fastq안에 서로 다른 lane number를 갖을 수 있음.

samtools pileup format

과거 포스팅을 보면..
http://graphy2111.blogspot.kr/2011/07/sanger-fastq-file-format-for-sequences.html



<conversion illumina 1.3+ quality ascii to phred score>

illumina 1.3+ fastq의 ascii를 quality로 바꿔주기 위해서는 
python의 경우)

• ord('문자') - 64   # 64가 offset이기 때문에

<pileup format>
http://samtools.sourceforge.net/pileup.shtml

예)

seq1 272 T 24  ,.$.....,,.,.,...,,,.,..^+. <<<+;<<<<<<<<<<<=<;<;7<&
<chrmosome> <1-based coordinate> <reference base> <# of reads> <read bases> <base qualities>

<base qualities> :: 위 illumina phred score 확인.

<read bases> :: 

"." :: match to ref base on the forward strand
"," :: match to ref base on the reverse strand
"ACGTN" :: mismatch on the forward
"acgtn" :: mismatch on the reverse
"\+[0-9]+[ACGTNacgtn]+" :: 그 위치와 다음 위치의 base 사이에 insertion이 있음을 의미
"\-[0-9]+[ACGTNacgtn]+" :: 그 위치와 다음 위치의 base 사이에 deletion이 있음을 의미

"^" :: start of a read segment, "^" 뒤의 문자는 mapping quality를 나타내는 ascii code (-33 필요)
"$" :: end of a read 


중요한건 insertion, deletion, start, end를 나타내는 기호 다음에는 반드시 ".", ",", "ACGTN", "acgtn" 가 나타난다. 

예)


seq2 156 A 11  .$......+2AG.+2AG.+2AGGG    <975;:<<<<< 
(.)($.)(.)(.)(.)(.)(.)(+2AG.)(+2AG.)(+2AGG)(G)  를 의미함.



위에서 보면 +n 다음 나오는 n개의 BASE(혹은 character)는 insert 된 base를 나타내고 그 다음 base는 그 위치에서는 base에 대한 정보. 곧 +2AG. 이라 함은 2개의 base AG가 insertion 되었고 .은 그 위치에서 base가 reference와 동일하다를 의미 

                     |  <===== 이 위치에서의 pileup을 생각해보면

reference :: AGTTCG--G

read      :: AGTTCGAGG

이는 +2AG. 으로 표현됨