时间戳

首页 > 解决方案 > 使用 txt 文件中列出的序列 ID（不带版本号）提取 FASTA 序列（带版本号）

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

问题描述

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

我的文件.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

成绩单_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

标签： unixawkgrepextractfasta

解决方案

第一个解决方案：您能否尝试以下操作。用 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.

推荐阅读