宏基因组单个样本数据处理流程笔记

前言

衷心感谢这个开源互助的网络时代，能够让我在36岁高龄重新回到程序员的行列。虽然人生会因为各种原因起起落落落落落落，但还好家学渊源不至于无事可做。我在落落落落落落期间接触生物信息基础，竟然发现：这不就是当年奥赛做错的动态规划大题么？真是山不转水转……在此特别感谢陈同和刘永鑫两位老师的生物信息交流平台，以及GitHub上各个开源工具的作者和用户们的无私分享。
本篇博文将总结分享我在进行宏基因组数据处理中的一些流程和经验，并根据工作进度和涉猎范围进行完善补充，尤其欢迎各位朋友能够在留言区进行交流。本文不单独介绍各个软件的安装调试过程，而是集中在具体使用。硬件使用的是自费购买搭建的Intel NUC10i7FNH平台，系统是Ubuntu 20.04LTS版，环境控制是Anaconda。有Dell的运算服务器当然更好，但有也不等于就一定能用，没有也不等于不能开工。代码只要想写，就一定能写，这跟让不让你做手术是不太一样的。

数据预处理

跟若干家公司打过交道，要么拒绝提供数据，要么仅仅提供处理过的CleanData。为了能够确保流程和数据的可靠性，特别是对于以后要撰写论文过程中的严谨性，我要求提供下机时生成的RawData，以及Adaptor和Barcode等配套数据。
这一部分内容主要参考羊城迷鹿的文章宏基因组测序流程，涉及数据质量控制、去除接头和去除宿主几个步骤。

质量控制

使用的是FastQC进行质量控制，使用也很方便：

# 这里会直接在当前目录下生成html报告文件，--threads 12参数表示使用12个线程。
$ fastqc --threads 12 seqs.fastq.gz

其实对任何生成的fastq文件，无论是RawData、CleanData还是中间处理过程中生成的文件，都最好能够随时进行质控，做到心中有数。

去除接头序列

目前我还没有真实的执行过去除接头序列的步骤，因为前面几家公司都没有提供这些原始数据供验证。我目前备用的工具是trimmomatic，等我实际操作之后再进行补充。

去除宿主序列

我使用bwa去除宿主序列。之前用过别的一些软件，但感觉bwa去除得干净一些。比对用的数据库是人类基因组的hg38数据库，貌似这个应该是最新的版本了。
注意在使用之前需要先对下载的数据库文件hg38.fa建立索引。

$ bwa index -a bwtsw ~/hg38/hg38.fa

建立索引之后就可以开始去除宿主序列，基本思路就是先进行数据比对，并将比对上的序列进行反转去除，剩下未能比对成功的，就是去除宿主之后的目标序列：

# 使用12个线程将fastq文件在hg38数据库中进行比较，并转换成sam文件。
$ bwa mem -t 12 -M ~/hg38/hg38 seqs.fastq.gz > seqs.sam
# 将sam文件转换成bam，方便后续进行处理。
$ samtools view -bS seqs.sam > seqs.bam
# 筛选出没有在hg38当中被