下面是基于组装的宏基因组数据分析流程

目录

基本流程介绍

megahit组装

什么是N50?

基于拼接结果的基因预测

cdhit去冗余

功能注释

宏基因组的分箱操作

分箱的目的：

分箱的原理：

---

单独对每个样本进行基因集组装，得到genome1,2,3，然后对组装的基因集进行编码基因预测，得到gene set1,2,3，然后把预测的基因集进行聚类合并得到unigenes，也就是已经去冗余的基因集，利用unigenes作为参考来做后续的基因功能的注释和丰度的估计。

简而言之，基因组的组装就是先把短的reads拼接成连续的序列contig，再把contig拼接成scaffold

"reads"、"contigs" 和 "scaffolds" 是描述数据组装过程中不同阶段的标准术语：

1. Reads：

1. Contigs：

1. Scaffolds：

简单来说，可以将这个过程比喻为拼图游戏：

• Reads 就像是拼图的小块，每个小块只包含有限的信息。
• Contigs 相当于将多个小块拼接起来形成的较大片段，但这些片段可能还不包括整个图像。
• Scaffolds 则是在 contigs 的基础上进一步确定它们在整幅拼图中的位置，形成的更大、更完整的结构，尽管其中可能还存在一些未知的空隙。

在基因组组装的过程中，科学家们的目标是尽可能地将所有的 reads 组装成更长的 contigs 和 scaffolds，以获得对基因组的完整理解。

组装使用的软件是megahit，输入是双末端的测序数据，megahit可以对指定的正向和反向测序数据进行组装， megahit 相对于其他的组装软件metasapdes来说运行速度更快，占用内存较小，组装准确率不会差特别多，主要是我们服务器内存太小了，只能用这个。一般服务器内存超过1T才推荐用metasapdes。

注：megahit输出的结果是contigs，metasapdes输出的结果是scaffold

megahit在运行的时候会分别取不同的kmer进行组装，然后整合组装结果。

这个过程是比较慢的

date="231228" input_path=~/Metagenome/clean_data/no_homo_genes/"${date}" output_path=~/Metagenome/megahit_out tmp_path=~/Metagenome/tmp #S1_1.clean.fq.gz cd "${input_path}" for file1 in *_1.clean.fq.gz;do file2="${file1/_1/_2}" sample_id=$(echo "$file1" | cut -d'_' -f1 | sed 's/S//') megahit \ -1 "$file1" \ -2 "$file2" \ --min-contig-len 500 \ --tmp-dir "${tmp_path}" \ --memory 0.8 \ --num-cpu-threads 20 \ --out-dir "${output_path}"/"${date}-${sample_id}" \ --out-prefix "S${sample_id}" done cd ~/Metagenome

解释一下megahit的参数

-1和-2用于指定双末端测序数据文件

--min-contig-len用于指定最小的拼接长度，因为太小的拼接长度是没有意义的，通常会设置成500或者1k

--tmp-dir用于指定临时文件输出到哪里去，因为运行megahit的时候会产生一些临时文件，这个参数如果不指定会输出到服务器的公共tmp里面去，那个文件很小，如果满了之后megahit就会报错退出，所以通常需要使用 --tmp-dir参数来指定临时文件输出到某个自己具有可写权限的路径下去，我们服务器就输出到磁盘3

--memory 0.8的意思是如果megahit运行的时候占用内存到了整个服务器内存的百分之八十，就会自动把这个任务kill掉，所以通常可以设置的高一点

--num-cpu-threads指定线程数

--out-dir和--out-prefix分别指定输出路径和输出文件名的前缀，比如我上面代码中的S${sample_id}，当sample_id是4的那次循环，输出结果就是这样一个文件夹

我们用的是S4.contigs.fa，长这样子，每一段时一个拼接好的序列

组装完成之后可以使用quest来对组装拼接的质量进行评估，运行的主程序放在了容器里面，可以直接调用

