Align Followed by Samtools

Posted on Edited on

出发点：
1.对原序列一段区域感兴趣，将这一小段序列取出来保存为fasta文件作为reference
2.fastq文件align到reference，大部分是map不上的，只要map上的
3.比对结果可视化

用bowtie2进行align，用samtools去掉比对不上的

1
2
3
4
5
6
7
8
9
10
bowtie2-build test2.fasta index
bowtie2 -p 8 -x index -1 H11_01_1.fastq -2 H11_01_2.fastq -S out.sam
# 用samtools将sam文件转为bam文件（二进制）
samtools view -bS out.sam > tmp.bam
# 用samtools去掉比对不上的
samtools view -F 0x04 -b tmp.bam > outtmp.bam
# 以上两个命令可以合并
samtools view -bS -F 0x04 out.sam > outtmp.bam
# sort，按位置排序，且能进一步压缩存储
samtools sort outtmp.bam > tmp.sort.bam

主要是因为，我感兴趣的序列是原来assembly出来的fasta文件的很小一部分，所以比对上的比例是很小的，Bowtie2会把所有的reads都输出包括那些unmapped and redundant data，根据我的文件这些冗余数据占了99.97%. 而我想看的只是0.03%比对上的部分。

这个地方还可以有其它方法，在比对输出的时候就输出比对上的reads，可以看bbmap.

  • 关于samtools

SAM(Sequence Alignment/Map) 格式是用于存储大核苷酸序列比对large nucleotide sequence alignments的通用格式。

  • 具有足够的灵活性来存储由各种alignment programs生成的所有alignment信息；
  • 足够简单，可以由alignment programs轻松生成或从现有alignment format转换而来；
  • 文件尺寸紧凑；
  • 大多数操作都可以 work on a stream ，而无需将整个文件加载到内存中；
  • 允许通过基因组位置对文件进行索引，以有效地检索与位点对齐的所有读取。

SAM Tools提供了各种实用程序，用于处理SAM格式的比对，包括按位置格式排序，合并，建立索引和生成alignments in a per-position format.

比对结果+可视化

vcf文件

1
2
3
4
5
# bcf->vcf
samtools mpileup -gSDf test2.fasta tmp.sort.bam > var.raw.bcf
bcftools view var.raw.bcf|less -S
bcftools view var.raw.bcf > var.raw.vcf
less -SN var.raw.vcf

在上面的命令行中，samtools收集输入的BAM中的摘要信息，计算给定每种可能基因型的数据的可能性，并以BCF格式存储可能性。It does not call variants.

Bcftools applies the prior and does the actual calling. 它还可以连接BCF文件，为快速随机访问索引BCF，并将BCF转换为VCF。 此外，bcftools可以在某些VCF上运行（例如，从带有GL标签的VCF调用SNP），但不能在所有VCF上使用； VCF到BCF的转换目前也不起作用。

更多见
http://samtools.sourceforge.net/mpileup.shtml

pileup文件

1
2
3
# pileup格式文件
samtools mpileup -f test2.fasta -Q 15 -q 20 tmp.sort.bam -o tmp.pileup
less -S tmp.pileup

关于参数-Q和-q，在mpileup里-Q控制碱基的质量,-q控制比对的质量。也就是看到的pileup文件实际上是经过筛选的文件。

1
2
3
4
5
6
7
8
9
10
samtools depth |& grep quality
   -q <int>            base quality threshold [0]
   -Q <int>            mapping quality threshold [0]

samtools mpileup |& grep smaller
  -q, --min-MQ INT        skip alignments with mapQ smaller than INT [0]
  -Q, --min-BQ INT        skip bases with baseQ/BAQ smaller than INT [13]

samtools view | grep qua
  -q INT   only include reads with mapping quality >= INT [0]

The pileup format有多种变体。 SAMtools的默认输出如下所示：

1
2
3
4
5
6
7
8
seq1 272 T 24  ,.$.....,,.,.,...,,,.,..^+. <<<+;<<<<<<<<<<<=<;<;7<&
seq1 273 T 23  ,.....,,.,.,...,,,.,..A <<<;<<<<<<<<<3<=<<<;<<+
seq1 274 T 23  ,.$....,,.,.,...,,,.,...    7<7;<;<<<<<<<<<=<;<;<<6
seq1 275 A 23  ,$....,,.,.,...,,,.,...^l.  <+;9*<<<<<<<<<=<<:;<<<<
seq1 276 G 22  ...T,,.,.,...,,,.,....  33;+<<7=7<<7<&<<1;<<6<
seq1 277 T 22  ....,,.,.,.C.,,,.,..G.  +7<;<<<<<<<&<=<<:;<<&<
seq1 278 G 23  ....,,.,.,...,,,.,....^k.   %38*<<;<7<<7<=<<<;<<<<<
seq1 279 C 23  A..T,,.,.,...,,,.,..... ;75&<<<<<<<<<=<<<9<<:<<

默认输出的六列分别代表：

  • 染色体 chromosome
  • 碱基坐标 1-based coordinate
  • 参考碱基 reference base
  • 覆盖该位点的reads数 the number of reads covering the site
  • read碱基 read bases
  • 碱基质量 base qualities

对于read碱基列，会按照以下顺序进行标注：

  • 如果这是read覆盖的起始位点，会用^符号并跟着比对质量
  • 一个字符表示碱基以及比对上的正反链
正链 反链 含义
. dot , comma 碱基比对上参考序列该位点
ACGTN acgtn 碱基发生突变mismatch
> < Reference skip (due to CIGAR “N”)
* */# Deletion of the reference base (CIGAR “D”)

删失正反向都会用\*除非使用--reverse-del参数，此时反链比对显示#。

  • 模式“\ + [0-9] + [ACGTNacgtn*#] +”表示此参考位置和下一个参考位置之间有插入。插入的长度由模式中的整数给出，后跟插入的序列。
  • 模式“- [0-9] + [ACGTNacgtn*#] +”表示此参考位置和下一个参考位置之间有删失。
  • 如果这是read覆盖的终止位点，跟上$符号。

更多见
http://samtools.sourceforge.net/pileup.shtml
and here

IGV可视化

为bam文件建立索引，生成一个.bai文件，两者放入同一文件夹下，它根据文件名自动和.bam关联。

1
samtools index tmp.sort.bam

打开IGV，File->load from file将bam导入IGV中；Genomes->load from genome导入参考序列fasta文件。

samtools tview可视化

1
samtools tview tmp.sort.bam test2.fasta

.表示和参考序列一致，,表示比对到参考序列互补链，字母代表mismatch。颜色标注比对质量或碱基质量，30-40白色，20-30黄色，10-20绿色，0-10蓝色。
按?显示帮助菜单。

  • H(左）J（上）K（下）L（右）移动显示界面。大写字母移动快，小写字母移动慢。或者上下左右键。
  • 点号.切换显示碱基和点号，r切换显示read name.
  • g快速定位，如我们定位到NODE_10_length_213724_cov_21.069325:100位置，按下enter就可以到达该位置。
  • q退出。