短读长序列比对软件 Bowtie2 的安装与使用

菌心说 2021-10-17

展开全文

Bowtie2 是一款经典的短读长序列（ 50-100 bp，最多可到1000 bp ) 比对软件，支持单端测序（unpaired) 和双端测序的比对。支持全局比对（end-to-end align ) 和局部比对( local align )。

Bowie2 安装

Conda 安装 bowtie2（需要先装好Anaconda）

终端输入命令即可完成安装

   conda install bowtie2

官网下载安装
官网：http://bowtie-bio./bowtie2/index.shtml

Github: https://github.com/BenLangmead/bowtie2/releases
下载对应平台的安装包安装即可。

源码编译安装

官方说明 Building from source

至少需要先装好 zlib 库

输入如下命令
```
cd ~
git clone https://github.com/BenLangmead/bowtie2
cd bowtie2
make 
make install 
```
如果报错，可能有依赖库的头文件没有安装或找到。没有安装先安装缺失的库。安装好了，没有找到对应的库请正确设置其查找路径

Bowtie2 使用

建立索引

Bwotie2 的使用需先对要比对的物种的 genome 建立索引（BWT index）。

如果该物种有建好的index（比如人的），可从官网上下载 http://bowtie-bio./bowtie2/index.shtml。

如果没有，需要自己构建。

先下载准备好该物种的基因组的 fastq 文件。

此处使用 bowtie2 的 example 里的 lambda virus 的 genome 文件 lambda_virus.fa

输入如下命令可以建立索引

# Building a small index
bowtie2-build example/reference/lambda_virus.fa example/index/lambda_virus

建立好索引后会得到6个文件

example/index/lambda_virus.1.bt2, example/index/lambda_virus.2.bt2, example/index/lambda_virus.3.bt2, example/index/lambda_virus.4.bt2,
example/index/lambda_virus.rev.1.bt2, and example/index/lambda_virus.rev.2.bt2.

这6个文件里面存储了基因组的信息，BWT index。

其中 .1.bt2 .2.bt2 文件村存储了 forward strand 的 BWT index.rev.1.bt2 .rev.2.bt2 文件存储了 reverse complement strand 的 BWT index。.3.bt2 .4.bt2 压缩存储了genome序列信息。

后续比对分析只用到这6个文件。

如果建立长的 index，使用如下命令。

# Building a large index
bowtie2-build --large-index example/reference/lambda_virus.fa example/index/lambda_virus

输出文件的后缀为 .bt2l.

比对

比对的过程是将 read 和 genome reference 对齐的过程，如下图所示

  Read:        CTGGCCGATCTCGACTT
               |||  |||||||||| |
  Reference: GACTG--CGATCTCGACATCG

其中竖线表示二者匹配，没有竖线表示有替换 , reference 中 A 被替换成了 T。短横线表示有空位（gap), 即 read 中插入了 GC 序列。

如果是单端测序的reads. -x 指定的建立的index，使用 -U 选项指定 unpaired align mode, 使用 -S 指定输出sam文件

# Aligning unpaired reads
bowtie2 -x example/index/lambda_virus -U example/reads/longreads.fq -S example/reads/longreads_align.sam

如果是双端端测序的 reads. 使用 -1 -2 选项, 分别制定两端测得的 Reads 文件

# Aligning paired reads
bowtie2 -x example/index/lambda_virus -1 example/reads/reads_1.fq -2 example/reads/reads_2.fq -S example/reads/reads_align.sam

比对好之后可使用 samtool 进行后续分析。

算法细节

Bowtie2 使用 Borrow Wheller Transorm 建立索引。对给定的一条 read 将其分成一些小片段（seeds），然后使用 BWT index 来进行 seeds 的精确比对（通过backtrap 的方法，可允许少量的 mismatch（替换））。然后使用动态规划建这些seeds拼接起来（允许有 gap)，看是否能够比对到 reference genome 上。如果比对上了，返回其在genome上的位置信息，和相关信息。如果没有比对上。则报告其结果

人的基因组，注释文件，Bowtie2 建好的index下载与分析

人的基因组下载

https://hgdownload.soe./goldenPath/hg38/bigZips/

 wget https://hgdownload.soe./goldenPath/hg38/bigZips/latest/hg38.fa.gz
 # 解压
 gzip -d hg38.fa.gz