singularity exec ~/Metagenome/software/MetaGenome.sif quast.py ./*.fa

调用quast.py程序即可完成汇总，传入的参数非常简单，就是要汇总的文件，我们这里就是以.fa结尾的文件（需要提前把所有拼接结果也就是那些.fa结尾的文件放到一个文件夹下）,运行结束之后会生成一个名为quast_resault的文件夹（目录），里面有这些文件

如果想把我们的total length和N50写到文章里面就可以去这个html

或者pdf里面去找。

在基因组学中，N50是一个衡量DNA测序和组装质量的指标。它代表了所有contigs（或scaffolds）按长度从长到短排列后，达到总长度中位数的最短长度。具体来说，N50的计算方式是将所有的contigs按长度降序排列，然后累加它们的总长度，直到累加的总长度达到所有contigs总长度的50%，此时对应的那个contig的长度就是N50值。如果这个值较大，通常意味着组装的质量和连续性较好，因为较长的contigs表示较少的片段间隙。

基因预测使用的软件是prokka，prokka是一个整合了许多其他软件的流程，可以做古菌，细菌，线粒体，病毒的基因预测

代码如下，这一步也比较慢

prokka \ --outdir prokka_out \ --prefix tumor \ --addgenes \ --addmrna \ --locustag A1 \ --kingdom Bacteria \ --metagenome \ --cpus 8 \ ./A1.contigs.fa

--outdir指定输出到哪个文件夹去，因为要输出的文件比较多。

--prefix指定输出文件的前缀，到时候prokka会输出很多文件到参数outdir指定的文件夹去

--addgenes 和--addmrna 表示把基因信息和RNA信息加上去

--locustag A1表示在预测的每个基因名前面加上A1，这是为了后续的区分。如果不写，会导致每个样本预测的基因很多都是一样的名字。

--kingdom Bacteria 表示要预测细菌的基因

--metagenome 表示用宏基因组模式运行

--cpus 8 指定线程数是8

最后加上输入文件，也就是上一步megahit 组装的.fa文件

输出结果有很多格式的文件，不过重要的只有两个

分别是：

.faa : 预测的氨基酸序列文件，重要，后面要做功能注释

.ffn ：预测出的转录本序列

重复上述过程，完成其他样品的基因预测。

使用prokka完成多个样本的基因预测之后会生成很多.faa和.ffn文件，我们把这些文件合并起来，进行如下处理

#合并所有预测的转录本文件 cat ~/Metagenome/prokka_out/*.ffn > all.trans.fa #合并所有预测的蛋白文件 cat ~/Metagenome/prokka_out/*.faa > all.pep.fa #统计转录本和氨基酸序列的条数，转录本的应该是要比氨基酸的多 transcript_count=$(grep -c ">" all.trans.fa) echo "基于拼接片段预测到的转录本数量为：$transcript_count" pep_count=$(grep -c ">" all.pep.fa) echo "基于拼接片段预测到的氨基酸序列数量为：$pep_count_count" #把所有蛋白的ID抓出来存到all.cds.id sed -n "s/^>\(\S\+\).*$/\1/p" all.pep.fa > all.cds.id #根据前面的id把编码蛋白的RNA序列提取出来存到all.cds.fa里面 seqtk subseq all.trans.fa all.cds.id > all.cds.fa

seqtk subseq all.trans.fa all.cds.id > all.cds.fa表示从 all.trans.fa文件里面提取all.cds.id对应的序列

运行上面这段代码就拿到了一个名为all.cds.fa的文件，它包含了所有编码蛋白的RNA序列，但是由于我们是把所有的样本预测结果直接合并起来了，而不同样本中肯定会预测到很多相同的基因，所以我们要去冗余，这是为了得到一个作为参考的预测编码序列来对每个样本进行定量，也就是看看每个样本里面这些基因的表达量有多少。去冗余使用的工具是cd-hit-est，代码如下

#使用cd-hit-est进行编码基因序列的去冗余 cd-hit-est \ -i all.cds.fa \ -o all.cds.cdhit \ -c 0.95 \ -aS 0.9 \ -M 10000 \ -T 4

其中-i指定输入文件名，-o指定输出文件名，-c 0.95意思是两条序列重复部分达到了百分之95就会把这两条序列聚类成同一个cluster，-aS 0.9意思是聚类要求相同序列至少要占到比对的较短的序列的百分之90才会把两条序列聚成一个cluster，也就是说只有两条序列重复区域达到百分之95且重复区域占到较短的那条序列的百分之90以上才会把这两条序列聚成一个cluster

-M参数不用管，这是和内存相关的一个参数，-T用来指定线程数

cd-hit-est运行结束之后会生成两个文件，一个是我们用-o参数指定的文件，另一个是-o参数指定的文件名加上一个cluster的文件

其中-o参数指定的文件也就是all.cds.cdhit长这样，这就是去冗余之后的基因序列了，可以看到有A1样本作为代表的，下面还有C1样本作为代表的，不过没截图，这是因为我们在前面运行prokka的时候使用了样本名作为区分，如果一个基因在样本只在A1里面预测到了，那没问题，但是如果一个基因在很多样本里面都被预测到了，我们前面如果不使用样本名来进行区分，就不知道这个基因是在哪个样本里面被预测到的了

all.cds.cdhit.cluster这个文件长这样，去冗余之后会生成很多cluster，如果一个cluster里面有多条序列，就选取最长的那条作为代表序列，代表序列后面有一个星号

之后我们使用命令

cp all.cds.cdhit unigene_cds.fasta

对去冗余之后的文件名进行一个修改，改成unigene_cds.fasta，这样可以更清楚的知道这个文件是干啥的。

然后我觉得前面生成的all.cds.cdhit文件（也就是刚才重命名生成的unigene_cds.fasta文件）里面预测到的基因有各种样本名开头的，不太规范，所以我想要修改一下这个文件里面的基因名，让他看起来更规范一些，这里我每个基因名换成对应的cluster名，这样就避免了基因名里面带有不同的样本名。

首先把去冗余之后的基因名抓取出来存到unigene.id里面，因为后面要做蛋白注释，所以顺手就用这个unigene.id把对应的蛋白序列提取出来了，存到了unigene_pep.fasta，顺手的事，和改基因名没关系，只是后面要用，且正好改基因名的时候生成的unigene.id文件正好可以用来提取蛋白，就顺带着做了

awk '$1 ~/^>/{print $1}' all.cds.cdhit | sed 's/^>//' > unigene.id seqtk subseq all.pep.fa unigene.id > unigene_pep.fasta

使用下面的代码来处理all.cds.cdhit.clstr文件

awk '{if($1~/^>/){printf $1 $2} else if($NF ~ /*$/){print "\t"$3}}' all.cds.cdhit.clstr |sed 's/>//g; s/...$//'|awk '{print $2"\t"$1}' > map_id.txt

最终可以得到一张这样的表也就是文件map_id.txt，第一列是旧id，第二列是新id，如图所示

有了这张表之后运行下面这两句代码即可完整两个fasta文件中基因序列id的替换，map_fasta_ids是一个用perl写的脚本

~/Metagenome/my_script/map_fasta_ids map_id.txt unigene_cds.fasta ~/Metagenome/my_script/map_fasta_ids map_id.txt unigene_pep.fasta

以上就是宏基因组的数据拼接和基因预测的过程，我们得到了两个文件，分别是预测到的且去冗余的编码基因序列文件 unigene_cds.fasta和蛋白序列文件unigene_pep.fasta ，后续分析都是基于这两个文件。

接下来进行预测到的基因的功能注释与丰度估计。

功能注释用的是eggnog数据库，推荐使用eggnog-mapper网页版进行在线注释，要比本地快很多，可以登陆 http://eggnog-mapper.embl.de 来进行在线注释，提交数据之后需要在邮件里面确认一下才会开始进行注释，不然就只是提交了数据。分析结束之后会收到一封邮件，里面有一个以emapper.annotation结尾的文件，在线注释每一次上传reads的上限是10w条，但是通常我们预测得到的去冗余序列文件要高于10w条，所以需要先做一个序列的切分。上传的时候我们直接使用去冗余的蛋白序列，这样就可以跳过翻译等步骤，让结果出的快一点。切分使用的工具是seqkit里面的split工具，代码如下

seqkit split \ --by-size 100000 \ --out-dir split \ unigene_pep.fasta

这样就可以把切分好的序列输出到split这个目录里面去

然后 登录 http://eggnog-mapper.embl.de/ 上传刚才切分好的蛋白序列进行注释，之后要对返回的结果进行合并，使用下面的代码

awk '{if( NR==FNR || $1 !~ /^#/){print $0} }' ./unigene*.emapper.annotations > unigene.annotations

把eggnog网站返回的结果合并成一个名为unigene.annotations的文件，他包含了预测基因的注释结果。

然后对结果进行统计绘图并构建orgdb，orgdb是使用R包clusterProfiler进行富集分析的时候要传的东西。运行过程中需要调用一个已经写好了的脚本emapperx.R，可以使用下面的代码直接运行，也可以手动去运行这个R脚本，这个脚本运行所需要的文件就是unigene.annotations和unigene_pep.fasta

## 对结果进行统计绘图并构建orgdb Rscript ~/Metagenome/my_script/emapperx.R unigene.annotations unigene_pep.fasta

emapperx.R脚本的输出结果有这些，对应着不同数据库的注释结果，比如GO，cog，KEGG（pathway）

值得注意的是结果中除了不同数据库的注释结果和图之外还有一个anno_stat.txt文件，如图，他展示了有多少基因注释到了对应的数据库，其中第一行Eggnog显示一共有百分之63的基因注释到了，其中注释到GO数据库的只有百分之5.76，这显然是太少了。kegg数据库有22.59，这是正常的，因为kegg数据库比较严格，研究明确的通路数量也没那么多。

也就是说对于原核生物来讲，不推荐使用eggnog输出的GO注释结果。推荐使用uniplot数据库来做GO的注释，需要先构建一个数据库，使用diamond模式，diamond 是一个用于蛋白质和翻译后DNA序列比对的软件，跟blast做的事情差不多，但是速度要比blast快很多，内存占用也小。

运行下面的代码构建比对数据库

## 构建diamond database 运行一次即可 uniref90_path=/data3/sunjs/DataBase/uniref90.fasta.gz diamond makedb --in "${uniref90_path}" --db uniref90

构建diamond数据库需要提前从uniplot下载一个uniref90的数据库，我下好了已经，- --db uniref90表示构建的数据库的前缀是uniref90，运行结束之后会生成一个名为uniref90.dmnd的文件，这个就是diamond比对所需的数据库了。

然后开始进行diamond比对

## 进行diamond比对 diamond blastp \ --db ../../Database/uniprot/uniref90 \ --query unigene_pep.fasta \ --out unigene_pep.uniref90.m6 \ --threads 20 \ --outfmt 6 \ --max-target-seqs 5 \ --evalue 1e-5

其中--db指定刚才构建好的数据库，但是注意只写uniref90就行，不用写dmnd这个后缀。

--query指定要比对的文件也就是我们前面的去冗余的蛋白序列。

--out指定输出文件名。

--outfmt 6 表示输出表格模式。

--max-target-seqs 5 输出比对的最好的5条，这参数不用动。

--evalue 1e-5 不用动。

有了diamond blastp 的输出文件也就是unigene_pep.uniref90.m6，就可以根据从uniref下载的另一个数据库idmapping_selected.tab.gz来把GO注释信息提取出来，用的是下面这段代码

## 基于idmapping 提取GO注释信息 perl ~/Metagenome/my_script/uniref90_idmapping.pl \ unigene_pep.uniref90.m6 \ /data3/sunjs/DataBase/idmapping_selected.tab.gz \ > unigene_pep.uniref90.GOanno

其中uniref90_idmapping.pl是一个提前准备好的脚本，该脚本需要的文件是diamond blastp 的输出文件unigene_pep.uniref90.m6，以及从uniref下载的数据库idmapping_selected.tab.gz，最终生成一个名为unigene_pep.uniref90.GOanno的文件，如图

第一列是我的序列ID，第二列是该序列对应的GO注释信息。

（课程里面没有将怎么把这个结果可视化，可能就就是挑选注释到的最多的一些GO编号用柱状图展示一下注释到了多少条基因吧，跟前面emapperx.R生成的pdf文件类似，可以参考一下）

所谓功能注释也不过是在数数，数的是有多少基因在某条通路富集，除了功能注释，prokka输出的两个结果（也就是去冗余的基因序列文件和蛋白序列文件）还可以进行丰度估计，注意这里所说的丰度估计可不是前面kraken2+bracken的那个丰度估计了，前面说的是物种的丰度估计，现在是基因的丰度估计，也就是基因的表达量，虽然不一定是同一物种在表达这个基因。

基因丰度估计使用的工具是 Salmon，使用salmon定量之前需要构建一个索引，代码如下

## Salmon mkdir salmon_out # 构建salmon index⽤于⽐对，使用的是前面去冗余过的cds文件 salmon index \ -t ~/Metagenome/prokka_out/unigene_cds.fasta \ -i my_unigene_index

其中-t是构建索引的参考文件也就是前面去冗余的基因序列文件，-i指定生成的索引文件名。

接下来运行salmon中的quant来进行丰度估计

salmon quant \ --validateMappings \ --meta \ -p 30 \ -i my_unigene_index \ -l IU \ -1 ../../dataFQ/A1_1.fq.gz \ -2 ../../dataFQ/A1_2.fq.gz \ -o quants/A1.quant

其中--validateMappings是一个增加准确性的参数

--meta表示是宏基因组的数据

-p 30 表示使用30个线程

-i my_unigene_index是上一步生成的索引

-l IU 不用改，常规双末端测序这个参数都是IU，表示双末端测序，非链特异性

-1，-2指定双末端测序的两个文件

-o指定输出目录的名称，因为salmon的输出结果是一个目录名为quants，目录里面有子目录，-o指定的就是这个目录中子目录的名字，quants是自动生成的。

运行结束之后-o指定的子目录中有一个文件叫quant.sf，这个文件长这样子，有用的主要是TPM和NumReads两列

如果有多个样本，就每个样本都跑一遍然后把.sf表格整合起来，就可以拿到一张基因的丰度表。其中TPM这列代表的就是基因丰度了他是 是一种用于评估基因表达水平的标准化指标，它能够反应基因的丰度。

合并每个样本的salmon结果的代码如下

cd path_to_qunants #这地方要改，改成实际的quants路径 quants=` ls quants/ | awk '{print "quants/"$1 }' | tr '\n' ' ' ` names=` ls quants/ | sed 's/.quant$//'| tr '\n' ' ' ` #合并所有样本的TPM列 salmon quantmerge \ --quants $quants \ --names $names \ --column tpm \ -o unigenens.tpm #合并所有样本的NumReads列 salmon quantmerge \ --quants $quants \ --names $names \ --column numreads \ -o unigenens.count #unigenens.tpm_fig是生成文件的前缀 Rscript draw_abundance.R unigenens.tpm ../../data/sample.txt unigenens.tpm_fig

解释一下quants=` ls quants/ | awk '{print "quants/"$1 }' | tr '\n' ' ' `

首先是列出quants文件夹里面的所有文件，然后把这些文件通过管道符传递给awk，awk '{print "quants/"$1 }'可以给这些文件每一个前面加上quants/

tr '\n' ' ' 的作用是：将输入中的所有换行符转换为空格

也直接自己手动一个个写，比如

#合并所有样本的TPM列 salmon quantmerge \ --quants quants/A1.quant quants/A2.quant quants/A3.quant \ --names A1 A2 A3 \ --column tpm \ -o unigenens.tpm

这跟前面的效果是一样的，因为

quants=` ls quants/ | awk '{print "quants/"$1 }' | tr '\n' ' ' `

names=` ls quants/ | sed 's/.quant$//'| tr '\n' ' ' `

这两句看着很复杂的代码做的也就是生成quants/A1.quant quants/A2.quant quants/A3.quant这样的东西和生成A1 A2 A3这样的东西。

通过整合TPM这列，我们可以得到一张基因的丰度表，比如下面这张图就整合了两个样本的salmon结果中的TPM列

为了做差异分析，我们还需要整合最后一列也就是NumReads这列，因为后续做差异分析的DEseq2需要不同基因序列的对应的总reads数量。

合并完成之后生成了两个文件，也就是前面salmon quantmerge中

-o参数指定的unigenens.tpm和unigenens.count，其中tpm文件是用来做一些基因表达量的相关性分析或者绘制热图用的，绘图用的脚本在/data1/sunjs/Metagenome/my_script这个路径下，绘图用到的文件有unigenens.tpm，前面kraken2物种鉴定用的那个分组文件sample.txt，就是这个文件，第一列是分组，第二列是样本

绘图可以直接打开R脚本来一步步运行。主要生成的图有

箱线图：

x轴是不同的样本，y轴是log10(TPM+1)，反映了每个样本中基因表达量的整体水平，我画的时候要改成每个组别的，不然样本太多了，不过这个图不是很重要。+1是因为有的基因在某些样本里面TMP是0，没法取对数，为了防止报错，加了1。箱线图中百分之九十以外的数据是用点来表示的，可以看到有一个个的点。在这里每一个点代表一个基因，也就是矩阵中某一行。

热图：

这个热图是基于TMP矩阵绘制的，反映了不同样本之间的相关性。可以看到每个样本跟自己组别里面其他样本的相关性是很高的。截图没截完整，下面x轴也是A1，A2.....C3

上面那个热图不带聚类树，下面这个带上了，和上面那个含义一样。可以看到AB两组是比较相近的。

PCA图

也是基于TPM矩阵来做的，主要反映了不同样本的基因表达趋势，基因表达趋势差不多的会被聚在一起，可以看到不同组别被区分的很明显。

TMP的密度图

反映的是不同表达量的基因的数量

count文件是用来做差异分析用的，差异分析的脚本如下

#处理一下unigenens.count的表头格式，让其符合后续输入格式 sed '1s/^Name\t//' unigenens.count > unigenens.count.matrix ## 差异基因分析，这个脚本在/data1/sunjs/Metagenome/my_script里面 #和找差异物种的时候用的是同一个程序 run_DE_analysis.pl \ --matrix unigenens.count.matrix \ --method DESeq2 \ --samples_file sample.txt \ --contrasts contrasts.txt \ --output DE_out ## A vs C组差异基因进行筛选 ## 标准为|log2FC| >1 && FDR <0.05 ## 链接eggnog注释时创建的R_Library/ ln -s ../05.uniprot/R_Library/ ln -s ../04.eggnog/unigene.annotations Rscript \ ../script/enrich_analysis.R \ ./DE_out/unigenens.count.matrix.A_vs_C.DESeq2.DE_results \ ./unigene.annotations \ ./DE_out/unigenens.count.matrix.A_vs_C.DESeq2.DE_results.enrich

其中有一个参数--contrasts，这个参数可以不指定，不指定的话默认就是不同分组两两进行比较，比如我有三个分组分别是A,B,C，但是我只想要比较AC两组之间的差异，我就可以创建一个名为contrasts.txt的文件，他的内容如下

这样就只会比较AC两组的差异了，标准为|log2FC| >1 && FDR <0.05。不过我们的课题都是来源于病人和健康人的两个组别的，所以这个参数通常不需要添加。

差异基因分析用的脚本和我们找差异物种的时候用的脚本是一样的，最后会输出差异分析的结果和一个巨丑的火山图，如图

这个火山图我直接不要，我们选择用输出的差异分析的表格导入到R里面自己画火山图。DEseq2输出结果会包含所有的基因，我们可以在R里面筛选出显著差异的基因（P_adj<0.05,logFC>1）然后去做富集分析。做富集分析需要orgdb，我们的orgdb就是前面eggnog注释的时候创建的，有一个文件叫R_Library，把他先链接过来，这里面包含的是GO的注释，同时把unigene.annotations也链接过来，这个是我们用eggnog网站在线注释之后返回结果合并起来的文件，这里面包含了KEGG的注释。有了这两个文件以及前面DEseq2结果里面的*DESeq2.DE_results结果这三个文件，就可以进行富集分析了，使用的脚本是enrich_analysis.R，这个脚本也放在了/data1/sunjs/Metagenome/my_script里面。

分箱的对象可以是reads，contigs，或者gene，但是一般情况下都是做contigs的分箱，目的是获得一个单菌的基因组，然后就可以从基因组的角度去做这个单菌的研究。

1.不同物种的核酸组成是不一样的，比如GC的含量等等，我们可以把核酸组成相似的那些分类单元聚在一起认为他们同一物种

2.同一物种的测序深度理论上来讲是相同的，假如某个物种有10个，那我这10个个体的基因组打碎测序，应该是每条基因都能测到10次，但是由于测序的随机性，实际上不是都能测到10次，可能是8次9次之类的，一条reads被测到的次数就是测序深度，如果一些基因的测序深度相近且核酸组成也相似，我们就怀疑他们是同一个物种。