基因组测序深度的估算方法之一—bwa/samtools
Published on 2020-05-19 10:48
Author: zhixy
2774 views
方法介绍
基因组de novo测序，没有参考序列，如要计算测序深度，可以将reads回贴到拼接好的contigs/scaffolds序列上，近而统计每个核酸位点上被reads覆盖的次数，即测序深度。
回贴mapping短序列reads到基因组上，成熟的软件有bwa、bowtie2等，下面以bwa为例介绍一种测序深度估算方法。
方法流程
1. 用短序列拼接结果，建立bwa索引。
(base) [user@server ~]# bwa index contigs.fa
2. 短序列mapping
(base) [user@server ~]# bwa mem -t 60 -M contigs.fa strain_trim1.fastq strain_trim2.fastq > bwa.sam
-t 60，计算核心/线程数
3. sam转bam
(base) [user@server ~]# samtools view -@ 60 -bS bwa.sam > bwa.bam
-@ 60，计算核心/线程数
4. bam结果排序
(base) [user@server ~]# samtools sort -@ 60 bwa.bam > bwa.sort.bam
5. 计算每个位点的覆盖度
(base) [user@server ~]# samtools depth bwa.sort.bam > bwa.depth
6. 统计平均覆盖度
(base) [user@server ~]# cat bwa.depth | awk '{ sum ＋= ＄3; } END { print "coverage = " sum/NR }'
对第三列（$3）求和，并除以行数（NR），行数即基因组序列长度。