0. 结论

在使用TopHat2匹配双端测序结果后,建议根据成对reads的map基因组位置唯一、方向正确和距离合适的标准,筛选得到的匹配结果。比如,TopHat2可能生成accepted_hits.bam文件,处理方法如下:

Filter TopHat Outputs
1
2
3
4
5
6
7
8
9
10
11
# 首先查看bam文件头部有多少行
$ samtools view -H accepted_hits.bam | wc -l
86

# 筛选成对且成功map到基因组唯一位置的reads,按照上一条输出结果,调整“NR <= 86”
$ samtools view -h accepted_hits.bam | \
    awk '{if (NR <= 86) print $0}; {if($5 == 50 && ($2 == 163 || $2 == 147 || $2 == 83 || $2 == 99)) print $0}' | \
    samtools view -b - > accepted_filtered.bam

$ samtools view accepted_filtered.bam | wc -l
79143942

1. TopHat输出sam文件的第五列

TopHat文档没有解释其输出bam文件(比如accepted_hits.bam)的第五列的意义。按照Bowtie2输出结果来看,是表示映射质量指标Mapping Quality scores(MAPQ),具体计算参考公式$\eqref{eq:1}$。

MAPQ值越大,表示对应的read的alignment质量越高。然而,在TopHat输出结果中,MAPQ所代表意义略有不同。

“MAPQ” distribution from TopHat2 accepted mapping reads
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
# 查看accepted_hits.bam文件的MAPQ数值,并统计出现频数
$ samtools view accepted_hits.bam | awk '{print $5}' | sort --parallel=4 -n | uniq -c
5057430 0
3117731 1
8058500 3
93044727 50

# 查看前100位MAPQ数值和NH:i:n分布
$ samtools view accepted_hits.bam | awk '{print $5}' | head -100 | sort --parallel=4 -n | uniq -c
35 0
42 1
22 3
1 50

$ samtools view accepted_hits.bam | awk '{print $20}' | head -100 | sort --parallel=4 -n | uniq -c
1 NH:i:1
22 NH:i:2
3 NH:i:20
13 NH:i:3
29 NH:i:4
5 NH:i:5
21 NH:i:6
1 NH:i:7
2 NH:i:8
3 NH:i:9

# 查看前200位MAPQ数值和NH:i:n分布
$ samtools view accepted_hits.bam | awk '{print $5}' | head -200 | sort --parallel=4 -n | uniq -c
60 0
43 1
67 3
30 50

$ samtools view accepted_hits.bam | awk '{print $20}' | head -200 | sort --parallel=4 -n | uniq -c
30 NH:i:1
8 NH:i:12
7 NH:i:14
67 NH:i:2
3 NH:i:20
13 NH:i:3
30 NH:i:4
5 NH:i:5
31 NH:i:6
1 NH:i:7
2 NH:i:8
3 NH:i:9

首先,我们可以看到TopHat输出的MAPQ只有四个数值,分别为50310。根据sam文件标准NH:i:n表示含有查询序列的alignment的数量。因此,通过上述前100位和前200位分析可以发现,MAPQ并不是按照公式$\eqref{eq:1}$计算,而有可能是以下关系:


MAPQ (tophat) Tag 描述
50 NH:i:1 map至唯一位置
3 NH:i:2 map至2个位置
1 NH:i:3/NH:i:4 map至3个或4个位置
0 NH:i:n (n > 4) map到多余4个位置

展示一个NH:i:1的例子,注意Illumina双端测序平台fr-unstranded

Two pairs of unique mapped reads from Illumina HiSeq2000
1
2
3
4
5
HISEQ2000-02:436:C2PG3ACXX:3:2313:10972:95322   163     chr1    637224  50      100M    =       637339  215     AAATGATCTGTACAATAACCCCCTGTGACACCAGTCTACCTATGTAACAAATGCCCCTAAACTTAAAATAAAAGTTAAAAAAAAAGAAAATTAAAATCTC  <@@BDABBDFBCDGEGHIGIIGIABAFHBGDFGGGHIIIFIEIGGGGIIIFFDHIGIIIIIIIICEHIIIIIIHEHHDHFFCBCCB@BCAACCCCCECCC    AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:100     YT:Z:UU  NH:i:1
HISEQ2000-02:436:C2PG3ACXX:3:2313:10972:95322   83      chr1    637339  50      100M    =       637224  -215    GTAATATGAAAAACACAAATCTTTCATTCATTCCTTTCAACTGATGAGGAAAATGAGGCATCGGGAGTTAGTAAAAGTCCACATTGAGATATGAGACCCA  CCADDDCCCCCCCDEEEECAEHEEGGIIHGFAAGGGHEF=IGGGIIGGHGCGIEIIIGIIIIIIIFIHIIIIIIGIIHFEAIGGIIGFFDFHDDDAD@@@    AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:100     YT:Z:UU  NH:i:1

