Bioinformatic Data Skills 学习专题(7) 比对数据的小实操

CHAPTER 12: Working with Alignment Data

前言

真是好久不见啦！在本该撰写专题内容的整个八月和九月前半部分，小编经历了离职、旅游以及西湖到复旦的开学典礼、寝室大扫除搬迁、生病养病、选课抢课等等诸多磨难终于在中秋节来临之前又完成了一章的内容（“一把辛酸泪”TnT），目测即将进入持续产出和摄入的状态。

这讲内容主要涉及常用常见的一些测序相关知识：测序数据比对结果的解读和量化分析，其实也是分析NGS数据的基础步骤了。

由于长时间没有更新此专题，先回顾一下这个专题已经讲了什么吧（从五月开始到现在了嗯~ o(￣▽￣)o）：

1. 常见的linux操作：指令、文本（正则表达式）、git协作

Bioinformatic Data Skills 学习专题(4) git

Bioinformatic Data Skills 学习专题(3) 文本操作

Bioinformatic Data Skills 学习专题(2)：一些基础的shell脚本和服务器指令

2. 常见生物数据，genomic range

Bioinformatic Data Skills 学习专题(6) Genomic Range之三

Bioinformatic Data Skills 学习专题(6) Genomic Range之二

Bioinformatic Data Skills 学习专题(6) Genomic Range之一

Bioinformatic Data Skills 学习专题(5) Bioinformatics Data

本章主要内容

• SAM/BAM的数据格式组成：header/section/flags/CIGAR/Mapping quality

• SAM/BAM文件处理的命令行工具: samtools

• 用IGV可视化

• 用Pysam自己做SAM/BAM处理工具（难度较大，日后再讲）

这些内容集中回答了如下几个问题……

Q1: 比对后的文件格式：多了什么？

从上一章我们讲到的原始测序数据用软件进比对和回帖（主流比对软件如bowtie, hisat，以及最近流行的超快转录组比对软件salmon和kallisto），接下来就要解释比对结果，它是以BAM/SAM文件格式呈现的，主要包括：

1. 头信息（Header），以@开头

1. @SQ, reference sequences，通常是染色体名，主要字段为SN(sequence name )，后面跟着染色体号；LN是则是sequnce length

2. @RG, 样本名和组名(read group and sample metadata)， ID必须是唯一的，作者建议the name of the sequencing run and lane，可以用于后期区分batch effect；同时SM号码可以用于存储样本的metadata信息(individual, treatment group, replicate)；PL 这个域可以存储测序平台信息，比如ILLUMINA，PACBIO等等。

3. @PG，存储比对软件的信息，ID是必须的；VN为版本号；CL则是command line的缩写。通常比对软件会自动加上这些信息

4. 第一行的比对信息不以@开头

2. 每个reads的具体信息（Alignmet Section）

1. 1    QNAME   Query template/pair NAME

2. 2    FLAG    bitwise FLAG

3. 3    RNAME   Reference sequence NAME

4. 4    POS 1-based leftmost POSition/coordinate of clipped sequence 位置信息

5. 5    MAPQ    MAPping Quality (Phred-scaled) 比对碱基质量，详见上一节

6. 6    CIGAR   extended CIGAR string, matching bases, insertions/deletions, clipping

