Megahit是一個(gè)簡(jiǎn)單易用，且內(nèi)存需求低，對(duì)micro variation recovery rate較高，assembly N50表現(xiàn)都比較優(yōu)異的組裝軟件，EBI-metagenomic 在線分析平臺(tái)目前就是使用metaSPAdes（默認(rèn)），Megahit和Minia三個(gè)作為候選組裝軟件
然而，在拿到contig之后，我們有時(shí)還想估計(jì)read mapping back的情況，以及coverage和未mapping read的情況。下面的內(nèi)容引自Megahit在Github上的講解
4. 計(jì)算contig coverage和提取未用于組裝的reads
Calculate contig coverage and extract unassembled reads
下面的分析主要使用 BBMap 和samtools。我們可以使用相似的代碼分析我們自己的數(shù)據(jù)
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接受來(lái)源于多個(gè)測(cè)序文庫(kù)的數(shù)據(jù). 比如兩個(gè)雙端測(cè)序文庫(kù)和一個(gè)單端測(cè)序文庫(kù)可以使用這個(gè)代碼分析：in=pe1_1.fq,pe2_1.fq,se.fq in2=pe1_2.fq,pe2_2.fq -
kfilter是連續(xù)匹配堿基數(shù)的最小值，通常設(shè)置為k_min + 1 -
maxindel限制最大插入indel的長(zhǎng)度 -
subfilterread和contig比對(duì)的最大錯(cuò)配數(shù)
4.3 用pileup.sh輸出每個(gè)contig的coverage cov.txt with:
bbmap/pileup.sh in=aln.sam.gz out=cov.txt
4.4 提取unmapped reads (單端數(shù)據(jù)用unmapped.se.fq，雙端數(shù)據(jù)用 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