# 统计 chars 
wc hg38.fa
#  65147824 65147824 3322506630 hg38.fa 
#  约 3.3G nucleotides
more hg38.fa
# 按回车向下划动会看到如下图所示，其中 ... 表示省略掉了很多行

>chr1
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
... many lines N 
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
taaccctaaccctaaccctaaccctaaccctaaccctaaccctaacccta
accctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaac
...
ggatttttgccagtctaacaggtgaAGccctggagattcttattagtgat
ttgggctggggcctggccatgtgtatttttttaaatttccactgatgatt
ttgctgcatggccggtgttgagaatgactgCGCAAATTTGCCGGATTTCC
TTTGCTGTTCCTGCATGTAGTTTAAACGAGATTGCCAGCACCGGGTATCA
TTCACCATTTTTCTTTTCGTTAACTTGCCGTCAGCCTTTTCTTTGACCTC
TTCTTTCTGTTCATGTGTATTTGCTGTCTCTTAGCCCAGACTTCCCGTGT
CCTTTCCACCGGGCCTTTGAGAGGTCACAGGGTCTTGATGCTGTGGTCTT
CATCTGCAGGTGTCTGACTTCCAGCAACTGCTGGCCTGTGCCAGGGTGCA
AGCTGAGCACTGGAGTGGAGTTTTCCTGTGGAGAGGAGCCATGCCTAGAG
TGGGATGGGCCATTGTTCATCTTCTGGCCCCTGTTGTCTGCATGTAACTT
AATACCACAACCAGGCATAGGGGAAAGATTGGAGGAAAGATGAGTGAGAG
CATCAACTTCTCTCACAACCTAGGCCAGTAAGTAGTGCTTGTGCTCATCT
CCTTGGCTGTGATACGTGGCCGGCCCTCGCTCCAGCAGCTGGACCCCTAC
CTGCCGTCTGCTGCCATCGGAGCCCAAAGCCGGGCTGTGACTGCTCAGAC
CAGCCGGCTGGAGGGAGGGGCTCAGCAGGTCTGGCTTTGGCCCTGGGAGA
GCAGGTGGAAGATCAGGCAGGCCATCGCTGCCACAGAACCCAGTGGATTG
GCCTAGGTGGGATCTCTGAGCTCAACAAGCCCTCTCTGGGTGGTAGGTGC

注释文件下载

基因注释文件注释了基因组上的功能片段的位置信息。比如基因，外显子等。

下载地址

http://ftp./pub/databases/gencode/Gencode_human/latest_release/

wget http://ftp./pub/databases/gencode/Gencode_human/release_38/gencode.v38.annotation.gtf.gz
# 解压
gzip -d gencode.v38.annotation.gtf.gz
# 查看表头
head gencode.v38.annotation.gtf

##description: evidence-based annotation of the human genome (GRCh38), version 38 (Ensembl 104)
##provider: GENCODE
##contact: gencode-help@
##format: gtf
##date: 2021-03-12
chr1    HAVANA  gene    11869   14409   .   +   .   gene_id 'ENSG00000223972.5'; gene_type 'transcribed_unprocessed_pseudogene'; gene_name 'DDX11L1'; level 2; hgnc_id 'HGNC:37102'; havana_gene 'OTTHUMG00000000961.2';
chr1    HAVANA  transcript  11869   14409   .   +   .   gene_id 'ENSG00000223972.5'; transcript_id 'ENST00000456328.2'; gene_type 'transcribed_unprocessed_pseudogene'; gene_name 'DDX11L1'; transcript_type 'processed_transcript'; transcript_name 'DDX11L1-202'; level 2; transcript_support_level '1'; hgnc_id 'HGNC:37102'; tag 'basic'; havana_gene 'OTTHUMG00000000961.2'; havana_transcript 'OTTHUMT00000362751.1';
chr1    HAVANA  exon    11869   12227   .   +   .   gene_id 'ENSG00000223972.5'; transcript_id 'ENST00000456328.2'; gene_type 'transcribed_unprocessed_pseudogene'; gene_name 'DDX11L1'; transcript_type 'processed_transcript'; transcript_name 'DDX11L1-202'; exon_number 1; exon_id 'ENSE00002234944.1'; level 2; transcript_support_level '1'; hgnc_id 'HGNC:37102'; tag 'basic'; havana_gene 'OTTHUMG00000000961.2'; havana_transcript 'OTTHUMT00000362751.1';

# 查看注释种类
cat gencode.v38.annotation.gtf | awk '{print $3}' | head -n 30000 |  sort | uniq

