上次内容我们介绍了一个人全基因组完整分析脚本，这次我们来介绍一个完整的细菌基因组分析。本次实验，有一个大小为4.5M的细菌基因组，采用illumina双末端测序，测序读长为pair-end 90bp，测序深度为100X。通过拼接得到该物种的全基因组序列，并进行基因预测，基因注释等分析。
分析流程图
一、数据质控
利用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/联系客服
