一个细菌基因组完整分析脚本

分析流程图

一个细菌基因组完整分析脚本

利用fastqc软件对原始测序reads进行指控，生成网页帮统计报告，根据报告内容对数据机型过滤。

mkdir result   fastqc -f fastq -o result _I249_FCC2BDTACXX_L4_SZAIPI030696-112_1.fq.gz _I249_FCC2BDTACXX_L4_SZAIPI030696-112_2.fq.gz

利用fastp软件对原始测序数据进行过滤，去取低质量，adapter，N碱基以及Duplication等，得到cleandata。

#利用fastp进行数据过滤  
fastp -i 130801_I249_FCC2BDTACXX_L4_SZAIPI030696-112_1.fq.gz -I 130801_I249_FCC2BDTACXX_L4_SZAIPI030696-112_2.fq.gz -o clean.1.fq.gz  -O clean.2.fq.gz  -z 4 -q 20 -u 30 -n 10 -h clean.html

利用fastp过滤得到的cleandata直接进行序列拼接。序列拼接过程中可以选择多种拼接软件进行拼接，如果是二代测序数据，可以一次选择多个测序文库进行拼接。也可以结合二代测序与三代测序得到更好的拼接结果。

#SOAPdenovo   mkdir kmer35   SOAPdenovo-63mer all -s lib.list -K 35 -o kmer35/kmer35 -D 1 -d 1 -u 2 -p 2>kmer35.log

#SPAdes   python /ifs1/Software/biosoft/SPAdes-3.12.0-Linux/bin/spades.py -o illumina_result -1 clean.1.fq.gz  -2 clean.2.fq.gz -t 2

补洞可以直接使用二代测序数据，也可以采用更长的sanger数据或者三代数据等。通过补洞可以得到更加完整的基因组，便于后续分析。

#补洞     GapCloser -a kmer35.scafSeq -b lib.list -o kmer35.fill.fa -l 100 -p 25 -t 1

对序列拼接结果进行统计，首先看拼接出来的基因组大小是否与真实基因组大小接近，然后在看拼接条数，最大长度，最小长度，N50等指标。如果拼接完整性和准确性较好，就可以用于后续分析了，否则还需要重新进行拼接。

#拼接结果统计   fa_stat kmer35.fill.fa Statistics:        Scaffold    Contig Total Number (#):    212     463 Total length (bp):          Gap(N)(bp):        26367       0 Average Length (bp):    20735.99        9437.71 N50 Length (bp):          37293 N90 Length (bp):    27653       9001 Maximum Length (bp):           Minimum Length (bp):    100     42 GC content :        65.58%      65.58%

#将最终拼接好的样品修改名字为样品名 mv kmer35.fill.fa sample1.fna

原核生物基因预测比较容易，可以利用自身基因组序列作为训练集，直接将拼接得到的基因组序列sample1.fna输入给prodigal软件，即可得到基因的核酸序列以及基因的氨基酸序列等，非常方便。

#原核生物基因预测 prodigal -a sample1.pep -d sample1.cds -f gff -g 11  -o sample1.gff -p single -s sample1.stat -i sample1.fna >prodigal.log

利用prodigal软件得到的基因的氨基酸序列sanple.pep文件与各种蛋白数据库进行比对即可得到基因功能注释信息，如果数据库还包括聚类以及富集信息，同理，可以注释得到相应的内容。这里利用eggnog-mapper软件一次可以完成多个功能数据库的注释。

mapper.py -i test/polb.fa --output polb_bact -d bact --data_dir /ifs1/Software/biosoft/eggnog-mapper-1.0.3/data/

将拼接得到的基因组sample1.fna直接输入给软件rnammer，几秒钟之后即可得到细菌的核糖体结果，包括5S，16S，23S等。其中16S可以用于构建系统发育树。

rnammer -S bac -m tsu,lsu,ssu  -gff sample1.gff -f sample1.frn sample1.fna

将拼接得到的基因组sample1.fna直接输入给软件rRNAscan，几秒钟之后即可得到细菌的装运RNA结果，包括转运RNA的二级结构等。

tRNAscan-SE  -B  -o tRNAScan.out -f tRNAScan.out.structure -m stat.list sample1.fna

将拼接得到的基因组sample1.fna直接输入给软件RepeatMasker，几分钟之后即可得到细菌的重复序列结果，

mkdir repeat
RepeatMasker -pa 2 -q -species tuber -html -gff -dir repeat sample1.fna

将拼接得到的基因组sample1.fna直接输入给软件trf，几秒之后即可得到细菌的串联重复序列结果，也就是vntr区域，vntr可以用于多样品的MLST分析，也可以用于CNV比较。

trf sample1.fna 2 7 7 80 10 50 500 -f -d -m

以上完成了一个原核生物常见的基因，核糖体，重复序列分析之后，可以进行统计，这些内容占据一个细菌基因组90%以上的内容，更多内容可以通过在线分析完成。如下所示。

#在线分析
基因预测：http://prodigal.ornl.gov/
http://ccb.jhu.edu/software/glimmer/index.shtml
RNAmer  http://www.cbs.dtu.dk/services/RNAmmer/
tRNAscan http://lowelab.ucsc.edu/tRNAscan-SE/
Go功能富集DAVID：http://david.abcc.ncifcrf.gov/
功能注释RAST：http://rast.nmpdr.org/
trf预测：http://tandem.bu.edu/trf/trf407b.linux.download.html
CPG岛预测：
http://www.ebi.ac.uk/Tools/emboss/cpgplot/index.html
http://www.ebi.ac.uk/Tools/seqstats/emboss_cpgplot/
TADB：http://bioinfo-mml.sjtu.edu.cn/TADB/
必须基因注释：http://tubic.tju.edu.cn/deg/
操纵子预测：
http://csbl.bmb.uga.edu/DOOR/index.php
http://operondb.cbcb.umd.edu/cgi-bin/operondb/operons.cgi
基因岛预测：http://www.islander.com/
IslandViewer：http://www.pathogenomics.sfu.ca/islandviewer/query.php
MobilomeFINDER：http://db-mml.sjtu.edu.cn/MobilomeFINDER/
CRISPR预测：http://crispr.u-psud.fr/Server/Advanced.CRISPRfinder.php
预测启动子：http://www-bimas.cit.nih.gov/molbio/proscan/
预测致病岛：https://www.gem.re.kr/paidb/about_paidb.php
预测分析转录终止信号
http://linux1.softberry.com/berry.phtml?topic=polyah&group=programs&subgroup=promoter
预测噬菌体：http://phaster.ca
预测信号肽：http://www.cbs.dtu.dk/services/SignalP/
Galaxy：https://usegalaxy.org/

系列教程：微生物组入门 Biostar 微生物组宏基因组

专业技能：学术图表高分文章生信宝典不可或缺的人

一文读懂：宏基因组寄生虫益处进化树

必备技能：提问搜索 Endnote

文献阅读热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读分析流程统计绘图

16S功能预测 PICRUSt FAPROTAX Bugbase Tax4Fun

在线工具：16S预测培养基生信绘图

科研经验：云笔记云协作公众号

编程模板: Shell R Perl

生物科普: 肠道细菌人体上的生命生命大跃进细胞暗战人体奥秘

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。

一个细菌基因组完整分析脚本