宏基因组数据二+三代混合组装并计算Read对Contig的深度

二+三代混合组装：OPERA-MS

对于我这一批样品，分别使用了Illumina NovaSeq PE250、插入文库约400bp的双端测序，另外对于每一组平行样进行了Oxford Nanopore的单分子测序。

metaSPAdes组装短读长

由于混合组装软件OPERA-MS的运行中首先要做的步骤就是利用Megahit或者metaSPAdes进行短读长的组装，我正好直接先自己组装好，该软件可以接受短读长的组装contig文件作为输入跳过该软件内置的短读长拼接步骤。metaSPAdes这一步我用的是默认参数。

OPERA-MS混合组装

上一步得到的metaSPAdes组装结果为contig_spades.fasta。
OPERA-MS混合组装command：

perl ../OPERA-MS/OPERA-MS.pl \
	--short-read1 Illumina_R1.fq \
	--short-read2 Illumina_R2.fq \
	--contig-file contig_spades.fasta \
	--long-read nanopore_seq.fastq \
	--out-dir sample_name_dir \
	--no-ref-clustering \
	--num-processors 64 \
	--no-polishing

即使提供了预组装好的contig文件，参数中仍需提供质控后用于预组装的短读序列文件。另外，如果不写--no-polishing参数的话，组装结束后会花费巨额的时间使用Pilon软件对序列进行polish，不是很需要的可以使用该参数跳过该步骤，节约时间。

相比于metaSPAdes的耗时，OPERA-MS的速度并不算慢。

计算Contig丰度以及Read对Contig的深度

我们想将所有二代reads给map到拼接好的contig上，每一条contig中能map上的全部read数目作为该contig的丰度。在一篇公众号上看到说全部read的总Nt数除以该contig长度的值作为该contig的覆盖度。但是覆盖度这个概念我总觉得应该以百分比的形式展示，这个定义我觉得更偏向于是：该contig上每个Nt平均被read所测到的次数，这个定义我觉得更像是“测序深度”的定义？

同样，这一部分的计算我也是从上述的公众号上学来的，可能是bowtie2和samtools版本问题，我使用的参数形式可能和该公众号里的不太一样，简单记录一下。

我使用的版本：
Bowtie 2 version 2.4.1
samtools: 0.1.19-44428cd
CheckM v1.1.3

第一步：bowtie-build

bowtie2-build --thread 64 1A/contigs.fasta ~/hotspring/NGS/read_mapping/1A/1A_contig_build

会生成这些文件：

第二步：bowtie2

bowtie2 -p 64 \
	-x ~/hotspring/NGS/read_mapping/1A/1A_contig_build \
	-1 ~/hotspring/NGS/dedupe/1A_dedupe_R1.fq \
	-2 ~/hotspring/NGS/dedupe/1A_dedupe_R2.fq \
	-S ~/hotspring/NGS/read_mapping/1A/1A_contig.sam

生成的就是sam文件

第三步：samtools view将sam文件转化为bam文件 然后用samtools sort排序

#sam to bam
samtools view -bS -@ 64 ~/hotspring/NGS/read_mapping/1A/1A_contig.sam > ~/hotspring/NGS/read_mapping/1A/1A_contig.bam
#sort
samtools sort ~/hotspring/NGS/read_mapping/1A/1A_contig.bam \
~/hotspring/NGS/read_mapping/1A/1A_contig_sorted -@ 64 -m 1024M

这一步sam文件就变成了bam文件之后排序！

注意 我自己按照samtools的说明文档操作的时候使用了-n这个参数，运行下一步的时候就会报错！我去官方manual上查看了，这么写的：

Note that if the sorted output file is to be indexed with samtools index, the default coordinate sort must be used. Thus the -n and -t options are incompatible with samtools index.

所以这一步就不要使用-n啦！
这一步会生成sorted的bam文件

第四步：samtools index

samtools index ~/hotspring/NGS/read_mapping/1A/1A_contig_sorted

这一步生成了一个sorted.bam.bai文件

第五步：使用CheckM软件计算coverage

这一步之前先建一个folder，把sorted的bam文件，bam.bai文件和contig文件扔到一个文件夹里~

mkdir checkm_bins
cp 1A/contigs.fasta checkm_bins
cp 1A_contig_sorted.bam checkm_bins
cp 1A_contig_sorted.bam.bai checkm_bins 

然后就可以使用checkm了：

checkm coverage -x fasta -m 20 -t 64 checkm_bins contigs_coverage.out checkm_bins/1A_contig_sorted.bam

结果输出在contigs_coverage.out文件里。

查看结果文件

结果文件大概是这样的：

Sequence Id	Bin Id	Sequence length (bp)	Bam Id	Coverage	Mapped reads
opera_contig_1	contigs	359737	1A_contig_sorted	31.891943	47476
opera_contig_2	contigs	301539	1A_contig_sorted	86.816714	108701
opera_contig_3	contigs	266462	1A_contig_sorted	35.674667	39262
opera_contig_4	contigs	167722	1A_contig_sorted	87.231794	60715
opera_contig_5	contigs	156833	1A_contig_sorted	36.943685	23955
opera_contig_6	contigs	149703	1A_contig_sorted	29.145709	18003
opera_contig_7	contigs	142800	1A_contig_sorted	56.074216	33284
opera_contig_8	contigs	137842	1A_contig_sorted	77.849074	44505
opera_contig_9	contigs	132689	1A_contig_sorted	42.780042	23506

由于该软件实际上是查看分箱结果的，所以他的header是这样的。在这里，sequence id实际上就是我们所有的contig，后面显示了该contig的长度，map到该contig的read数目以及计算的“coverage”。我大概验算了一下，比如按照第一行，用这里的“coverage”乘以contig长度除以map到的read数：
31.891943*359737/47476 ≈ 241.65
我基本上推测他这个coverage真就这么算出来的~

至此，这一部分的计算就完成了。

待解决的问题

其实我对于sam，bam文件的内容和操作还不是很熟悉。之后进行基因预测等操作的时候也会用到（还会用到很多别的格式），这几天学习一下，可能的话再写一篇学习笔记。