使用 txt 文件中列出的序列 ID（不带版本号）提取 FASTA 序列（带版本号）

Question

我想根据 transcript_id.txt 文件中列出的 id 从myfile.fasta 中提取特定序列。 我的主要问题是我的 transcript_id.txt 文件只列出了成绩单 id，而 fasta 文件也有成绩单版本，并且 transcript_id.txt 中列出的成绩单在 fasta 文件中可以有多个版本。 我尝试了几种方法（如下所列），但无法得到我需要的。

myfile.fasta

>transcriptA.t1
ATGGAAGTAGGAAGTGGTTTGGGGAATGATGAACGACGCACCTGGCCAGCCGATGTTCTA
GCAGCAGCAGATGCCTACGTTGTCCTCGCTGCAGAGTACAACCACAGCCTTCCCCCGGCA
AATCAAAGCATCCTAAATGAAGCATGTGTGACAGCAGGATCATCTGCCAGGGTTTCCAAA
CTCACCAACCTAATGGACCACTTCTTCTCGTACAAGTATCGTCCATCTGGCATTGTCTGC
TTGGCTGTTGAGACTTATGTCTGCTGGGGAGCTGAAGCACGCCTACATCTTTCACGCCTT
>transcriptA.t2
GTGCCAAATTGGGATCTGGAGAAACCTGCAGCTTTTGATCTTACTGTTGCATCACTGCTT
TTCATCAAGAAG
>transcriptB.t1
GGAGTGGTGAGTCTTCTGCTCACTGCACTAAACCTTAATGACACTGGGACCTACAGCTGC
GTTGATGTGCCTCAATCGATTGATGAATTTGCTCGGCGTCATCCTCGGCGATTGCAATTA
GTAGATATTCTTAACGATTGA
>transcriptC.t1
AAAAGGCTCTGGGAGTTTGAGGCCAACGGGGGGGGAGGCCCCATTACCTCCAGCGTCAGG
TTCGCTGATGTGTACAACGATGGAACCCTAGACGTGATCTTTGTCACCCTCACCGGAACC
TTCTGGGTCCTGGAGGGGCTCACTGGA
>transcriptD.t1
CAAAGGAAGCATGCCTCTAATGATGCTAAGTGCTCAGAGCTAGGTTGGTCATGCATACCA
GCCATGGGAGACCCGTCCCCATCCATCCAATGGAGCTCCACTCCGCAACTTTAG
>transcriptD.t2
ATGCCTCAGGTGAATGTGGCCCCACCCACTGCCAAGGTCAAGGGGGCGTGTAGGGTTGTA
TGTGGTGCGTTGTCTCTCCAACAATTCATTATGCCCGACCAAGAGGTGTCACCTGTTCAG
CAAGGAGAATCTGACCATTTGCACATTGAAGCTTTCACTCTGGTGTCTGGAGGCACGAGT
ACGGATGTCGTAACTTTGCTGCAAGAGCAATACAAAGTGCTAAGCTGA
>transcriptE.t1
TCTATTCCAGTAGTCTACAAGGCACTGACCCCTGAGGGAGTGCCATCCAACAAGGAAGAT
GTCATTGGAGTGGTGCCAGACAAGGTTGTTGGAACACCAGATGTGAAGCCAACTGGAAGT
GTAGCTGCTGCTGCTGCCTTGGGCGTGTGCAAAGTGACTCT

transcript_id.txt

transcriptA
transcriptC
transcriptD

目标

>transcriptA.t1
ATGGAAGTAGGAAGTGGTTTGGGGAATGATGAACGACGCACCTGGCCAGCCGATGTTCTA
GCAGCAGCAGATGCCTACGTTGTCCTCGCTGCAGAGTACAACCACAGCCTTCCCCCGGCA
AATCAAAGCATCCTAAATGAAGCATGTGTGACAGCAGGATCATCTGCCAGGGTTTCCAAA
CTCACCAACCTAATGGACCACTTCTTCTCGTACAAGTATCGTCCATCTGGCATTGTCTGC
TTGGCTGTTGAGACTTATGTCTGCTGGGGAGCTGAAGCACGCCTACATCTTTCACGCCTT
>transcriptA.t2
GTGCCAAATTGGGATCTGGAGAAACCTGCAGCTTTTGATCTTACTGTTGCATCACTGCTT
TTCATCAAGAAG
>transcriptC.t1
AAAAGGCTCTGGGAGTTTGAGGCCAACGGGGGGGGAGGCCCCATTACCTCCAGCGTCAGG
TTCGCTGATGTGTACAACGATGGAACCCTAGACGTGATCTTTGTCACCCTCACCGGAACC
TTCTGGGTCCTGGAGGGGCTCACTGGA
>transcriptD.t1
CAAAGGAAGCATGCCTCTAATGATGCTAAGTGCTCAGAGCTAGGTTGGTCATGCATACCA
GCCATGGGAGACCCGTCCCCATCCATCCAATGGAGCTCCACTCCGCAACTTTAG
>transcriptD.t2
ATGCCTCAGGTGAATGTGGCCCCACCCACTGCCAAGGTCAAGGGGGCGTGTAGGGTTGTA
TGTGGTGCGTTGTCTCTCCAACAATTCATTATGCCCGACCAAGAGGTGTCACCTGTTCAG
CAAGGAGAATCTGACCATTTGCACATTGAAGCTTTCACTCTGGTGTCTGGAGGCACGAGT
ACGGATGTCGTAACTTTGCTGCAAGAGCAATACAAAGTGCTAAGCTGA