HISEQ2000-06:325:C2RC0ACXX:5:2205:6961:88285    99      chr1    643662  50      100M    =       643707  145     CCTATCAAAATCTTAGCATTCCTCTTAGCCCTCAACAAAGCATTTCTAAAATGTGTATAGAAGACCAAAGGGCCAAAAGAGTCAACTTCTGAAGAAGCGC  CCCFFFFFHHHHHJJJJIJJHJJJJJJJIIJJJJJJJJJJJJIJJJJJIIIJJIIIJJJJJJJJJIJJJJIJJJIHHHHFFEFFFEEEEEEEDDCDDCDD    AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:100     YT:Z:UU  NH:i:1
HISEQ2000-06:325:C2RC0ACXX:5:2205:6961:88285    147     chr1    643707  50      100M    =       643662  -145    CTAAAATGTGTATAGAAGACCAAAGGGCCAAAAGAGTCAACTTCTGAAGAAGCGCAAAAAGAAAGTTGAGGAAATCTTAAAACATGTTATTGAGCTTAAA  CEEEDDDDEDFEEEEEEEBFFFFFHHHHEJJJJJJJJIJJJJGJIIIIHFJJJIIJJJJJIJJJJJJJJJJJJGHJJJJIJJJJJJJGHHHHFFFFFCCC    AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:100     YT:Z:UU  NH:i:1

2. sam文件的第二列

sam文件中的第二列提供了具体的map情况,下列表格摘自sam文件标准,sam/bam文件中第二列各种条件求和十进制标识,快速解释第二列的bitwise FLAG


Decimal Hexadecimal Description
1 0x1 template having multiple segments in sequencing
2 0x2 each segment properly aligned according to the aligner
4 0x4 segment unmapped
8 0x8 next segment in the template unmapped
16 0x10 SEQ being reverse complemented
32 0x20 SEQ of the next segment in the template being reverse complemented
64 0x40 the first segment in the template
128 0x80 the last segment in the template
256 0x100 secondary alignment
512 0x200 not passing quality controls
1024 0x400 PCR or optical duplicate
2048 0x800 supplementary alignment

在map完成双端测序序列中,我们感兴趣的是一对reads都正确align到基因组上,而且方相匹配又距离合适。符合这样条件的reads,对应的第二列数值为99、147、83和163,具体图示参考Directional RNA-seq— Part 1: SAM file flags。下面表格解释四个数值的具体意义,其中1标识双端测序,2表示一对reads正确地map到基因组合适位置,表格中着重陈述643212816


Flag Composition Explanation
99 64+32+2+1 一对引物中第一个map到基因组正义链;第二个反方向map到基因组正义链
147 128+16+2+1 一对引物中第二个反方向map基因组正义链;第一个map到基因组正义链
83 64+16+2+1 一对引物中第一个map到基因组反义链;第二个反方向map到基因组反义链
163 128+32+2+1 一对引物中第二个反方向map基因组正义链;第一个map到基因组正义链

之后,我们需要筛选含有这些flags的reads。由于我们通常需要操作bam文件,也希望输出是bam文件,中间过程不希望再重新生成sam文件。那么,就需要结合使用awk进行筛选,具体方法见本篇文章开头所示。当然,如果是只是查看,可以使用下面例子中的 samtools view -f 0x2

Count number of unique pair-mapped alignments
1
2
3
4
5
6
7
8
9
# map到基因组上唯一位置的reads数目
$ samtools view -q 50 accepted_hits.bam | wc -l
93044727
# 成对reads都map到基因组对应位置的reads数目
$ samtools view -f 0x2 accepted_hits.bam | wc -l
88793640
# 成对且唯一mapped的reads数目
$ samtools view -q 50 -f 0x2 accepted_hits.bam | wc -l
79143942

3. 操作TopHat输出的文件命令集锦

Useful bash for bam files from TopHat
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# samtools的view -c命令,其实就是输出sam文件有多少行
$ samtools view accepted_hits.bam | wc -l
109278388
$ samtools view -c accepted_hits.bam
109278388

# 查看bam文件中mapped的reads长度分布
# 第二种方法是运行FastQC,输出结果中也有显示
samtools view accepted_hits.bam | awk '{print length($10)}' | sort -n | uniq -c

# 查看bed文件前几行
$ head junction.bed
# 统计bed文件有多少行,需要去除第一行注释
# 以下两种方式相同,但不够完美
$ wc -l junction.bed
220648 junctions.bed
$ awk 'END {print NR}' junctions.bed
220648

4. DNA链和mRNA链的称呼总结

双链DNA和单链mRNA,对每条链都有特定的称呼。总结如下:

DNA/RNA strands
1
2
3
4
5
6
3'~~~~~UCUGAU~~~~~ 5' mRNA的对应基因信息在reverse strand

5'-----AGACTA----------ATTGTT----- 3'
3'-----TCTGAT----------TAACAA----- 5'

                5'~~~~~AUUGUU~~~~~ 3' mRNA的对应基因信息在forward strand

对于一条双链DNA,称呼列表如下:

方向a 名称1 名称2
从左至右 forward plus
从右至左 reverse minus

a:方向是指5’至3’的阅读方向,用于区分两条DNA链条


对于一条RNA链,其对应的双链DNA称呼如下:

mRNA方向a 名称1 名称2 名称3 名称4 名称5
同向 coding nontemplate sense positive +
反向 template noncoding antisense negative -

a:mRNA方向是指5’至3’。


参考网址

更新记录

2015年5月25日

Comments