CDS
UTR
annotation
exon
gene
start_codon
stop_codon
transcript

Bowtie2 index 下载

官网 http://bowtie-bio./bowtie2/index.shtml

 wget --timestamping  https://genome-idx.s3./bt/GRCh38_noalt_as.zip 
 # 解压到 GRCh38_noalt_as
 unzip GRCh38_noalt_as.zip
 # 查看 index 的 6 个文件

小测试

将上面的 hg38.fa chromosome 1 的序列 copy 一段 50 bp 使用 bowtie2 比对

# -x 索引名称（目录/base_name=bowtie2_index/GRCh38_noalt_as/GRCh38_noalt_as) 
# -c 命令行（commande line ）输入序列, 多个序列以英文 ，分隔。
bowtie2 -x bowtie2_index/GRCh38_noalt_as/GRCh38_noalt_as -c ACAGTACTGGCGGATTATAGGGAAACACCCGGAGCATATGCTGTTTGGTC

输出为：

1 reads; of these:
  1 (100.00%) were unpaired; of these:
    0 (0.00%) aligned 0 times
    0 (0.00%) aligned exactly 1 time
    1 (100.00%) aligned >1 times
0   16  chr15   101979467   1   50M *   0   0   GACCAAACAGCATATGCTCCGGGTGTTTCCCTATAATCCGCCAGTACTGT  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII  AS:i:0  XS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:50 YT:Z:UU

可以发现有多个匹配，aligned >1 times。比对到了 chr15, 并且完全比对上了， 50M （M：match)。

因为匹配到了chr15, 说明序列同时在chr1 诃 chr 15 存在，可见基因组中重复还是挺多的。

将上面的 hg38.fa chromosome 1 的序列连续取五行粘在一起，一共 300 bp，使用 bowtie2 比对

bowtie2 -x bowtie2_index/GRCh38_noalt_as/GRCh38_noalt_as -c CATCAACTTCTCTCACAACCTAGGCCAGTAAGTAGTGCTTGTGCTCATCTCCTTGGCTGTGATACGTGGCCGGCCCTCGCTCCAGCAGCTGGACCCCTACCTGCCGTCTGCTGCCATCGGAGCCCAAAGCCGGGCTGTGACTGCTCAGACCAGCCGGCTGGAGGGAGGGGCTCAGCAGGTCTGGCTTTGGCCCTGGGAGAGCAGGTGGAAGATCAGGCAGGCCATCGCTGCCACAGAACCCAGTGGATTGGCCTAGGTGGGATCTCTGAGCTCAACAAGCCCTCTCTGGGTGGTAGGTGC

输出结果为：

1 reads; of these:
  1 (100.00%) were unpaired; of these:
    0 (0.00%) aligned 0 times
    0 (0.00%) aligned exactly 1 time
    1 (100.00%) aligned >1 times
100.00% overall alignment rate
0   0   chr1    12201   1   300M    *   0   0   CATCAACTTCTCTCACAACCTAGGCCAGTAAGTAGTGCTTGTGCTCATCTCCTTGGCTGTGATACGTGGCCGGCCCTCGCTCCAGCAGCTGGACCCCTACCTGCCGTCTGCTGCCATCGGAGCCCAAAGCCGGGCTGTGACTGCTCAGACCAGCCGGCTGGAGGGAGGGGCTCAGCAGGTCTGGCTTTGGCCCTGGGAGAGCAGGTGGAAGATCAGGCAGGCCATCGCTGCCACAGAACCCAGTGGATTGGCCTAGGTGGGATCTCTGAGCTCAACAAGCCCTCTCTGGGTGGTAGGTGC    IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII    AS:i:0  XS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:300    YT:Z:UU

可以发现还是有多个匹配，aligned >1 times。但是得分最高的匹配比对到了 chr1, 并且完全比对上了， 300M （M：match)。

上面我们是将 chromesome 1 上的序列序列开头随机取了几段来比，发现开头序列在多个在基因组中存在多个匹配。

下面我们将基因注释文件 gencode.v38.annotation.gtf 文件里的 gene 'DDX11L1' 从 hg38.fa 里取出来

