shell實現鹼基/氨基酸序列提取（sed,grep,awk三大利器走向天下）

假定存在file.fa文件，其序列如下所示：

現在有四個任務，分別如下所示：

任務一：提取包含“RBFOX1”字符的序列

第一步：sed -i '/>/i\@' file.fa

file.fa指的是包含所有序列的文件；這一步的意思是在以">"開頭的字符串前一行加上"@",方便后續序列的識別和匹配;

修改后的文件效果如下所示：

第二步：sed -n '/RBFOX1 /,/^@/p' file.fa > RBFOX1.fa

匹配“RBFOX1”和"@"字符的行，並打印包含“RBFOX1”和"@"字符之間的行

效果如下：

第三步：sed -i '/^@/d' RBFOX1.fa

刪掉包含"@"所在的行

以上是提取RBFOX1序列的工作。

從上面可以看出，提取出來的RBFOX1序列不止一條，但很多時候我們只需要其中最長的一條，由此，引申了任務二。

任務二：從多條“RBFOX1”序列中提取最長的一段序列

第一步：need1=`sed -n "/RBFOX1 /,/^@/p" file.fa | grep -En "@|>" | awk -F ":" 'NR==1{tmp=$1}NR>=0{print $1-tmp,$2;tmp=$1}'`

計算不同序列的長度，這里統計序列的行數(下圖紅色方框)，行數越多，表明序列越長；

效果如下：

第二步：need2=`echo "$need1" | awk 'BEGIN {max = 0} {if ($1+0 > max+0) max=$1} END {print max}'`

統計最長的序列；

結果如下：

第三步：need3=`echo "$need1" | grep -B 1 "^$need2" | head -1 | awk '{print $2}'`

輸出最長序列的ID；

結果如下：

第四步：提取最長序列

sed -n  "/^$need3/,/^@/p" file.fa > RBFOX1.fa

sed -i '/^@/d' RBFOX1.fa

結果如下：

以上四步，即可完成從多條“RBFOX1”序列中提取最長序列的工作。

但在實際工作中，我們需要提取多個基因的序列，而非單個基因（比如只提取RBFOX1基因）的序列，由此引出任務三。

任務三：提取多個基因，且每條轉錄本的序列最長

首先，准備包含多個基因的文件genename,一個基因一行，其內容如下所示：

隨后，輸入如下命令：

while read line;do
need1=`sed -n "/${line} /,/^@/p" file.fa | grep -En "@|>" | awk -F ":" 'NR==1{tmp=$1}NR>=0{print $1-tmp,$2;tmp=$1}'`
need2=`echo "$need1" | awk 'BEGIN {max = 0} {if ($1+0 > max+0) max=$1} END {print max}'`
need3=`echo "$need1" | grep -B 1 "^$need2"  | head -1 | awk '{print $2}'`
sed -n  "/^$need3/,/^@/p" file.fa > ${line}.fa
sed -i '/^@/d' ${line}.fa
done < genename

以上命令即可一次性提取多個基因。

除了一次性提取多個基因外，我們還想一次性提取多個物種的多個基因，這個時候依舊可以使用循環命令完成工作。

任務四：提取多個物種、多個基因，且每條轉錄本的序列最長

假定有5個物種的序列存儲在file.fa,file1.fa,file2.fa,file3.fa,file4.fa中，每一個fa文件代表一個物種。

此外，通過命令cat *.fa >> all.fa 生成 all.fa 文件,該文件包含5個物種的所有序列。

按照慣例，先在每個fa文件中加入"@",方便后續序列的識別和匹配，命令如下：

for j in *.fa; do
echo $j
sed -i '/>/i\@' $j 
done

准備好以上文件后，輸入如下命令即可完成多物種多基因，且每條轉錄本的序列最長的提取工作：

while read line;do
       for i in file*.fa;do
                 need1=`sed -n "/${line} /,/^@/p" "$i" | grep -En "@|>" | awk -F ":" 'NR==1{tmp=$1}NR>=0{print $1-tmp,$2;tmp=$1} '`
                 need2=`echo "$need1" | awk 'BEGIN {max = 0} {if ($1+0 > max+0) max=$1} END {print max}'`
                 echo "$need1" | grep -B 1 "^$need2"  | head -1 | awk '{print $2}' >> ${line}name
       done

sed -i '/^\s*$/d' ${line}name

      while read line1;do
                 echo ${line1}
                 sed -n -e "/^${line1}/,/^@/p"  all.fa >>  ${line}.fa 
      done < ${line}name

sed -i '/^@/d' ${line}.fa

done < genename

---

此推文感謝PCC同學的建議~

---