组装细菌基因组

一、找到合适的基因序列

1.在Genome Announcements网站找一篇细菌基因组文章，并找到文章记载的SRA号；

• 找了一个一篇才发布不久的关于Mycobacterium orygis Strain全基因组测序的文章https://mra.asm.org/content/8/40/e01080-19

image.png

• 在文章中可以看到文章RUN的数据为在SRR9157804

image.png

• 点击进去会跳转到ncbi的SRA数据库，可以看到具体的Run的信息，只有300M左右，数据量不大（https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR9157804）

image.png

二、下载sra序列

• 使用下列命令下载下来

prefetch SRR9157804

• 成功下载，下载后的内容放在~/ncbi/public/sra路径下

  
  image.png

三、Fastq-dump解压

1.新建一个文件夹存放该实验所有操作的结果

 mkdir ~/ncbi/public/sra/xijun
 mv  SRR9157804.sra  ~/ncbi/public/sra/xijun/
 cd  ~/ncbi/public/sra/xijun/

2.解压SRA文件为fastq格式

• 使用下列命令解压

fastq-dump --gzip --split-files  SRR9157804.sra

• 成功转换成以fastq.gz结尾的两个文件，因为是双端测序，有正向和反向两个文件

  
  image.png

三、Fastqc质控

• 输入下列命令

fastqc SRR9157804_1.fastq.gz
fastqc SRR9157804_2.fastq.gz

• 得到下列结果

  
  image.png

• 网页上查看一下序列质量

  
  image.png image.png
  

  可以看出这个样本正反序列都质量很不错，报告中只有一个被警告，其他的都很成功

四、数据过滤

• 输入以下命令

 mkdir trim_out #用来存放过滤后的文件
 java -jar ~/Biosofts/Trimmomatic-0.38/Trimmomatic-0.38/trimmomatic-0.38.jar PE -phred33 SRR9157804_1.fastq.gz  SRR9157804_2.fastq.gz ./trim_out/output_forward_paired.fq.gz ./trim_out/output_forward_unpaired.fq.gz ./trim_out/output_reverse_paired.fq.gz ./trim_out/output_reverse_unpaired.fq.gz ILLUMINACLIP:/home/zhouqi/Biosofts/Trimmomatic-0.38/Trimmomatic-0.38/adapters/TruSeq2-PE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:20 TRAILING:20 MINLEN:60

• 文章中过滤参数为minimum length, >60 bp
• Trimmomatic 过滤数据的步骤与命令行中过滤参数的顺序有关，通常的过滤步骤如下：
  ILLUMINACLIP: 过滤 reads 中的 Illumina 测序接头和引物序列，并决定是否去除反向互补的 R1/R2 中的 R2。
  SLIDINGWINDOW: 从 reads 的 5' 端开始，进行滑窗质量过滤，切掉碱基质量平均值低于阈值的滑窗。
  MAXINFO: 一个自动调整的过滤选项，在保证 reads 长度的情况下尽量降低测序错误率，最大化 reads 的使用价值。
  LEADING: 从 reads 的开头切除质量值低于阈值的碱基。
  TRAILING: 从 reads 的末尾开始切除质量值低于阈值的碱基。
  CROP: 从 reads 的末尾切掉部分碱基使得 reads 达到指定长度。
  HEADCROP: 从 reads 的开头切掉指定数量的碱基。
  MINLEN: 如果经过剪切后 reads 的长度低于阈值则丢弃这条 reads。
  AVGQUAL: 如果 reads 的平均碱基质量值低于阈值则丢弃这条 reads。
  TOPHRED33: 将 reads 的碱基质量值体系转为 phred-33。
  TOPHRED64: 将 reads 的碱基质量值体系转为 phred-64。
  <链接：https://www.jianshu.com/p/7a248999063f>

• 成功运行

  
  image.png

• 生成以下四个文件

  
  image.png
• 对正反匹配上的序列用fastqc看下质量变化如何

fastqc output_forward_paired.fq.gz output_reverse_paired.fq.gz

• 在网页上查看一下过滤后结果

  
  image.png image.png
  

  可以看出过滤掉了318024个碱基

五、组装基因组草图

只看正反配对上的，不考虑未配对上的
可以采用Spades和velvet组装，我这里使用后者

1.velveth利用数据构建一个hash表

• 输入下列命令

velveth velvet_out 31 -shortPaired -fastq -separate output_forward_paired.fq.gz output_reverse_paired.fq.gz  

31代表kmer值

• 成功创建

  
  image.png

2.velvetg进行序列拼接

• 输入下列命令

 velvetg velvet_out -exp_cov auto -cov_cutoff auto -very_clean yes

• 生成的contigs.fa ，为我们最终需要的拼接结果文件

  
  image.png

七、Quast评价组装的基因组效果

• 输入下列命令

quast.py  contigs.fa  -o  ./quast_out

• 生成了一个report.html，可以下载下来查看拼接结果

  
  image.png

• 拼接结果如下

  
  image.png image.png

可以看出拼接后contigs数有146个，序列长度为4217673bp,N50为67893，GC含量为65.41%，这些数据和文章中有点差距，可能是文章所用的是spades组装的，它用的K值是自动选择的，我用的velvet，且K值设置为31

image.png

• 我把K值调高到77后，发现和文章中的数据没有很大差别

  
  image.png