chr1    HAVANA  gene    11869   14409   .   +   .   gene_id 'ENSG00000223972.5'; gene_type 'transcribed_unprocessed_pseudogene'; gene_name 'DDX11L1'; level 2; hgnc_id 'HGNC:37102'; havana_gene 'OTTHUMG00000000961.2';

 # chr1 HAVANA  gene    11869   14409 计算其在 hg38.fa 行大概位置
 echo $((11869/50+1)) $((14409/50+1+1))
 # 238  290 

 # gene 'DDX11L1' chr1 start 11869  end 14409
 # locate at line 238 - 290 of hg38.fa
 # number lines 290 - 237 = 53
 head -n 290 hg38.fa | tail -n 53

# 基因的序列为
TTTGCTGTTCCTGCATGTAGTTTAAACGAGATTGCCAGCACCGGGTATCA
TTCACCATTTTTCTTTTCGTTAACTTGCCGTCAGCCTTTTCTTTGACCTC
TTCTTTCTGTTCATGTGTATTTGCTGTCTCTTAGCCCAGACTTCCCGTGT
...
GGTAGAACGGAGCAGCTGGTGATGTGTGGGCCCACCGGCCCCAGGCTCCT
GTCTCCCCCCAGGTGTGTGGTGATGCCAGGCATGCCCTTCCCCAGCATCA
GGTCTCCAGAGCTGCAGAAGACGACGGCCGACTTGGATCACACTCTTGTG
AGTGTCCCCAGTGTTGCAGAGGTGAGAGGAGAGTAGACAGTGAGTGGGAG
TGGCGTCGCCCCTAGGGCTCTACGGGGCCGGCGTCTCCTGTCTCCTGGAG
AGGCTTCGATGCCCCTCCACACCCTCTTGATCTTCCCTGTGATGTCATCT
GGAGCCCTGCTGCTTGCGGTGGCCTATAAAGCCTCCTAGTCTGGCTCCAA
GGCCTGGCAGAGTCTTTCCCAGGGAAAGCTACAAGCAGCAAACAGTCTGC
ATGGGTCATCCCCTTCACTCCCAGCTCAGAGCCCAGGCCAGGGGCCCCCA
AGAAAGGCTCTGGTGGAGAACCTGTGCATGAAGGCTGTCAACCAGTCCAT
AGGCAAGCCTGGCTGCCTCCAGCTGGGTCGACAGACAGGGGCTGGAGAAG
GGGAGAAGAGGAAAGTGAGGTTGCCTGCCCTGTCTCCTACCTGAGGCTGA
GGAAGGAGAAGGGGATGCACTGTTGGGGAGGCAGCTGTAACTCAAAGCCT
TAGCCTCTGTTCCCACGAAGGCAGGGCCATCAGGCACCAAAGGGATTCTG
CCAGCATAGTGCTCCTGGACCAGTGATACACCCGGCACCCTGTCCTGGAC
ACGCTGTTGGCCTGGATCTGAGCCCTGGTGGAGGTCAAAGCCACCTTTGG
TTCTGCCATTGCTGCTGTGTGGAAGTTCACTCCTGCCTTTTCCTTTCCCT
...
GCCCCCACTGCCTAGGGACCAACAGGGGCAGGAGGCAGTCACTGACCCCG
AGACGTTTGCATCCTGCACAGCTAGAGATCCTTTATTAAAAGCACACTGT
TGGTTTCTGCTCAGTTCTTTATTGATTGGTGTGCCGTTTTCTCTGGAAGC

取中间6行来比对，-x 以逗号分隔

bowtie2 -x bowtie2_index/GRCh38_noalt_as/GRCh38_noalt_as -c TAGCCTCTGTTCCCACGAAGGCAGGGCCATCAGGCACCAAAGGGATTCTG,CCAGCATAGTGCTCCTGGACCAGTGATACACCCGGCACCCTGTCCTGGAC,ACGCTGTTGGCCTGGATCTGAGCCCTGGTGGAGGTCAAAGCCACCTTTGG,TTCTGCCATTGCTGCTGTGTGGAAGTTCACTCCTGCCTTTTCCTTTCCCT,AGAGCCTCCACCACCCCGAGATCACATTTCTCACTGCCTTTTGTCTGCCC,AGTTTCACCAGAAGTAGGCCTCTTCCTGACAGGCAGCTGCACCACTGCCT

6 reads; of these:
  6 (100.00%) were unpaired; of these:
    0 (0.00%) aligned 0 times
    0 (0.00%) aligned exactly 1 time
    6 (100.00%) aligned >1 times
