宏基因组-(1)基因预测和基因相对丰度计算

我用的脚本都可以从+Q 840842906下载。

简单介绍一下方向就可以了，说不定以后还能互相帮助呢~

请说明你的目的。你的真名？学校+名称

-2.两组样本，目的是通过宏基因组学找到两组之间的差异以及差异产生的原因。

-1.原则

省略

0.超级打击/黑桃组合

省略

1.基因预测

1.1软件安装-meta gene mark/浪子

省略

参考:/blog-2379401-1085480 . html？-/hyattpd/浪子

1.2

Genemark参考资料:/a/207036304_99908466

为我在{B2A，B3A......} ?(下面省略，$i表示括号内的内容。)

做

浪子-a $ I/$ I . pro . fa-d $ I/$ I . gene . fa-f gff-g 11-I $ I/$ I . fasta-p meta-s $ I/$ I . stat _ potential . gene-m-q

完成的

-i是输入文件，即组装的重叠群。这里不同的文章会筛选出不同长度的contig和不同长度预测结果的CDS，大家可以自行查阅文献。

我举两篇论文为例:(当然，我没有混编。我每个样本组装一次，预测一次，然后把所有样本的预测结果放在一起。

1.3然后我们需要给预测的基因重新命名，过滤长度。

(1)sed '/& gt；/s/gene/B2A _ gene/' $ I/$ I . gene . fa | sed '/& gt；/s/| gene mark . hmm//' | awk ' { print $ 1 } ' & gt；$i/$i.gene.fa2

(2)awk-F " | " ' { print $ 1 } ' $ I/$ I . gene . fa2 & gt；$i/$i.gene.fa3

(3)bioawk -c fastx '{print $name，length($ seq)} ' $ I/$ I . gene . fa3 & gt；$i/$i .长度

(4)awk“{ if(＄2 & gt；99)print $ 1 } ' $ I/$ I . length & gt；$i/$i.keep

(5)Sam tools faidx $ I/$ I . gene . fa3-r $ I/$ I . keep & gt；$i/$i.filter.gene.fa

(6)RM $ I/$ I . gene . fa2；RM $ I/$ I . gene . fa3；RM $ I/$ I . gene . fa3 . fai；RM $ I/$ I . keep；rm $i/$i .长度

(1)和(2)正在重命名，每个基因命名为B2A _基因_1和B2A _基因_2。......

(3)获取每个基因的长度(4)获取长度大于100的基因名(5)用samtools提取之前长度大于100的基因名(6)作为清晰的中间文件。

丑，见谅。

在过滤浪子的时候，突然发现bbmap有一个现成的工具:

(其实我是在运行atlas的时候看atlas的命令行发现的。atlas是宏基因组分析的snakemake过程:/metagenome-atlas/atlas)。

java -ea -Xmx1g -cp。/opt/bbmap-37.99/current/ jgi。RenameReads in=。/浪子/$ I/$ I . gene . fa out = $ I . rename . gene . fa ow = t prefix = $ I . gene mins caf = 100

Prefix=$i.gene表示将基因重命名为B2A.gene.1？B2A .基因.。。。。。

Minscaf=100，即过滤长度小于100？In=即输入文件？Out=输出文件。

1.4?重命名/长度所有样本，过滤基因. fa？一起猫

2.CD-Hit冗余删除

CD-hit-I input . fa-o out . fa-c 0.95-d 0-aL 0.9-uL 0.05-aS 0.9-T 80-M 15000

我用的是-T 80，去学校的超算，一个下午就做完了；；以前用的-T 14，集群需要一周时间。

参数解读/a/207036304_99908466，个人感觉参数可以从上面复制。

3.相对丰度和长度标准化的计算

3.1 ref:/p/a28c 1729168 e？ivk_sa=1024320u？这个环节相对丰度的计算用的是SOAP，但是我没有找到SOAP相关的资料，所以bwa mem包括Salmon和bowtie2+bbmap也可以做同样的事情。

去除冗余后的Bwa mem -p结果样本1.1.fq.gz样本1.2.fq.gz >样本1.sam

获得样品1上每个基因的比对。

3.2

grep-v ^@ $ I . AlN . Sam | cut-f 3 | sort | uniq-c | awk ' begin { fs = " "；OFS= "，" }{print $2，$1}' | awk 'BEGIN{ FS= "，"；OFS= "，" } { if(＄2 & gt；2)打印$1，$ 2；否则打印$1，" 0 " } ' & gt$i.count.csv #获取每个基因匹配的阅读数，直接将小于等于2的阅读数标记为0，然后标准化长度。

3.3 ref:/p/a28c 1729168 e？ivk_sa=1024320u？这个环节相对丰度的计算是基于SOAP的，但是我没有找到SOAP相关的资料，所以用了bwa mem。

然后你需要得到非冗余基因集中每个基因的长度(也就是CDHIT的结果):

Bioawk-c fastx' {print $ name，length ($ seq)} '冗余结果| tr' \ t ' '，' >；基因.长度. csv？第一列是基因名称/序列名称，第二列是序列长度。

然后得到gene.count.csv中每个基因的长度

一般来说，gene.count.csv的第一列与gene的第一列匹配。长度和csv。如果匹配，则每列gene.count.csv(基因名称，计数)和第二列gene。长度和csv(基因长度)输出。

如果不匹配，则输出gene.count.csv的第一列，基因的计数为0，长度为0。

python jiyinfengdu . py $ I . count . CSV $ I & gt；$i.tmp

现在按照上面标准化的公式计算非常简单，很多方法都可以实现。

3.4:然后结合每个样本的数据。

粘贴b2a . tmp b3a . tmp b4a . tmp & gt；end.tsv

可视化、功能注释等。可以稍后执行。

待续

使用的脚本

2022.4更新

关于基因相对丰度的计算，已经发现了一些其他的方法(对于宏基因组的量化还没有好的权威的方法)

本文使用的方法是12岁(虽然只有少数文章还在使用)。

1.MOCAT2？(发表于《生物信息学》16)

好处是方便多了，组装和丰度预测都在一个包里。

2.文，郑，赵，邵铁涛，刘，李，谢，，...& amp李，H. (2017)。定量宏基因组学揭示强直性脊柱炎独特的肠道微生物生物标志物。基因组生物学

这篇文章比较好，也给出了方法。