Megahit 组装结果进一步分析

Megahit是一个简单易用，且内存需求低，对micro variation recovery rate较高，assembly N50表现都比较优异的组装软件，EBI-metagenomic 在线分析平台目前就是使用metaSPAdes（默认），Megahit和Minia三个作为候选组装软件

然而，在拿到contig之后，我们有时还想估计read mapping back的情况，以及coverage和未mapping read的情况。下面的内容引自Megahit在Github上的讲解

4. 计算contig coverage和提取未用于组装的reads

Calculate contig coverage and extract unassembled reads

下面的分析主要使用 BBMap 和samtools。我们可以使用相似的代码分析我们自己的数据

4.1 下载 & 安装 BBMap 和 samtools:

# BBmap
wget http://downloads.sourceforge.net/project/bbmap/BBMap_35.34.tar.gz
tar zvxf BBMap_35.34.tar.gz
# samtools
wget https://github.com/samtools/samtools/releases/download/1.2/samtools-1.2.tar.bz2
tar xvjf samtools-1.2.tar.bz2
cd samtools-1.2 && make -j && cd ..

4.2 Align reads with bbwrap.sh:

bbmap/bbwrap.sh ref=SRR341725.megahit_asm/final.contigs.fa in=SRR341725_1.fastq.gz in2=SRR341725_2.fastq.gz out=aln.sam.gz kfilter=22 subfilter=15 maxindel=80

Tips:

• bbwrap 接受来源于多个测序文库的数据. 比如两个双端测序文库和一个单端测序文库可以使用这个代码分析：in=pe1_1.fq,pe2_1.fq,se.fq in2=pe1_2.fq,pe2_2.fq
• kfilter 是连续匹配碱基数的最小值，通常设置为k_min + 1
• maxindel 限制最大插入indel的长度
• subfilter read和contig比对的最大错配数

4.3 用pileup.sh输出每个contig的coverage cov.txt with:

bbmap/pileup.sh in=aln.sam.gz out=cov.txt

4.4 提取unmapped reads (单端数据用unmapped.se.fq，双端数据用 unmapped.pe.fq):

samtools-1.2/samtools view -u -f4 aln.sam.gz | samtools-1.2/samtools bam2fq -s unmapped.se.fq - > unmapped.pe.fq