利用宏基因组数据组装基因组-评估篇

前言：

最近组装了一种病原体的基因组，基因组大小为610kb，结果发现在300,000-400,000之间发现很多的Gap区域，需要找一下原因。因为是用二代数据测的，我先推测的原因是基因组这个区域有可能GC含量比较高，那下载一下它的基因组，看一下，找到了bedtools工具，发现这个软件功能十分强大，bedtools总共有二三十个工具/命令来处理基因组数据。比如：根据bed中的位置信息提取目标基因及其上下游序列；统计基因组不同区间的GC含量；提取gff文件的所有基因位置,并转换成bed格式，计算覆盖度（coverage）（coverageBed，genomeCoverageBed）等等。

一．Bedtools评估基因组二代数据覆盖度

分析思路：把这个基因组文件按照100,000bp大小，按照碱基位置分割成6个文件，然后分别计算不同区间的GC含量。

1.软件下载安装：

https://github.com/arq5x/bedtools2

wget https://github.com/arq5x/bedtools2/releases/download/v2.30.0/bedtools-2.30.0.tar.gz

tar -zxf bedtools-2.30.0.tar.gzcd bedtools2/make

然后用help看一下，安装成功。

 2.软件使用：

准备基因组文件，如果是多个序列文件，比如我的基因组文件，先计算各个contig/scaffold的长度：

使用python里面的pyfaidx模块的faidx命令，代码如下：

conda activate py36 (我自己的一个py36小环境)

pip install pyfaidx

faidxminia_k81.contigs.fa -i chromsizes > size.genome

   结果如下：

划分窗口：

/public/home/rp1016swf/rp1016swf/software/bedtools2/bin/bedtools makewindows -g sizes.genome -w 50000 > windows. Bed

-g sizes.genome是要划分的基因组，格式为两列：染色体、染色体长度-w 50000为窗口大小为5wwindows.bed为输出文件，格式为三列：染色体、区间开始位点、区间结束位点。

 统计窗口内的GC含量：

/bedtools2/bin/bedtools nuc -fi ViralProj237323_genomic.fna -bed windows.bed | cut -f 1-3,5 > 5w_gc.bed

统计窗口内的平均覆盖深度

bedtools coverage -a windows.bed -b SRR081241.sorted.bam > RR081241.depth.txt

bedtools coverage对划分好的每个滑动窗口进行reads数（depth)的统计。

-a windows为上一步划分好的区间

-SRR081241.sorted.bam为测序数据mapping到参考基因组的比对文件

HG00096.depth.txt为统计结果的输出文件，格式为7列：染色体、区间起始位点、区间结束位点、该区间内的reads数、该区间内的碱基数、区间大小、该区间的平均覆盖度

生成的txt文件共有7列，分别为序列编号、起始位置、结束位置、reads数、碱基数、区间大小、平均覆盖深度

二、Bowtie2+Samtools评估基因组二代数据比对率

Bowtie下载安装

wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.5.1/bowtie2-2.5.1-linux-x86_64.zip/

1.先是构建索引：

#bowtie2-build -f ./minia_k81.contigs.fa -p ./minia_k81.contigs -p索引文件前缀名

2.bowtie2比对及samtools转为bam文件,并根据比对情况进行提取

bwa比对生成的为sam（sequence Alignment mapping)文件，将SAM转换为二进制对应的BAM格式。二进制格式对于计算机程序来说更容易使用。要将SAM转换为BAM，我们使用samtools view命令。

bowtie2 -x ./minia_k81.contigs -p 20 -1 20220652_mapped_P1.fq -2 20220652_mapped_P2.fq -S  ./minia_k81.contigs.sam

3.准备序列比对后生成的 bam 文件或者 .sam 文件

samtools view -bS minia_k81.contigs.sam > minia_k81.contigs.bam

4.统计序列比对情况

samtools flagstat scaffolds.bam > flagstatstat.txt

本文使用 文章同步助手 同步