2020.12.30丨随机读取reads使用ncbi在线工具比对nt数据库

• 背景：最近有个项目，客户让比对到人的参考基因组，比对率却只有70%左右.为了搞清楚测序数据有没有被污染，我们需要随机读取一些reads，放到nt数据库去比对。之前都是一条一条提交，这种批量提取和提交都会遇到一些问题，因此，写这篇文章进行一个统计。
• 比对流程：随机提取序列→fastq转换fasta→提交序列→统计结果
• step.1 随机提取序列
  • 使用工具：seqtk
  • 安装方式：conda install seqtk
  • 使用代码：
    • seqtk sample -s 100 Carisma_control_ATCACGAT.fastq.gz 1000 > 1000reads_control.fastq #这里提取1000条
    • # -s 随机种子，默认值11，不清楚可以百度了解一下
  • 注意：一般进行随机比对会选择5000/10000条，不过ncbi在线平台对进程有控制，我试过，一次性跑不了5000条（PE100）。所以可以尝试跑1000或者2000条。
• step.2 格式转换
  • 这个就比较简单，稍微百度一下，出来的前几个链接都是同样的一行代码可以解决
  • 使用代码：
    • awk '{if(NR%4 == 1){print ">" substr($0, 2)}}{if(NR%4 == 2){print}}' 1000reads_control.fastq > 1000reads_control.fasta
  • 注意：awk是个相当强大的命令，很多时候处理文本都可以用它，建议可以认真学习。
  • 结果展示：
    • 
• step.3 提交序列（在线版）
  • 这个就是一个常规操作了，进入ncbi的Blast比对页面，选择文件进行上次，默认选择nr/nt数据库，直接提交。
  • 网站链接：Nucleotide BLAST: Search nucleotide databases using a nucleotide query (nih.gov)
• step.4 统计结果
  • 分析时间受当地作息影响，国内的话下午提交（美国是半夜）分析速度会比较快。分析结束后会跳转到分析界面。在Download All处可以下载到所有序列比对出来的结果，选择text可以看到比对详细结果。注意，建议下载数据前先在下方表格处调整展示条目，这样可以缩小文件大小，方便查阅。
  • 
    • 
  • 结果展示
    • 这里以text文本为例，可以看到序列大部分比对到载体上了，怪不得比对到人的参考基因组时比对率会比较低。
      • 
• 参考文章：
  • 从 fastq 文件中随机抽取 reads - 知乎 (zhihu.com)
  • FASTQ 文件格式转换为 FASTA 格式 - 遗世独立的愚公 - 博客园 (cnblogs.com)
•