7. 7    MRNM    Mate Reference sequence NaMe (`=' if same as RNAME)  MPOS   1-based Mate POSistion,下一个信息？

8. 8 TLEN    inferred Template LENgth (insert size)， template length for paired-end read

9. 9 SEQ    query SEQuence on the same strand as the reference

10. 10 QUAL    query QUALity (ASCII-33 gives the Phred base quality)

我们需要的信息主要从下面三个字段获得：

• bitwise flag 用16进制的数据表述read的信息， samtools flags可以用于显示各类flag代表的具体信息，当然SAM格式说明里又更详细的，可以参考官网说明，需要特定reads的时候就是通过这个条目进行筛选的

• CIGAR string: matches/mismatches, insertions, deletions, soft or hard clipped, and so on.根据描述，就是看看reads是否有前面这些特征

• Mapping qualities: Q = -10 log10P(incorrect mapping position). 碱基质量，于30就可以了。

相关代码

看SAM文件的header信息

1. samtools view -H celegans.sam

2. @SQ SN:I LN:15072434

3. @SQ SN:II LN:15279421

4. @SQ SN:III LN:13783801

看bitwise flag信息

1. $ samtools flags 147

2. 0x93 147 PAIRED,PROPER_PAIR,REVERSE,READ2

3. $ samtools flags 0x93

4. 0x93 147 PAIRED,PROPER_PAIR,REVERSE,READ2

Q2: 了解了数据格式。处理这些比对好的文件？换一种说法就是，怎么分析SAM文件以得到所需要的特定reads（正确的比对结果）？

1. 格式转换

可以参考samtools的详细说明书(详见http://www.htslib.org/)sam到bam文件格式转化： samtools view-b celegans.sam>celegans_copy.bam包含header信息的转换： samtools view-h celegans.bam>celegans_copy.sam对map好的文件进行排序： samtools sort celegans_unsorted.bam celegans_sorted由于体积小，一般sort好的bam文件被用于后续的分析，这样我们用 samtools view的速度就会快很多，一般情况下，sort好的文件分析都需要以fai结尾的index文件： samtools index celegans_sorted.bam

2. reads挑选实例1：使用 samtools view找特定区域上的比对情况

作者再github项目上分享了1000 Genomes Project Consortium的数据来实战，具体可以到https://github.com/vsbuffalo/bds-files/tree/master/chapter-11-alignment 下载

1. samtools index NA12891_CEU_sample.bam

2. # 看一下 chromosome 1, 215,906,469-215,906,652上的比对reads情况

3. samtools view NA12891_CEU_sample.bam 1:215906469-215906652 | head -n 3

4. SRR003212.5855757 147 1 215906433 60 33S43M = 215906402 [...]

5. SRR003206.18432241 163 1 215906434 60 43M8S = 215906468 [...]

6. SRR014595.5642583 16 1 215906435 37 8S43M * 0 [...]

7. # 我们可以把获得的信息重定向到别的文本

8. $ samtools view -b NA12891_CEU_sample.bam 1:215906469-215906652 >

9. USH2A_sample_alns.bam

10. # 如果你想获得多个区域的比对情况，就需要用这些区域的BED文件，同时加上参数 -L

11. $ samtools view -L USH2A_exons.bed NA12891_CEU_sample.bam | head -n 3

12. SRR003214.11652876 163 1 215796180 60 76M = 215796224 92 [...]

13. SRR010927.6484300 163 1 215796188 60 51M = 215796213 76 [...]

14. SRR005667.2049283 163 1 215796190 60 51M = 215796340 201 [...]

3. reads挑选实例2：使用 samtools view-f筛选特定的比对结果

比如我们想看看那些reads没有比对上哪~ samtools flags unmap告诉我们：0x4 4 UNMAP，4就是没有map上的bitflag，接下来就要用view -f的参数筛选出这些reads:

1. $ samtools view -f 4 NA12891_CEU_sample.bam | head -n 3

2. SRR003208.1496374 69 1 215623168 0 35M16S = 215623168 0 [...]

3. SRR002141.16953736 181 1 215623191 0 40M11S = 215623191 0 [...]

4. SRR002143.2512308 181 1 215623216 0 * = 215623216 0 [...]

参数 view-F的含义刚好相反，就是map好的reads信息了

1. samtools view -F 4 NA12891_CEU_sample.bam | head -n 3

2. SRR005672.8895 99 1 215622850 60 51M = 215623041 227 [...]

3. SRR010927.10846964 163 1 215622860 60 35M16S = 215622892 83 [...]

4. SRR005674.4317449 99 1 215622863 37 51M = 215622987 175 [...]

可以注意到这些bitflag有很多别的数字，想看69含义，如下：

1. $ samtools flags 69

2. 0x45 69 PAIRED,UNMAP,READ1

同时我们也可以发过来用 samtools flags

1. $ samtools flags READ1,PROPER_PAIR

2. 0x42 66 PROPER_PAIR,READ1

最后上一个-f -F的组合使用，需要paired且mapped的reads（有没有像逻辑学的感觉）

1. 0x1 1 PAIRED

2. $ samtools flags unmap,proper_pair

3. 0x6 6 PROPER_PAIR,UNMAP

4. $ samtools view -F 6 -f 1 NA12891_CEU_sample.bam | head -n 3

5. SRR003208.1496374 137 1 215623168 0 35M16S = 215623168 [...]

6. ERR002297.5178166 177 1 215623174 0 36M = 215582813 [...]

7. SRR002141.16953736 121 1 215623191 0 7S44M = 215623191 [...]`

不能盲目使用，需要检验：

1. $ samtools view -F 6 NA12891_CEU_sample.bam | wc -l # total mapped and paired

2. 233628

3. $ samtools view -F 7 NA12891_CEU_sample.bam | wc -l # total mapped, paired,

4. 201101 # proper paired

5. $ samtools view -F 6 -f 1 NA12891_CEU_sample.bam | wc -l # total mapped, paired,

6. 32527 # and not proper paired

7. $ echo '201101 + 32527' | bc

8. 233628

Q3: 用samtools过滤得到正确的比对结果，然后？

直接看fasta文件和bam/sam文件的关系可以用samtools的功能实现：

1. samtools tview -p 1:215906469-215906652 NA12891_CEU_sample.bam \

2. human_g1k_v37.fasta

到了这一步，我们已经能够保证获得的结果都是通过质量控制的，所以技术上的问题已经解决，接下来就要检验数据在生物学意义上是不是靠谱的。比如看看reads的分布是不是合理，一般都需要结合已有的生物学背景知识来进行检验，比如ChIPseq的话就需要看看input是不是均匀分布，IP的结果是不是真的蛋白富集在特定位点，甲基化数据就要看一些特异区域比如CpG岛是不是高甲基化的，重复序列是不是低甲基化等等。通常我们需要把这些结果可视化再做检验，所以数据量化必不可少，大致过程如下：

• 量化步骤：通过samtools或者bedtools把reads pileup，获得压缩的bedgraph/wiggle文件

• 可视化软件：IGV

step1. 量化步骤和脚本

• switch1: 这里墙裂推荐一个工具，deeptools，bam一键量化到bw格式，用于igv的可视化！

• switch2: 用samtools/bedtools pileup的功能获得bedgraph，再用UCSC上的wig2bigwig进行格式转换，也是用于浏览器的可视化

step2. IGV可视化的一些关键点

我再六月初的时候写了一篇关于所需要基因组的IGV参考格式制作，大家可以参考该文章: 何快速有效地将新拿到的参考基因组在IGV里可视化

讲到这里，本章的主要内容差不多结束了，首先介绍了比对数据SAM文件格式的信息组成，然后讲了处理和分析SAM文件的命令行工具samtools。那可不可以不用工具自己分析呢？大拿们当然早已推出各种包了，在此我仅仅摘录一些书中部分（这里摘录的目的其实是今后打算推出单独章节用于此类工具介绍），有兴趣的可以先结合本书学习。

补充部分可能会下期或者最后补讲：定制化的文件分析（Pysam）和自己编程创建SAM文件