通过Prokka对基因组数据进行注释并提取16S序列
直接通过conda安装prokka
$ conda create -n prokka
$ conda activate prokka
$ conda install prokka
直接使用prokka进行基因注释
prokka ~/project/ONT/basecalling_gpu_sup/6_flye/assembly.fasta
最后的输出结果：
屏幕截图 2021-12-12 205903.png
可以查看得到的gff文件是否有16S这个关键信息，发现我们得到的注释基因组里面有多个拷贝。
contig_1 barrnap:0.9 rRNA 889083 890619 0 - . ID=FHNIADIC_00822;locus_tag=FHNIADIC_00822;product=16S ribosomal RNA
contig_1 barrnap:0.9 rRNA 1614887 1616423 0 - . ID=FHNIADIC_01508;locus_tag=FHNIADIC_01508;product=16S ribosomal RNA
contig_1 barrnap:0.9 rRNA 1655057 1656593 0 - . ID=FHNIADIC_01547;locus_tag=FHNIADIC_01547;product=16S ribosomal RNA
contig_1 barrnap:0.9 rRNA 1794724 1796260 0 - . ID=FHNIADIC_01673;locus_tag=FHNIADIC_01673;product=16S ribosomal RNA
contig_1 barrnap:0.9 rRNA 1913066 1914602 0 - . ID=FHNIADIC_01788;locus_tag=FHNIADIC_01788;product=16S ribosomal RNA
contig_1 barrnap:0.9 rRNA 2449422 2450958 0 + . ID=FHNIADIC_02302;locus_tag=FHNIADIC_02302;product=16S ribosomal RNA
contig_1 Prodigal:002006 CDS 2904854 2905270 . - 0 ID=FHNIADIC_02743;eC_number=3.1.-.-;Name=yqgF;db_xref=COG:COG0816;gene=yqgF;inference=ab initio prediction:Prodigal:002006,similar to AA sequence:UniProtKB:P0A8I1;locus_tag=FHNIADIC_02743;product=Putative pre-16S rRNA nuclease
contig_1 barrnap:0.9 rRNA 3252096 3253632 0 + . ID=FHNIADIC_03075;locus_tag=FHNIADIC_03075;product=16S ribosomal RNA
enter image description here
对GFF注释文件的解析：
#### **1. column1**
第一列是`seqid`, 代表序列ID, 通常是染色体的ID, 每条染色体拥有一个唯一的ID。
#### **2. column2**
第二列是`source`, 代表基因结构的来源，可以是[数据库](https://cloud.tencent.com/solution/database?from=10680)的名称，比如来自`genebank`数据库，也可以是软件的名称，比如用`GeneScan`软件预测得到，当然，也可以为空，用`.`点号填充。
#### **3. column3**
第三列是`type`, 代表区间对应的特征类型，比如`gene`, `exon`等。
#### **4. column4**
第四列是`start`, 代表区间的起始位置。
#### **5. column5**
第四列是`end`, 代表区间的终止位置。
#### **6. column6**
第六列是`score`, 软件提供了统计值，如果没有，就用`.`填充。
#### **7. column7**
第七列是`strand`, 代表正负链的信息, `+`表示正链，`-`表示负链，`?`表示不清楚正负链的信息，当正负链信息没有意义时，可以用`.`填充。
#### **8. column8**
第八列是`phase`，当描述的是CDS区间信息时，需要指定翻译时开始的位置，取值范围包括0,1,2。
#### **9. column9**
第九列是`attributes`, 表示属性，每种属性采用`key=value` 的形式，多个属性之间用`;`分号分隔。
对于正负链的探讨：
测序之前，要先将序列打断，进行建库。在建库的时候是不区分正负链的，所以在后续测序过程中产生序列也不包含正负链信息。
但是，测的read要往参考基因组上比对，以便知道read在基因组中的位置。参考基因组中每条染色体只有一条序列，这条序列表示的序列即为正链序列，常用‘+’表示，与正链互补的序列为负链，用‘-’表示。
对一条染色体而言，只有一条正链序列，互补的即为负链。
文献或者数据库中提到的基因所在strand，即为这里所提的正负链。
用samtools来提取序列信息
samtools faidx PROKKA_12122021.fna
得到一个PROKKA_12122021.fna.fai文件，于是接下来便通过samtools工具截取gff文件中提示的16S的序列，以编号为ID=FHNIADIC_00822的数据集举例。
samtools faidx PROKKA_12122021.fna contig_1:889083-890619
最后得到对应位置的序列，可以另行保存为fasta格式，另外可以看一下序列的大小，显示为1537bp
屏幕截图 2021-12-12 215913.png
同样，我们可以再看一条显示为16S的序列，以编号为ID=FHNIADIC_01508举例
samtools faidx PROKKA_12122021.fna contig_1:1614887-1616423
对应的序列信息，同样序列大小为1537 bp
屏幕截图 2021-12-13 100957.png
通过比较发现，通过不同位置提取的16S序列都是一样的，当然也可以将contig里面所有的负链区域全部提取出来，然后再对其进行多序列比对进行比较。
将该序列中的所有负链16S序列全部提取出来然后对其进行序列比对
首先将对应的位置提取出来
cat PROKKA_12122021.gff | grep "16S" | cut -f 4,5 | sed -n '1,5p' > negetive_position.txt
接着提取出序列信息
cat negative_position.txt | while read line; do IFS=$'\t'; arr=($line); samtools faidx PROKKA_12122021.fna contig_1:${arr[0]}-${arr[1]}; done > contig_16s_all.fasta
然后进行多序列比对，发现对于来源于同一种的多拷贝序列中的16S序列也不是完全一样的。
屏幕截图 2021-12-14 155123.png
可以参考的链接：https://www.jianshu.com/p/8f4061c38508
http://tiramisutes.github.io/2016/03/18/bedtools.html
https://www.jianshu.com/p/11b8804e570d
主要看：http://tiramisutes.github.io/2016/03/18/bedtools.html
