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)