More about sequence alignments

发表于

序列比对PPT

链接:http://pan.baidu.com/s/1i5Ia9wD 密码:jsio

一、延伸材料

Blast使用
https://www.ncbi.nlm.nih.gov/books/NBK1734/
Needleman-Wunsch algorithmdemo
http://experiments.mostafa.io/public/needleman-wunsch/
Dynamic Programming
http://www.avatar.se/molbioinfo2001/dynprog/dynamic.html
FastQ格式
https://en.wikipedia.org/wiki/FASTQ_format
SAM/BAM格式说明
http://samtools.github.io/hts-specs/SAMv1.pdf
SAM flags解释
https://broadinstitute.github.io/picard/explain-flags.html

二、常用比对软件及简单描述

  • 多序列比对
    http://www.ebi.ac.uk/Tools/msa/
    ClustalW，最广、经典
    MUSCLE，快
    T-coffee，相比ClustalW准确度高
    Prank，相比ClustalW准确性高
    MAFFT，速度准确度都好过ClustalW

速度：MUSCLE > MAFFT ~ Prank> ClustalW >T-Coffee
准确性：MAFFT > Prank > MUSCLE > T-Coffee > ClustalW

BWA为最常用的序列比对软件，可以处理几十bp到1Mb长的reads；
Bowtie2常用于RNA-seq、Chip-seq等，速度比BWA快；
SOAP2输出的比对格式与BWA、Bowtie2不同（不是SAM/BAM格式），使用时需要注意；
Samtools 是比对文件（SAM/BAM）处理相关最常用的工具软件。

三、比对相关命令（BWA）：

  1. 为参考序列Reference建库：

    1
    $ bwa index –a bwtsw ref.fa

  2. 使用BWA-MEM算法比对，保存为SAM格式

    1
    $ bwa mem ref.fa read.1.fq.gz read.2.fq.gz > align.sam

  3. 使用BWA-MEM算法比对，保存为BAM格式（压缩的SAM格式）

    1
    2
    $ bwa mem ref.fa lane1.read.1.fq.gz lane1.read.2.fq.gz | samtools view -bS -o align1.bam –
    $ bwa mem ref.fa lane2.read.1.fq.gz lane2.read.2.fq.gz | samtools view -bS -o align2.bam -

  4. 合并多个比对文件（比如同一个样品，分别在两条lane测序，生成两对reads）

    1
    $ samtools merge merge.bam align1.bam align2.bam

  5. 将比对结果进行排序

    1
    $ samtools sort align.bam align.sort

  6. 去除PCR重复

    1
    $ samtools rmdup output.sort.bam output.dedup.bam

  7. 使用GATK套件Indel附近区域进行重新比对（Realignment）

    1
    2
    $ java -jar GenomeAnalysisTK.jar -R ref.fa -T RealignerTargetCreator -o output.realn.intervals -I output.dedup.bam
    $ java -jar GenomeAnalysisTK.jar -R ref.fa -T IndelRealigner -targetIntervals output.realn.intervals -o output.realn.bam -I output.dedup.bam

四、处理比对文件（SAM/BAM格式）

  1. 打开SAM文件

    1
    $ less align.sam

  2. 打开BAM文件

    1
    2
    3
    $ samtools view align.bam |less     # 查看比对信息（无Header信息）
    $ samtools view -h align.bam |less     # 查看Header和比对信息
    $ samtools view -H align.bam |less     # 查看Header信息

  3. 查看比对深度

    1
    2
    $ samtools depth align.bam | less    
    $ samtools depth -r chr1:1000-2000 | less      # 只查看染色体chr1上1000-2000位置

  4. samtools主要功能描述
    image

试一试：

  • 使用下列序列进行搜索
    MSTAVLENPGLGRKLSDFGQETSYIEDNCNQNGAISLIFSLKEEVGALAKVLRLFEENDVNLTHIESRPSRLKKDEYEFFTHLDKRSLPALTNIIKILRHDIGATVHELSRDKKKDTVPWFPRTIQELDRFANQILSYGAELDADHPGFKDPVYRARRKQFADAYNYRHGQPIPRVEYMEEEKKTWGTVFKTLKSLYKTHACYEYNHIFPLLEKYCGFHEDNIPQLEDVSQFLQTCTGFRLRPVALLSSRDFLGGLAFRVFHCTQYIRHGSKPMYTPEPDICHELLGHVPLFSDRSFAQFSQEIGLASLGAPDEYIEKLATIYWFTVEFGLCKQGDSIKAYGAGLLSSFGELQYCLSEKPKLLPLELEKTAIQNYTVTEFQPLYYVAESFNDAKEKVRNFAATIPRPFSVRYDPYTQRIEVLDNTQQLKILADSINSEIGILCSALQKIK

从结果中获得以下信息：
该序列在人类中的基因描述；
第一个匹配结果的bit score得分、E值、Identities、Gap数目；
上述指标第二个匹配结果如何；
下载前两个匹配结果，选择任一比对软件将与原序列与下载的软件一起进行多序列比对，它们之间的差异在什么地方，与blast结果比对的结果有什么异同。

  • 使用NCBI数据库获得不少于10个物种BRCA2基因氨基酸序列，并应用Clustal软件（http://www.ebi.ac.uk/Tools/msa/clustalo/）进行序列多重比对得到系统进化树，分析其进化关系