使用 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":""),[=15=]);N++;next;} {printf("%s",[=15=]);} 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[[=10=]]
  next
}
/^>/{
  found=""
  match([=10=],/^>[^.]*/)
  if(substr([=10=],RSTART+1,RLENGTH-1) in arr){
    found=1
  }
}
found
' transcript_id.txt myfile.fasta

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

awk '
FNR==NR{
  arr[[=11=]]
  next
}
/^>/{
  found=""
  if( 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[[=12=]]                          ##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([=12=],/^>[^.]*/)              ##using match function to match regex starting from > to before dot occured in current line.
  if(substr([=12=],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[[=13=]]                     ##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( 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[];next}
   substr(,1,index(,".")-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[[=11=]]; next }
    sub(/^>/, "") &&  in a
' transcript_id.txt myfile.fasta \
\
> transcript_ids.txt

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

Extract FASTA sequences (with version number) using sequence IDs (without version number) listed in txt file

unix

awk

grep

extract

fasta