- 背景:最近有个项目,客户让比对到人的参考基因组,比对率却只有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文本为例,可以看到序列大部分比对到载体上了,怪不得比对到人的参考基因组时比对率会比较低。
- 这里以text文本为例,可以看到序列大部分比对到载体上了,怪不得比对到人的参考基因组时比对率会比较低。
- 参考文章:
- 从 fastq 文件中随机抽取 reads - 知乎 (zhihu.com)
- FASTQ 文件格式转换为 FASTA 格式 - 遗世独立的愚公 - 博客园 (cnblogs.com)