BSA分析(四)——序列比对及比对信息统计

拿到质控后的数据就可以开始进行数据的比对了。(这里复习一下)需要确定好参考基因组,参考基因组选择标准:质量高、完整、和某一亲本相近(是亲本之一的基因组当然更好了)

比对软件

由于转录组存在可变剪切,所以和DNA重测序数据还是有区别的。简单来说就是DNA重测序数据包含基因组所有的碱基信息(内含子外显子都有),但是RNA数据只有表达数据(外显子,同时存在可变剪切)。所以比对时要注意这一点。DNA数据一般使用BWA,RNA除了BWA还可以使用Hisat2等。

BWA

下载安装

wget http://superb-sea2.dl.sourceforge.net/project/bio-bwa/bwa-0.7.12.tar.bz2
tar jxf bwa-0.7.12.tar.bz2
cd bwa-0.7.12 && make

#编译完成后将软件写入环境变量中
echo 'PATH=$PATH:~/softwarepath/bwa-0.7.12/setup' >> ~/.bashrc && source ~/.bashrc

#测试是否安装好
bwa 
Program: bwa (alignment via Burrows-Wheeler transformation)
Version: 0.7.12
Contact: Heng Li 

Usage:   bwa  [options]

Command: index         index sequences in the FASTA format
         mem           BWA-MEM algorithm
         fastmap       identify super-maximal exact matches
         pemerge       merge overlapping paired ends (EXPERIMENTAL)
         aln           gapped/ungapped alignment
         samse         generate alignment (single ended)
         sampe         generate alignment (paired ended)
         bwasw         BWA-SW for long queries

         fa2pac        convert FASTA to PAC format
         pac2bwt       generate BWT from PAC
         pac2bwtgen    alternative algorithm for generating BWT
         bwtupdate     update .bwt to the new format
         bwt2sa        generate SA from BWT and Occ