0   0   chr1    13201   1   50M *   0   0   TAGCCTCTGTTCCCACGAAGGCAGGGCCATCAGGCACCAAAGGGATTCTG  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII  AS:i:0  XS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:50 YT:Z:UU
1   16  chr15   101977667   1   50M *   0   0   GTCCAGGACAGGGTGCCGGGTGTATCACTGGTCCAGGAGCACTATGCTGG  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII  AS:i:0  XS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:50 YT:Z:UU
2   0   chr1    13301   30  50M *   0   0   ACGCTGTTGGCCTGGATCTGAGCCCTGGTGGAGGTCAAAGCCACCTTTGG  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII  AS:i:0  XS:i:-6 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:50 YT:Z:UU
3   16  chr2    113600039   1   50M *   0   0   AGGGAAAGGAAAAGGCAGGAGTGAACTTCCACACAGCAGCAATGGCAGAA  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII  AS:i:0  XS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:50 YT:Z:UU
4   16  chrX    156026608   1   50M *   0   0   GGGCAGACAAAAGGCAGTGAGAAATGTGATCTCGGGGTGGTGGAGGCTCT  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII  AS:i:0  XS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:50 YT:Z:UU
5   0   chr1    183970  1   50M *   0   0   AGTTTCACCAGAAGTAGGCCTCTTCCTGACAGGCAGCTGCACCACTGCCT  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII  AS:i:0  XS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:50 YT:Z:UU

6 个序列都比对上了但是有多个匹配，其中第 0, 2 个序列比对到了基因上, 其比对位置chr1 13201， 13301在 'DDX11L1' chr1 [11869 14409] 范围内，其他都没有。再次验证genome 里重复序列很多。

下面从 gencode.v38.annotation.gtf 文件里的随便挑一个 gene 'OR4F5' 的一个外显子，在chr1 [65520, 65573] 上。

chr1    HAVANA  exon    65520   65573   .   +   .   gene_id 'ENSG00000186092.7'; transcript_id 'ENST00000641515.2'; gene_type 'protein_coding'; gene_name 'OR4F5'; transcript_type 'protein_coding'; transcript_name 'OR4F5-201'; exon_number 2; exon_id 'ENSE00003813641.1'; level 2; protein_id 'ENSP00000493376.2'; hgnc_id 'HGNC:14825'; tag 'RNA_Seq_supported_partial'; tag 'basic'; tag 'Ensembl_canonical'; tag 'MANE_Select'; tag 'appris_principal_1'; havana_gene 'OTTHUMG00000001094.4'; havana_transcript 'OTTHUMT00000003223.4';

将其序列从 hg38.fa 里取出来

# 计算其所在的大概行
echo $((65520/50 +1)) $((65573/50 + 1 +1))
# 第 1311-1313行
# 测试比对后应该在1312-1314 行
head -n 1314  hg38.fa | tail -n 3 
# 以下下输出的是 chr1 65501-65650 的 序列

ACATATTATTCTCTTACAGTTTTTATGCCTCATTCTGTGAAAATTGCTGT
AGTCTCTTCCAGTTATGAAGAAGGTAGGTGGAAACAAAGACAAAACACAT
ATATTAGAAGAATGAATGAAATTGTAGCATTTTATTGACAATGAGATGGT

将上面 150 bp 序列使用 bowtie2 比对

bowtie2 -x bowtie2_index/GRCh38_noalt_as/GRCh38_noalt_as -c ACATATTATTCTCTTACAGTTTTTATGCCTCATTCTGTGAAAATTGCTGTAGTCTCTTCCAGTTATGAAGAAGGTAGGTGGAAACAAAGACAAAACACATATATTAGAAGAATGAATGAAATTGTAGCATTTTATTGACAATGAGATGGT

结果为

1 reads; of these:
  1 (100.00%) were unpaired; of these:
    0 (0.00%) aligned 0 times
    0 (0.00%) aligned exactly 1 time
    1 (100.00%) aligned >1 times
0   0   chr1    65501   1   150M    *   0   0   ACATATTATTCTCTTACAGTTTTTATGCCTCATTCTGTGAAAATTGCTGTAGTCTCTTCCAGTTATGAAGAAGGTAGGTGGAAACAAAGACAAAACACATATATTAGAAGAATGAATGAAATTGTAGCATTTTATTGACAATGAGATGGT  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII  AS:i:0  XS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:150    YT:Z:UU

可以发现其比对到了多个位置，最好的是在 chr1 65501 符合预期。