宏基因组数据二+三代混合组装并计算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-build

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

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  1A_contig_sorted 31. 47476 opera_contig_2 contigs  1A_contig_sorted 86.  opera_contig_3 contigs  1A_contig_sorted 35. 39262 opera_contig_4 contigs  1A_contig_sorted 87. 60715 opera_contig_5 contigs  1A_contig_sorted 36. 23955 opera_contig_6 contigs  1A_contig_sorted 29. 18003 opera_contig_7 contigs  1A_contig_sorted 56.074216 33284 opera_contig_8 contigs  1A_contig_sorted 77. 44505 opera_contig_9 contigs  1A_contig_sorted 42. 23506 

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

待解决的问题

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