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

本文转自“基因学苑”，已获授权

分析流程图

一、数据质控

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

mkdir result  
fastqc -f fastq -o result 130801_I249_FCC2BDTACXX_L4_SZAIPI030696-112_1.fq.gz 130801_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):    4396029     4369662
Gap(N)(bp):        26367       0
Average Length (bp):    20735.99        9437.71
N50 Length (bp):    115412      37293
N90 Length (bp):    27653       9001
Maximum Length (bp):    230197      198997
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/

八、核糖体RNA预测

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

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

九、转运RNA预测

将拼接得到的基因组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/

猜你喜欢

10000+：菌群分析 宝宝与猫狗 梅毒狂想曲 提DNA发Nature Cell专刊 肠道指挥大脑

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

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

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

必备技能：提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

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

16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun

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

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

编程模板: Shell  R Perl

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

写在后面

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

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”