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

Posted 2023-03-13

参考技术A -3.仅作笔记使用

我所用到的所有脚本可以+Q 840842906索

简单介绍一下方向也可以，说不定以后能互相帮忙～

注明来意，请实名  学校+名字

-2.两组样品，目的是通过宏基因组找到两组间的差异以及差异的成因

-1.原理

略

0. Megahit/SPAdes组装

略

1.基因预测

1.1 软件安装 - metagenemark/prodigal

略

参考：https://www.jianshu.com/p/b8d0cd39bb67    ---  genemark

          http://wap.sciencenet.cn/blog-2379401-1085480.html      ---    https://github.com/hyattpd/Prodigal

1.2

genemark参考：https://www.sohu.com/a/207036304_99908466

for i in B2A,B3A......          (下面会省略，$i的意思即中括号内的内容)

do

prodigal -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

done

-i是输入文件即组装好的contig，这里不同的文章会过滤掉不同长度的contig，过滤掉不同长度的预测结果的CDS，可以自行查阅文献

我这里以两篇学位论文为例：  (当然我没有混组装，我是每个样本组装一次，预测一次, 然后把所有样本的预测结果cat到一起

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

(1)sed '/>/s/gene/B2A_gene/' $i/$i.gene.fa | sed '/>/s/|GeneMark.hmm//' | awk 'print $1' > $i/$i.gene.fa2

(2)awk -F "|" 'print $1' $i/$i.gene.fa2 > $i/$i.gene.fa3

(3)bioawk -c fastx 'print $name, length($seq)' $i/$i.gene.fa3 > $i/$i.length

(4)awk 'if ($2 > 99) print $1' $i/$i.length > $i/$i.keep

(5)samtools faidx $i/$i.gene.fa3 -r $i/$i.keep > $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.length

(1)和(2)是重命名，把每个基因命名为 B2A_gene_1、B2A_gene_2......

(3)是获取每个基因的长度    (4)是得到长度大于100的基因名    (5)是用samtools提取前面长度大于100的基因名 (6)是清楚中间文件

比较丑陋，见谅

在过滤prodigal的时候猛然发现bbmap有现成的工具：

(其实是在运行atlas的时候翻看atlas的命令行发现的，atlas是一个宏基因组分析的snakemake流程：https://github.com/metagenome-atlas/atlas)

java -ea -Xmx1g -cp ./opt/bbmap-37.99/current/ jgi.RenameReads in=./prodigal/$i/$i.gene.fa out=$i.rename.gene.fa ow=t prefix=$i.gene minscaf=100

prefix=$i.gene就是把基因重命名为 B2A.gene.1  B2A.gene.2 。。。。。。

minscaf=100即过滤掉长度小于100      in=即输入文件  out=即输出文件

1.4  把所有样品，重命名/长度过滤后的gene.fa  cat到一起

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，在集群上很慢要一个周

参数的解释https://www.sohu.com/a/207036304_99908466, 个人感觉参数可以跟着上面照抄

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

3.1 ref : https://www.jianshu.com/p/a28c1729168e?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 > 2) print $1,$2;else print $1,“0”’ > $i.count.csv    #获取到每个gene比对到的reads数量,比对数小于等于2的直接标0, 接着进行长度标准化

3.3 ref : https://www.jianshu.com/p/a28c1729168e?ivk_sa=1024320u  这个链接里面相对丰度的计算是用的SOAP，但是我没有找到SOAP相关的任何信息，所以用的bwa mem

然后还需要获取非冗余基因集(即CDHIT的结果)每个gene的长度：

bioawk -c fastx ‘print $name, length($seq)‘ 去冗余的结果 | tr “\t” “,” > gene.length.csv      第一列是基因名/序列名，第二列是该序列长度

然后获取gene.count.csv中每个gene的长度

大体上就是用gene.count.csv的第一列和gene.length,csv的第一列匹配，匹配得上，就输出gene.count.csv的每一列(基因名，count)和gene.length,csv的第二列(基因长度)

如果匹配不上，就输出gene.count.csv的第一列, 该基因的count为0，长度

python jiyinfengdu.py $i.count.csv $i > $i.tmp

现在再根据上面那个标准化的公式来算就很简单了，很多方法都能实现

3.4：再将每个样本的数据合到一块即可

paste B2A.tmp B3A.tmp B4A.tmp > end.tsv

后续就可以进行可视化、功能注释等

待续

用到的脚本

2022.4更新

关于基因相对丰度的计算，又发现了一些其他的方法 (宏基因组的定量还没有一个很好很权威的方法

本文所用方法已经是12年的了(虽然才看到几篇文章还在用)

1.MOCAT2  (16年发在Bioinformatics)

好就好在方便很多，组装预测丰度一条龙全包了

2.Wen, C., Zheng, Z., Shao, T., Liu, L., Xie, Z., Le Chatelier, E., ... & Li, H. (2017). Quantitative metagenomics reveals unique gut microbiome biomarkers in ankylosing spondylitis. Genome biology

这篇发的好一点，也给出了一个方法