題:
通過ID從multifasta刪除/刪除序列
andresito
2018-03-20 23:08:44 UTC
view on stackexchange narkive permalink

我有一個這樣的fasta文件：

  >Id1ATCCTT>Id2ATTTTCCC>Id3TTTCCCCAAAA>Id4CCCTTTAAA

我想刪除具有以下ID的序列。 > Id2Id3

ID位於.txt文件中，文本文件將用於匹配和刪除這些序列。

我的輸出應該是這樣的 fasta 文件

  >Id1ATCCTT>Id4CCCTTTAAA

但是我想要 awk 和/或 sed 和/或 bash （沒有python或perl）。
我該怎麼做？

為什麼要限制語言？
確實：使用正確的工具完成工作。除非您可以對輸入文件做出非常嚴格的假設，否則請使用適當的FASTA解析器，而不要使用一些黑手黨的混蛋。
如果在截止日期之前完成工作，則使用您想使用的任何工具。與Unix工具相比，許多解析器的速度都慢得多，但是我承認，它們在理解事物如何工作方面花了很多精力。
五 答案:
Kamil S Jaron
2018-03-21 17:28:48 UTC
view on stackexchange narkive permalink

今天，我編寫了一個腳本來使用 Biopython 精確地執行此操作。如果原始過濾器中沒有要過濾的標題，我還會添加一個警告。因此python腳本 filter_fasta_by_list_of_headers.py 是：

 ＃！/ usr / bin / env python3from Bio import SeqIOimport sysffile = SeqIO.parse（sys.argv [1] ，對於ffile中的seq_record，header_set = set（line.strip（）for open（sys.argv [2]）中的行）：try：header_set.remove（seq_record.name），但KeyError：print（seq_record.format （“ fasta”））如果len（header_set）！= 0，則繼續：print（len（header_set），'在輸入的fasta文件中未標識列表中的頭。'，file = sys.stderr）

和用法：

  filter_fasta_by_list_of_headers.py input.fasta list_of_scf_to_filter >filtered.fasta

PS翻轉腳本以從列表中提取序列非常容易（只是打印行必須在 header_set.remove（seq_record.name ）行之後移動

他們說沒有python，但不知道為什麼。忽略這一點，我的直覺是，假設大多數記錄要​​打印，在這種情況下執行if if header_set：中的seq_record.name會更好。錯誤處理有點貴，LBYL似乎比EAFP快
嗯，我錯過了“ no python”請求。提速技巧聽起來不錯（儘管實際上並不是一段需要被激化的代碼，但它在大約2分鐘的時間內解析了大約十萬個重疊群（同時過濾了一些污染物）？）
+1用於忽略無意義的要求並使用適當的工具。也就是說，我不太確定這種方法的效果。以這種方式使用錯誤在Python中是慣用的，但仍然要付出一定的代價。一個簡單的“ in”測試可能會更有效率。甚至可能沒有必要“刪除”條目，這也可能非常昂貴。
@KonradRudolph好，所以我做了一個建議，建議在正常昆蟲基因組裝配上進行標記-1Gbp，300,000個支架，並給它列出5,000和100,000個支架進行過濾。無論列表的大小或使用的算法是29到30秒。換句話說，如果LBYL比EAFP快，那麼它的規模完全沒有關係。
@KamilSJaron不錯！
啊+1表示我的直覺是錯誤的，另一個想法是，一次寫所有記錄而不是一張一張地寫應該*更快*，這樣快嗎？ https://pastebin.com/n7RUgW5x對不起，如果看起來像微優化，我只是很好奇
Alex Reynolds
2018-03-21 01:03:21 UTC
view on stackexchange narkive permalink

如果您想學習如何使用命令行工具進行操作，則可以使用 awk 使FASTA線性化，通過管道傳輸到 grep 來過濾名為的感興趣項在 patterns.txt 中，然後通過管道傳輸到 tr 進行線性化：

  $ awk'{if（（（NR>1）&&（$ 0〜/ ^ > /））{printf（“ \ n％s”，$ 0）; }否則，如果（NR == 1）{printf（“％s”，$ 0）; } else {printf（“ \ t％s”，$ 0）; }} in.fa | grep -Ff patterns.txt-| tr“ \ t”“ \ n” > out.fa

這將適用於多行FASTA。

優秀的解決方案，除了它將保留patterns.txt中指定的順序，而不是將其刪除
好點子。我一定錯過了那部分問題。可以使用`grep -v`刪除與指定模式匹配的元素。
user172818
2018-03-21 07:09:00 UTC
view on stackexchange narkive permalink

假設您將序列名稱保留在 ids.txt 中，並將序列保留在 seq.fa 中：

  awk'BEGIN {while（（ getline<“ ids.txt”）>0）l [“ >” $ 1] = 1} / ^ > / {f =！l [$ 1]} f'seq.fa
terdon
2018-03-21 00:02:01 UTC
view on stackexchange narkive permalink

使用我之前發布的腳本和 TblToFasta 腳本，您可以執行以下操作：

  FastaToTbl file.fa | grep -vwf ids.txt | TblToFasta > file.2.fa
user1141
2018-03-20 23:41:16 UTC
view on stackexchange narkive permalink

如果sed，awk都假定grep正常。 -A（n）和-B（n）標誌在匹配後給出（n）行。假設您所有的fasta序列都是一行，這有點危險，但是對於您的示例而言，它是可行的。首先刪除要刪除的ID（在rmid.txt中），然後將grep與初始的fasta反向。

  grep -A1 -f rmid.txt fasta.fa > rmfile.fagrep -v- f rmfile.fa fasta.fa

真正的答案是使用一個腳本，該腳本定義除換行符“ \ n”以外的分隔符，然後解析ID，因此最好使用其中之一您不想使用的語言...

這假定每個ID僅包含一行序列，這對於fasta文件而言可能不太正確（但對於fastq通常是正確的）。
使用ip1的問題是子字符串，這也將捕獲id10等。您可以添加--word-regexp標誌，但這通常會確實損害性能。 awk可能是更好的工具
是的，兩個要點，示例仍根據請求完成。 @terdon我將行號指定為問題，\ @Chris_Rands不太可能會形成fasta ID來允許此類問題，但是要點表示讚賞。 FWIW，為了操作fasta，我不會使用grep，awk，sed，因為有很多基於腳本語言的良好軟件包，我更喜歡Bioperl。我是新手，因此無法向OP評論。
啊，你做到了。抱歉，我錯過了免責聲明，只看到了命令。但是，請注意，讓ID為其他ID的子字符串是絕對常見的。例如，標准人類基因組Fasta序列同時具有“ chr1”和“ chr11”。就是說，您的答案還有一個更嚴重的問題：rmfile.fa的每一行都將作為搜索模式傳遞給grep，因此包含該行的_any_行將被刪除。這意味著，如果序列“ id3”為“ TAT”，則將刪除與“ TAT”匹配的_all_行，而不考慮其ID。


ⓘ
該問答將自動從英語翻譯而來。原始內容可在stackexchange上找到，我們感謝它分發的cc by-sa 3.0許可。
Loading...