试过了：

https://bioinformaticsworkbook.org/dataWrangling/fastaq-manipulations/retrieve-fasta-sequences-using-sequence-ids.html#gsc.tab=0

1) ncbi-blast/makeblastdb

makeblastdb -in myfile.fasta -parse_seqids -dbtype nucl
blastdbcmd -entry "lcl|transcriptA.t1" -db myfile.fasta -target_only

这部分工作，但我只能通过输入准确的转录版本并添加“lcl|”来搜索它。 没有设法使用通配符或 transcript_id.txt。

https://www.biostars.org/p/319099/

2) grep

grep -w -A 2 -f transcript_id.txt myfile.fasta --no-group-separator

这很好用，但我必须将 -A 设置为一些数量，并且每个成绩单的行数各不相同。

3) 线性化 fasta 文件

awk '/^>/ {printf("%s%s\t",(N>0?"\n":""),$0);N++;next;} {printf("%s",$0);} END {printf("\n");}' < myfile.fasta > linear.fasta
awk '/^>transcriptA/' linear.fasta

再次，我不知道如何使用 awk 逐一搜索带有 transcript_id.txt 的线性化 fasta 文件。

4) seqkit

这仅在我 *将脚本版本添加到 transcript_id.txt 时才有效。任何使用 --id-regexp 的尝试都失败了。

seqkit grep -n -f transcript_id*.txt myfile.fasta

Answer 1

第一种解决方案：您能否尝试以下操作。使用 GNU awk 中的示例编写和测试。

awk '
FNR==NR{
  arr[$0]
  next
}
/^>/{
  found=""
  match($0,/^>[^.]*/)
  if(substr($0,RSTART+1,RLENGTH-1) in arr){
    found=1
  }
}
found
' transcript_id.txt myfile.fasta

第二个解决方案： 使用多个字段分隔符选项。

awk '
FNR==NR{
  arr[$0]
  next
}
/^>/{
  found=""
  if($2 in arr){ found=1 }
}
found
' transcript_id.txt FS="[>.]" myfile.fasta

---

---

第一种解决方案的解释：

awk '                              ##Starting awk program from here.
FNR==NR{                           ##Checking condition FNR==NR which will be TRUE when transcript_id.txt file is being read.
  arr[$0]                          ##Creating arr with index of current line.
  next                             ##next will skip all further statements from here.
}
/^>/{                              ##Checking condition if line starts from > then do following.
  found=""                         ##Nullifying found here.
  match($0,/^>[^.]*/)              ##using match function to match regex starting from > to before dot occured in current line.
  if(substr($0,RSTART+1,RLENGTH-1) in arr){ ##Checking condition if sub string which we get after above succesul matched regex is present in arr
    found=1                        ##Then setting found to 1 here.
  }
}
found                              ##Checking condition if found is SET then print that line.
' transcript_id.txt myfile.fasta   ##Mentioning Input_file names here.

第二种解决方案的解释：

awk '                         ##Starting awk program from here.
FNR==NR{                      ##Checking condition FNR==NR which will be TRUE when transcript_id.txt is being read.
  arr[$0]                     ##Creating arr with index of current line here.
  next                        ##next will skip all further statements from here.
}
/^>/{                         ##Checking condition if line starts from > then do following.
  found=""                    ##Nulifying found here.
  if($2 in arr){ found=1 }    ##Checking condition if 2nd field is present in arr then set found.
}
found                         ##Checking condition if found is set then print line.
' transcript_id.txt FS="[>.]" myfile.fasta  ##Mentioning Input_file names and setting FS as > or dot for Input_file myfile.fasta here.

Answer 2

awk '
   FNR==NR {a[$1];next}
   substr($1,1,index($1,".")-1) in a
' transcript_id.txt RS='>' FS='\n' ORS='>' myfile.fasta

Answer 3

部分问题在于文本文件中列出的成绩单 ID 不包含版本号。如果您的成绩单 ID 文件确实包含带有版本号的成绩单 ID，则您可以为此使用 samtools faidx。使用最新版本的 samtools（当前为 v1.11），您只需要：

samtools faidx -r transcript_ids.txt myfile.fasta

那么实际的问题就变成了：如何获取带有版本号的成绩单ID列表？您可以为此使用任何文本处理工具。例如，使用awk:

awk -F "." '
    FNR==NR { a[$0]; next }
    sub(/^>/, "") && $1 in a
' transcript_id.txt myfile.fasta \
\
> transcript_ids.txt