Note: To use BWA, you need to first index the genome with `bwa index'.
      There are three alignment algorithms in BWA: `mem', `bwasw', and
      `aln/samse/sampe'. If you are not sure which to use, try `bwa mem'
      first. Please `man ./bwa.1' for the manual.

使用

bwa使用手册

bwa包含三种算法:

BWA-backtrack: 比对Illumina 序列,read length: ≤100bp。

BWA-SW: 比对 long read,也支持剪切比对。read length:70bp-1Mbp。

BWA-MEM: 比对 long read,也支持剪切比对。read length:70bp-1Mbp。(该算法更新,准确性和速度都更优)

一般使用步骤:

#ref.fa即为参考基因组,genome为建立索引库的前缀,不加-p genome,默认建立索引以ref.fa为前缀
bwa index ref.fa -p genome

#因为BSA分析中的重测序数据大多是150 PE,这里就直接介绍mem算法的用法,别的算法不再详细介绍
bwa mem [-aCHMpP] [-t nThreads] [-k minSeedLen] [-w bandWidth] [-d zDropoff] /
        [-r seedSplitRatio] [-c maxOcc] [-A matchScore] [-B mmPenalty] [-O gapOpenPen] /
        [-E gapExtPen] [-L clipPen] [-U unpairPen] [-R RGline] [-v verboseLevel] /
        db.prefix reads.fq [mates.fq]

常用参数注释:
-t INT  线程数,默认是1。
-M      将 shorter split hits 标记为次优,以兼容 Picard’s markDuplicates 软件。
-p      若无此参数:输入文件只有1个,则进行单端比对;若输入文件有2个,则作为paired reads进行比对。若加入此参数:则仅以第1个文件作为输入(输入的文件若有2个,则忽略之),该文件必须是read1.fq和read2.fa进行reads交叉的数据。
-R STR  完整的read group的头部,可以用 '/t' 作为分隔符,sam文件会识别。需要指定RG表头是因为同一样品可包括多个不同lane、不同文库的测序结果,此外不同样品比对结果合并时,也需要通过RG进行标记区分。
-T INT  当比对的分值比 INT 小时,不输出该比对结果,这个参数只影响输出的结果,不影响比对的过程。
-a      将所有的比对结果都输出,包括 single-end 和 unpaired paired-end的 reads,但是这些比对的结果会被标记为次优。

SE(单):
bwa aln [options] ref.fa read.fq > aln_sa.sai
bwa samse [options] ref.fa aln_sa.sai read.fq > aln-se.sam
PE(双):
bwa aln [options] ref.fa read1.fq > aln_sa1.sai
bwa aln [options] ref.fa read2.fq > aln_sa2.sai
bwa sampe [options] ref.fa aln_sa1.sai aln_sa2.sai read1.fq read2.fq > aln-pe.sam

#实例(双端),sam文件即为我们的比对结果文件。
##不指定样品名
bwa mem ref.fa read1.fq read2.fq > out.sam
##指定样品名,sample即为样品名
bwa mem -M -R "@RG/tID:sample/tLB:sample/tSM:sample"  ref.fa read1.fastq read2.fastq > out.sam 2> out.sam.log

当我们想要提高比对的速率的时候,可以直接使用shell的基本命令split,或者写python脚本对测序数据进行分割再并行比对(每四行为一组),最后再将比对文件合并起来。

注意:1、因为涉及到后面的变异检测,所以当基因组存在染色体长度大于500M时,一定要对染色体拆分成小于500M的序列再进行比对,在变异检测完毕后再进行合并(拆分基因组的脚本我会放在流程脚本里,这里不再详细记录);2、一定不能是基因组拆分进行比对,因为数据是无序的,且涉及到多重比对的问题!!!

比对结果格式解读

@HD VN:1.5  SO:coordinate
@SQ SN:Chr01    LN:91897250
@RG ID:sample   LB:sample   SM:sample 
@PG     ID:bwa  CL:~/softwarepath/bwa-0.7.12/setup/bwa mem -M -R "@RG/tID:sample/tLB:sample/tSM:sample"  ref.fa read1.fastq read2.fastq > out.sam 2> out.sam.log      PN:bwa  VN:0.7.17-r1188
E150015117L1C019R0062820889     433     Chr01       186     0       72M78H  ChrD02_S6       89876638        0       AAAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCTA        ?A+<80-?C5-@A5BC<7>A<@CC-2A9CCB67@2BCC/;8:@CCC=66@CCA<>=BCCCAB68B@C)@AB<8C4=?2A=A76=6/867/C=4>A--=(2B=A1;@@/9.#@'4C>).@?92/<.@@CB9C.C?CC#CBBA6BC:C6ACCBCC>CC=@CCCA7A>CACAC6CCACC?1?CCC=B9CCCCC.2-CCC<;C4CCCC;5?CCCCC-CCCAC&(CCC;A*C5CAC7C'9BBCCC?7CCC>8B  NM:i:8  MD:Z:5C11C6C13^C0C81C0C13C14    MC:Z:3M1D147M   AS:i:108        XS:i:106        RG:Z:CY619

#@开头的列都是表头信息行,不能少
HD:(排序类型)
头部区第一行:VN是格式版本;SO表示比对排序的类型,有unknown(default),unsorted,queryname和coordinate几种。
SQ:染色体号/t染色体长度
RG:样本信息
PG: 比对命令行,可以直接由此查到比对方式,参考基因组,测序数据,样本名等信息

#比对内容行每列内容注释:
第一列:QNAME。比对的序列名称,类似于E150015117L1C019R0062820889。
第二列:FLAG Bwise FLAG,比对类型:paring,strand,mate strand等。(这个值需要进行计算)
第三列:RENAME 比对上的参考序列名,类似于Chr01
第四列:POS 1-Based的比对起始物理位置
第五列:MAPQ 比对质量
第六列:CIGAR Extended CIGAR string(操作符:MIDNPSH=X)比对结果信息,包含数字和字符。
第七列:MRNM 相匹配的另外一条序列,比对上的参考序列名,“=”:比对到一个位置;“*”:无匹配;染色体号:比对多个位置。
第八列:MPOS 1-Based leftmost Mate Position 与该reads对应的mate pair reads的比对位置
第九列:ISIZE 插入片段长度
第十列:SEQ 和参考序列在同一个链上比对的序列(若比对结果在负义链上,则序列是其反向重复序列,反向互补序列)
第十一列:QUAL 比对序列的质量(ASCII-33=Phred base quality)reads碱基质量值
可选列:以TAG:TYPE:VALUE的形式提供额外的信息。

比对文件格式处理

samtools

下载安装

# 建议直接conda安装
conda create -n samtools
conda activate samtools 
conda install -c bioconda samtools -y

#测试是否安装成功
samtools 

Program: samtools (Tools for alignments in the SAM format)
Version: 1.3.1 (using htslib 1.3.1)

Usage:   samtools  [options]

Commands:
  -- Indexing
     dict           create a sequence dictionary file
     faidx          index/extract FASTA
     index          index alignment

  -- Editing
     calmd          recalculate MD/NM tags and '=' bases
     fixmate        fix mate information
     reheader       replace BAM header
     rmdup          remove PCR duplicates
     targetcut      cut fosmid regions (for fosmid pool only)
     addreplacerg   adds or replaces RG tags

  -- File operations
     collate        shuffle and group alignments by name
     cat            concatenate BAMs
     merge          merge sorted alignments
     mpileup        multi-way pileup
     sort           sort alignment file
     split          splits a file by read group
     quickcheck     quickly check if SAM/BAM/CRAM file appears intact
     fastq          converts a BAM to a FASTQ
     fasta          converts a BAM to a FASTA

  -- Statistics
     bedcov         read depth per BED region
     depth          compute the depth
     flagstat       simple stats
     idxstats       BAM index stats
     phase          phase heterozygotes
     stats          generate stats (former bamcheck)

  -- Viewing
     flags          explain BAM flags
     tview          text alignment viewer
     view           SAM<->BAM<->CRAM conversion
     depad          convert padded BAM to unpadded BAM

使用

这里主要介绍如何将sam文件处理成二进制bam文件,建立索引,排序,去重,以及统计比对信息

#sam文件查看
less -S out.sam
#bam文件查看
samtools view -hS out.sort.rmdup.bam |less -S 
##以上二者区别就是是否为二进制文件

#sam转bam
samtools view -hb out.sam >out.bam

#拆分并行比对后会涉及到bam合并的问题,merge也可以合并
samtools cat 1.bam 2.bam ... -o out.bam

#bam排序,先建索引(生成*.bam.bai文件,当染色体长度过长时则会出现索引无法建立的情况)再排序
samtools index out.bam
samtools sort out.bam -o out.sort.bam

#去重,s:SE,S:PE,默认双端
samtools rmdup -S out.sort.bam out.sort.rmdup.bam 

#统计比对信息
samtools flagstat out.sort.rmdup.bam >  flagstat.txt

#flagstat.txt文件内容
7524 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
7523 + 0 mapped (99.99% : N/A)
7524 + 0 paired in sequencing
3760 + 0 read1
3764 + 0 read2
7510 + 0 properly paired (99.81% : N/A)
7522 + 0 with itself and mate mapped
1 + 0 singletons (0.01% : N/A)
12 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

第一行:QC pass的reads的数量为7514,未通过QC的reads数量为0,意味着一共有7524*2条reads;
第四行:重复reads的数量,QC pass和failed,这里都是0条
*第五行:比对到参考基因组上的reads数量,这里是共7523条比对上了参考基因组,比对率为99.99%;
第六行:paired reads数据数量,这里是7524对;
第七行:read1比对上的数量;
第八行:read2比对上的数量;
*第九行:双端都正确匹配到参考序列的reads数量,这里是7510,比对率为99.81%;
第十行:双端都匹配到了参考序列上的reads数量(但是并不一定比对到同一条染色体上);
*第十一行:双端中只有一条与参考序列相匹配的数量;
*第十二行:双端比对到不同染色体的数量;
第十三行:双端比对到不同染色体的且比对质量值大于5的数量。
## *为主要数据关注行

#覆盖率统计
/work1/Software/samtools-1.4/setup/bin/samtools depth -a out.sort.rmdup.bam > out_depth.txt 2>out_depth.statlog
#out_depth.txt文件内容
#染色体号   物理位置    覆盖碱基数
Chr01      1       0
Chr01      2       10
Chr01      3       0
Chr01      4       0

## a为有碱基覆盖的位点总长度
a=`awk '$3!=0' out_depth.txt|wc -l `
## b是基因组总长度
samtools faidx ref.fa
b=`awk '{sum+=$2}END{print sum} ref.fai`
## a/b 即为覆盖率了

往期回顾

BSA分析(一)——原理及发展史

BSA分析(二)——分析准备工作

BSA分析(三)——测序数据的质控

参考资料

https://blog.csdn.net/genome_denovo/article/details/78712972

版权声明:
作者:admin
链接:https://www.techfm.club/p/81251.html
来源:TechFM
文章版权归作者所有,未经允许请勿转载。

THE END
分享
二维码
< <上一篇
下一篇>>