當前位置: 華文星空 > 知識

轉錄組比對步驟中的指令碼執行問題與echo進行debug方法

2021-03-16知識

學徒在完成轉錄組實戰任務的時候,出現了一個小bug,自顧自的糾結了一晚上不過來尋求我的幫助。讓我非常生氣,因為確實是很簡單的兩個小技巧,主要是理解全域變量、環境變量。

就是hisat2進行批次比對程式碼失效

範例:當前目錄檔如下:

1.2G 10月 30 10:14 SRR1039510_1_val_1.fq.gz 1.2G 10月 30 10:14 SRR1039510_2_val_2.fq.gz 1.1G 10月 30 10:14 SRR1039511_1_val_1.fq.gz 1.1G 10月 30 10:14 SRR1039511_2_val_2.fq.gz 1.5G 10月 30 10:18 SRR1039512_1_val_1.fq.gz 1.5G 10月 30 10:18 SRR1039512_2_val_2.fq.gz 816M 10月 30 10:12 SRR1039513_1_val_1.fq.gz 820M 10月 30 10:12 SRR1039513_2_val_2.fq.gz 2.1G 10月 30 10:27 SRR1039514_1_val_1.fq.gz 2.1G 10月 30 10:27 SRR1039514_2_val_2.fq.gz 253M 10月 30 10:05 SRR1039515_1_val_1.fq.gz 253M 10月 30 10:05 SRR1039515_2_val_2.fq.gz 131M 10月 30 10:05 SRR1039516_1_val_1.fq.gz 135M 10月 30 10:05 SRR1039516_2_val_2.fq.gz 1.8G 10月 30 10:23 SRR1039517_1_val_1.fq.gz 1.8G 10月 30 10:23 SRR1039517_2_val_2.fq.gz 1.6G 10月 30 10:20 SRR1039518_1_val_1.fq.gz 1.6G 10月 30 10:20 SRR1039518_2_val_2.fq.gz 451M 10月 30 10:07 SRR1039519_1_val_1.fq.gz 458M 10月 30 10:07 SRR1039519_2_val_2.fq.gz 645M 10月 30 10:09 SRR1039520_1_val_1.fq.gz 648M 10月 30 10:09 SRR1039520_2_val_2.fq.gz 856M 10月 30 10:12 SRR1039521_1_val_1.fq.gz 859M 10月 30 10:12 SRR1039521_2_val_2.fq.gz 1.6G 10月 30 10:20 SRR1039523_1_val_1.fq.gz 1.6G 10月 30 10:20 SRR1039523_2_val_2.fq.gz

他首先定義index變量,為ref檔路徑

index=/refdir/server/reference/index/hisat/hg38/genome

然後執行指令碼進行hisat2比對,程式碼如下

for id in {08..23} do hisat2 -t -x $index \ -1 /home/project/rna/clean/SRR10395${id}_1_val_1.fq.gz \ -2 /home/project/rna/clean/SRR10395${id}_2_val_2.fq.gz \ -S /home/project/rna/align/SRR10395${id}.sam done

送出指令碼後出現Exit,也就是說報錯了!

  1. 此時如何進行故障排查?

迴圈中加echo,檢視指令碼中所執行的真實程式碼

for id in {08..23} do echo hisat2 -t -x $index \ -1 /home/data/jjl/project/rna/clean/SRR10395${id}_1_val_1.fq.gz \ -2 /home/data/jjl/project/rna/clean/SRR10395${id}_2_val_2.fq.gz \ -S /home/data/jjl/project/rna/align/SRR10395${id}.sam done

執行後真實程式碼打印出來如下

"hisat2 -t -x $index \ -1 /home/data/jjl/project/rna/clean/SRR1039508_1_val_1.fq.gz \ -2 /home/data/jjl/project/rna/clean/SRR1039508_2_val_2.fq.gz \ -S /home/data/jjl/project/rna/align/SRR1039508.sam"

發現$id已經成功呼叫08,09,10.......

但是 $index 沒有成功呼叫,echo出來還是 $index 本身,也就是說這個變量並沒有被解析啊!

故障原因:在指令碼中執行的程式碼如果要呼叫變量index, 指令碼外 的環境定義的變量index是 無法被成功呼叫的

簡言之: "全域" 定義的變量不能有效覆蓋 "局部" 指令碼內變量的呼叫需求!

  1. 應該如何解決故障?

必須在 指令碼內 現場定義變量進行賦值index,然後 指令碼內 才能成功呼叫index。

因此調整 指令碼全部的程式碼 如下(即添加下面第一行)

# 添加這一行 index=/refdir/server/reference/index/hisat/hg38/genome # 後面程式碼不變 for id in {08..23} do echo hisat2 -t -x $index \ -1 /home/project/rna/clean/SRR10395${id}_1_val_1.fq.gz \ -2 /home/project/rna/clean/SRR10395${id}_2_val_2.fq.gz \ -S /home/project/rna/align/SRR10395${id}.sam done

執行成功!

兩個小tips

  • 寫好的指令碼需要多echo看看。
  • 先執行測試數據看自己的指令碼是否正確,然後上真實的大量樣品數據專案。